数字信号处理的matlab实现

数字信号处理的matlab实现简谐振动的特性完全取决于振幅、频率、初相位角。

x1 = (0.5).^t;x1 = 0.5* sin(2*pi*f*t+pi/4);x1 = [(n - n0) >= 0]; %阶跃信号x1 = [(n-n0) == 0]; %脉冲信号y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除n = [-10:10];alpha = -0.1 + 0.3 * j;x = exp(alpha * n);%复指数信号rx = real(x);ix = imag(x);mx = abs(x); %振幅px = (180/pi) * angle(x); % 相位，并转换为度x = rand(1,10);x = randn(1,10); %guass series numbersxp = [x x x x x];x2 = fliplr(x);%左右折叠第二章信号%test sampling ruledt = 0.01; %samping frequence for draw 100Hzt = 0:dt:1;f = 10;x = sin(2 * pi * f * t + 0.3);dt1 = 0.1; t1 = 0:dt1:1;%10Hzx1 = sin(2 * pi * f * t1 + 0.3);dt2 = 0.05; t2 = 0:dt2:1;%10Hzx2 = sin(2 * pi * f * t2 + 0.3);subplot(3,1,1), plot(t, x);title('f = 10Hz, fs = 100Hz');subplot(3,1,2), plot(t1, x1), hold on, stem(t1, x1), plot(t, x);title('f = 10Hz, fs = 10Hz');subplot(3,1,3), plot(t2, x2), hold on, stem(t2, x2), plot(t, x);title('f = 10Hz, fs = 20Hz');基本信号：x1 = (0.5).^t; %指数信号x1 = 0.5* sin(2*pi*f*t+pi/4); %余弦信号x1 = [(n - n0) >= 0]; %阶跃信号x1 = [(n-n0) == 0]; %脉冲信号y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除，sinc函数信号n = [-10:10];alpha = -0.1 + 0.3 * j;x = exp(alpha * n);%复指数信号rx = real(x);ix = imag(x);mx = abs(x); %振幅px = (180/pi) * angle(x); % 相位，并转换为度x = rand(1,10); %随机信号x = randn(1,10); %高斯随机序列xp = [x x x x x]; %复制信号运算：时移，倒置，尺度改变x2 = fliplr(x);%左右折叠信号加信号微分和积分% test2 微分和积分clear all; clc;dt = 0.01; t = 0:dt:4*pi;y = sin(t);y1 = diff(y)/dt; % 微分，dt为采样间隔for ii = 1:length(y) %积分，dt为采样间隔y2(ii) = sum(y(1:ii)) * dt;endsubplot(3,1,1), plot(t, y);subplot(3,1,2), plot(t, [0 y1]); %微分后信号比原信号少一个元素，用零补subplot(3,1,3), plot(t,y2);grid on;xlabel('Time/s');ylabel('Amplitude');信号乘注意：若参与信号乘的两个信号长度不一样，则必须进行转换之后才能在matlab中进行操作。

% test3 信号乘，滤波器滤波。

注意振幅谱的绘制方法clear all; clc;dt = 0.02; df = 1/(6000 * dt);% 采样间隔为0.02s，则2min内数据个数为6000，df为信号频率分辨率n = 0:2999; %折叠频率之前（取前3000个）数据进行操作f = n * df; %给出频率序列sig = rand(1, length(n)); % 运用随机函数产生信号振幅谱filt = [ones(1, 5/df), zeros(1, (length(n) - 5/df))]; % 理想滤波器幅频响应函数subplot(3,1,1), plot(n * df, sig);subplot(3,1,2), plot(n * df, filt, 'LineWidth', 3);subplot(3,1,3), plot(n * df, sig .* filt);grid on;xlabel('Time/s');ylabel('Amplitude');信号奇偶性第三章Fourier变换一个带相位的余弦函数可以看成一个不带相位的正弦函数与一个不带相位的余弦函数的合成。

谐波函数是周期函数中最简单的函数，它描述的也是最简单的周期现象，在实际中遇到的周期现象往往比它复杂的多。

但这些复杂函数均在一定近似程度上可分解为不同频率的正弦函数和余弦函数。

可以将一种复杂的函数分解为一系列不同频率的正弦函数和余弦函数。

1 周期函数的Fourier变换连续傅里叶级数2 离散Fourier变换离散傅里叶级数对于一个有限长观测数据序列来说，使用Fourier级数法求得的各种周期或频率都是有限的。

????=1????=???????=，T为所取数据总的时间长度????Nyquist频率由取样间隔???来决定，为取样间隔2倍的倒数。

Ck表示了k次谐波的振幅大小。

观察哪些分量的振幅大，哪些分量的振幅小，看出所给波列种含有哪种频率成分，如果我们想滤掉某种频率的波，直接去掉该频率的谐波即可。

我们通常用频率作为横坐标，用振幅作为纵坐标绘出图形来研究频率与振幅的关系，简称为振幅谱。

3 信号的Fourier分解与合成举例% test3 DFTclear allN = 256; dt = 0.05; % data numbers and sampling intervel,sampling frequence is 20Hzn = 0:N-1; t = n * dt; % 序号序列和时间序列x1 = sin(2 * pi * t);x2 = 0.5 * sin(2 * pi * 5 * t);x = sin(2 * pi * t) + 0.5 * sin(2 * pi * 5 * t); %signals addm = floor(N/2) + 1; %down for integera = zeros(1, m);b = zeros(1, m);for k = 0:m-1for ii = 0:N-1a(k+1) = a(k+1) + 2/N * x(ii+1) * cos(2 * pi * k * ii/N);%matlab's array index must be increase from 1 b(k+1) = b(k+1) + 2/N * x(ii+1) * sin(2 * pi * k * ii/N);endc(k+1) = sqrt(a(k+1).^2 + b(k+1).^2);endif (mod(N, 2) ~= 1) a(m) = a(m)/2;endfor ii = 0:N-1xx(ii + 1) =a(1)/2;for k = 1:m-1;xx(ii + 1) = xx(ii + 1) + a(k+1) * cos(2 * pi * k * ii/N) + b(k+1) * sin(2 * pi * k * ii/N);endendsubplot(4,1,1), plot(t, x1); title('orignal signal One '),xlabel('Time/s');subplot(4,1,2), plot(t, x2); title('orignal signal Two'),xlabel('Time/s');%subplot(4,1,1), plot((0:N-1) * dt, xx);title('Composedsitgnal');subplot(4,1,3), plot(t, x); title('orignal signal'),xlabel('Time/s');subplot(4,1,4), plot((0:m-1)/(N * dt), c), title('Fouriertransform'),xlabel('Frequence/Hz'),ylabel('Amplitude'); 此处的1Hz和5Hz的振幅与原来信号振幅不完全一致，是由于数据采样点较少导致的，即N较小。

离散有限信号的频谱为周期谱。

方波及其Fourier分解系数和合成与原始方波信号对比% test4 方波的DFT分解和合成clear all; close allN = 200; dt = 4/N;for n = 1:N %方波信号if (n * dt >= 2)x(n) = 0.8;elsex(n) = -0.8;endendfigure(1)subplot(2,1,1), plot((1:N) * dt, x), hold on;plot((1:N) * dt, zeros(1, N), 'k'), xlabel('Time/s'), title('Orignal signal');a = zeros(1, N),b = zeros(1, N);nn = floor(N/2) + 1;for k = 0:nn - 1%DFT分解方波信号，得到a和ba(k+1) = 0;b(k+1) = 0;for ii = 0:N-1a(k+1) = a(k+1) + 2/N * x(ii + 1) * cos(2 * pi * k * ii /N);b(k+1) = b(k+1) + 2/N * x(ii + 1) * sin(2 * pi * k * ii /N);endc(k+1) = sqrt(a(k+1).^2 + b(k+1).^2);endsubplot(2,1,2), plot((0:nn - 1) / (N * dt), c);title ('DFT'), xlabel('Frequence/Hz'),ylabel('Amplitude');m = input('Input the max k value: ');%观察m的值不到N/2的时候的合成情况if (m > floor(N/2) + 1)error('Too large!');endif (mod(N, 2) ~= 1) a(nn) = a(nn)/2;endfor ii = 0:N-1xx(ii+1) = a(1)/2;for k = 1:mxx(ii+1) = x(ii+1) + a(k+1) * cos(2 * pi * k * ii /N) + b(k+1) * sin(2 * pi * k * ii /N);%合成信号endendfigure(2)plot((1:N) * dt,xx, (0:N-1) * dt, x)hold onplot((1:N) * dt, zeros(1, N), 'k'), xlabel('Time/s'), title('composed signal');最基本的Fourier变换的形式。