1、Harbin Institute of Technology大作业一课程名称： 试验方法与数字信号处理 院 系： 机械电子 班 级： 15S0825 学 号： 姓 名： 哈尔滨工业大学给出信号 ()=sin(210)+sin(280)+ sin(2200) 1. 绘出信号波形。利用 matla 软件，绘制出的原信号波形如图 1 所示。图 1 原波形信号2. 低通滤波，分别用 FIR，IIR 滤波器，保留 10Hz，去除 80Hz 和 200Hz，并画出波形，并与 10Hz 信号对比。解：原信号的最大 Fmax = 200Hz， 取：=103 12= 1400=0.0025此时，满足采样定理。（
2、1 ） 、用 FIR 滤波器（附录 1）选择低通滤波的截止频率为 50Hz，滤波器项数为 80，通过 FIR 滤波器公式，可得到滤波后的信号。编写 matlab 程序，对比滤波后信号和 10Hz 信号，如图 2 所示。图 2 FIR 滤波后信号与 10Hz 信号对比通过图 2 可以发现，滤波后的信号大致反应了 10Hz 信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了 80Hz，200Hz 的信号。为了进一步说明问题，绘制滤波后信号的频谱图，如图 3 所示。从图 3 可以看出，随着 N的增大，10Hz 信号幅值衰减的程度变小，会趋于至原幅值的一半，其余信号幅值衰减的程度变大，滤波
3、效果更加明显。图 3 FIR 滤波后频谱（N = 8，30， 80， 800） 10Hz尝试用汉宁窗口对泄漏进行修正，修正前后的波形如图 4 所示。图 4 采用汉宁窗口修正（2 ） 、用 IIR 滤波器（附录 2）选择低通滤波的截止频率为 50Hz 的二阶 IIR 滤波器，根据相关公式，可以得到 IIR 滤波器的滤波因子，进而可得到滤波后的信号。编写 matlab 程序，对比滤波后信号和 10Hz 信号，如图 5 所示。图 5 IIR 滤波后信号与 10Hz 信号对比通过图 5 可以发现，滤波后的信号大致反应了 10Hz 信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了 80Hz
4、，200Hz 的信号。在滤波信号开始阶段，会出现一较大的波动，该波动会随滤波的进行而消失。为了便于说明问题，绘制出滤波后信号的频谱，如图 6 所示。从图 6 可以看出，滤波后的信号幅值基本与原幅值一样，且高频信号衰减幅度比较大，滤波效果比 FIR 滤波效果好。图 6 IIR 滤波后频谱3、带通滤波，分别用 FIR，IIR 滤波器，保留 80Hz，去除 10Hz 和 200Hz，并画出波形，并与 10Hz 信号对比。解：原信号的最大 Fmax = 200Hz， 取：=103 12= 1400=0.0025满足采样定理。（1 ） 、用 FIR 滤波器（附录 3）选择带通频率为 40120Hz， 即
5、 F1 = 40Hz， F2 = 120Hz，滤波器项数为 80，根据公式，可得相应的滤波因子，编写相应的程序，可得到滤波后的信号，如图 7 所示。图 7 FIR 滤波后信号与 80Hz 信号对比通过图 7 可以发现，滤波后的信号大致反应了 80Hz 信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了 10Hz，200Hz 的信号。为了进一步说明问题，绘制滤波后信号的频谱图，如图 8 所示。从图 8 可以看出，随着 N的增大，80Hz 幅值衰减的程度变小，会趋于至原幅值的一半，10Hz 和 200Hz 信号幅值衰减程度变大，滤波效果更加明显。图 8 FIR 滤波后频谱（N = 8，
6、30， 80， 800） 80Hz(2)用 FIR 滤波（附录 4）选择带通频率为 40120Hz， 即 F1 = 40Hz， F2 = 120Hz， ，根据公式，可得相应的滤波因子，编写相应的程序，可得到滤波后的信号，如图 9 所示。图 9 IIR 滤波后信号与 80Hz 信号对比通过图 9 可以发现，滤波后的信号大致反应了 80Hz 信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了 10Hz，200Hz 的信号。在滤波信号开始阶段，会出现一较大的波动，该波动会随滤波的进行而消失。为了便于说明问题，绘制出滤波后信号的频谱，如图 10 所示。从图 10 可以看出，滤波后的信号幅值
7、基本与原幅值一样，10Hz 信号和 200Hz 信号的幅值衰减较大，滤波效果比 FIR滤波效果好。图 10 IIR 滤波后频谱（4 ） 、原信号波形加 5%的白噪声信号，进行滤波（附录 5）解：利用 matlab 的 awgn 函数，对原信号添加 50%的白噪声，命令如下：y = awgn(x,SNR) 在信号 x 中加入高斯白噪声。信噪比 SNR，本例中，SNR = 2。加入白噪声信号之后的信号波形如图 11 所示。图 11 添加白噪声信号之后的信号波形采用低通 IIR 滤波器，滤去 80Hz,200Hz 信号，保留 10Hz 信号，滤波后信号如图 12 所示。图 12 加白噪声之后滤波信号
8、与 10Hz 信号对比为了便于分析，绘制滤波后的频谱，如图 13 所示。图 13 加入白噪声滤波之后频谱将该频谱与未加白噪声的滤波之后的信号的频谱（图 6）对比可以发现，加入白噪声之后，滤波之后的信号同样被白噪声影响，并未滤去白噪声信号。附录 1%采用 FIR 滤波器 低通滤波器%滤波效果和 N，F 有关clc;clear;Dt = 0.0001;t = 0:Dt:0.5;xt = (t)sin(2*pi*10*t) + sin(2*pi*80*t) + sin(2*pi*200*t);F = 50; %低通滤波的频率；N = 80; %滤波器项数；fi_fir = sin(2*pi*F*(1
9、:N)*Dt)./(pi*(1:N); %滤波因子 f0_fir = 2*F*Dt;f_fir = f0_fir fi_fir; %得到的滤波因子序列for k = 1:length(t)k_t = Dt*(k-N):k);x_k_t = xt(k_t); w = conv(f_fir,x_k_t);y(k) = w(length(f_fir);endfigure;plot(t,y,r);hold on;plot(t,sin(2*pi*10*t);title(滤波后信号与 10Hz 信号对比);xlabel(时间 t);ylabel(xt);legend(滤波后,y = sin(2*pi*10
10、*t);% % 采用汉宁窗口对泄漏进行修正hold on;fi_hanning = 0.5*fi_fir.*(1 + cos(pi*(1:N)/N);f_hanning = f0_fir fi_hanning;for k = 1:length(t)k_t = Dt*(k-N):k);x_k_t = xt(k_t); w = conv(f_hanning,x_k_t);y_hanning(k) = w(length(f_hanning);endfigure;hold onplot(t,y,b-,t,y_hanning,g)title(采用汉宁窗口修正对比);xlabel(时间 t);ylabel
11、(xt);legend(未修正,修正后);% %频谱分析 幅值频谱subplot(4,1,4);N = length(t);Y = fft(y,N)/N*2; ff = 1/Dt/N*(0:1:N-1);plot(ff(1:N/20),abs(Y(1:N/20);title(滤波后频谱 N = 800)xlabel(频率(Hz)ylabel(H(f);附录 2%采用二阶 IIR 滤波器 低通滤波器clc;clear;%绘制信号波形Dt = 1/1000;t = 0:Dt:0.5;xt = sin(2*pi*10*t) + sin(2*pi*80*t) + sin(2*pi*200*t);F =
12、 50; %低通滤波的频率；omega = tan(pi*F*Dt);f0 = omega2/(1+sqrt(2)*omega+omega2);f1 = 2*omega2/(1+sqrt(2)*omega+omega2);f2 = omega2/(1+sqrt(2)*omega+omega2);g1 = -2*(1- omega2)/(1+sqrt(2)*omega+omega2);g2 = (1-sqrt(2)*omega+omega2)/(1+sqrt(2)*omega+omega2);y(1) = 0;y(2) = xt(2);for k = 3:length(t)%x_k = xt(k
13、); x_k_1 = xt(k-1); x_k_2 = xt(k-2); y(k) = f0*xt(k) + f1*xt(k-1) + f2*xt(k-2) - g1*y(k-1) - g2*y(k-2);endplot(t,y)hold on;plot(t,sin(2*pi*10*t);title(滤波后信号与 10Hz 信号对比);xlabel(时间 t);ylabel(xt);legend(滤波后,y = sin(2*pi*10*t);% %频谱分析 幅值频谱N = length(t);Y = fft(y,N)/N*2; ff = 1/Dt/N*(0:1:N-1);plot(ff(1:N
14、/2),abs(Y(1:N/2);title(滤波后频谱)xlabel(频率(Hz)ylabel(H(f);附录 3%fir 滤波器 带通clc;clear;Dt = 0.0001;t = 0:Dt:0.1;xt = (t)sin(2*pi*10*t) + sin(2*pi*80*t) + sin(2*pi*200*t);F1 = 40; F2 = 120;N = 800; %滤波器项数；f0 = 2*Dt*(F2-F1);fi = 2./(pi.*(1:N).*sin(pi*(F2-F1).*(1:N)*Dt).*cos(pi*(F2+F1).*(1:N)*Dt);f = f0 fi;for
15、 k = 1:length(t)k_t = Dt*(k-N):k);x_k_t = xt(k_t); w = conv(f,x_k_t);y(k) = w(length(f);endplot(t,y)hold on;plot(t,sin(2*pi*80*t);title(滤波后信号与 80Hz 信号对比);xlabel(时间 t);ylabel(xt);legend(滤波后,y = sin(2*pi*80*t);subplot(4,1,4);N = length(t);Y = fft(y,N)/N*2; ff = 1/Dt/N*(0:1:N-1);plot(ff(1:N/20),abs(Y(1
16、:N/20);title(滤波后频谱 N = 800)xlabel(频率(Hz)ylabel(H(f);附录 4%iir 滤波器 带通clc;clear;%绘制信号波形Dt = 1/1000;t = 0:Dt:0.2;xt = sin(2*pi*10*t) + sin(2*pi*80*t) + sin(2*pi*200*t);F1 = 40; F2 = 120; omega = tan(pi*(F2-F1)*Dt);beta = cos(pi*(F2+F1)*Dt)/cos(pi*(F2-F1)*Dt);K = 1+sqrt(2)*omega+omega2;f0 = omega2/K; f1
17、= 0; f2 = -2*f0; f3 = 0; f4 = f0;g1 = -2*beta*(2+sqrt(2)*omega)/K;g2 = 2*(1+2*beta2-omega2)/K;g3 = -2*beta*(2-sqrt(2)*omega)/K;g4 = (1-sqrt(2)*omega+omega2)/K;y(1) = xt(1);y(2) = xt(2);y(3) = xt(3);y(4) = xt(4);for k = 5:length(t)y(k) = f0*xt(k) + f1*xt(k-1) + f2*xt(k-2) + f3*xt(k-3) + f4*xt(k-4) -.
18、g1*y(k-1) - g2*y(k-2) - g3*y(k-3) - g4*y(k-4);endplot(t,y)hold on;plot(t,sin(2*pi*80*t);title(滤波后信号与 80Hz 信号对比);xlabel(时间 t);ylabel(xt);legend(滤波后,y = sin(2*pi*80*t);N = length(t);Y = fft(y,N)/N*2; ff = 1/Dt/N*(0:1:N-1);plot(ff(1:N/2),abs(Y(1:N/2);title(滤波后频谱)xlabel(频率(Hz)ylabel(H(f);附录 5%添加高斯白噪声，信号
19、比为 2%采用二阶 IIR 滤波器clc;clear;%绘制信号波形Dt = 1/1000;t = 0:Dt:0.2;xt = sin(2*pi*10*t) + sin(2*pi*80*t) + sin(2*pi*200*t);%添加白噪声yt = awgn(xt,2);F = 30; %低通滤波的频率；omega = tan(pi*F*Dt);f0 = omega2/(1+sqrt(2)*omega+omega2);f1 = 2*omega2/(1+sqrt(2)*omega+omega2);f2 = omega2/(1+sqrt(2)*omega+omega2);g1 = -2*(1- o
20、mega2)/(1+sqrt(2)*omega+omega2);g2 = (1-sqrt(2)*omega+omega2)/(1+sqrt(2)*omega+omega2);%y(1),y(2)如何赋值？y(1) = 0;% y(2) = sin(2*pi*10*Dt);y(2) = yt(2);for k = 3:length(t)%x_k = xt(k); x_k_1 = xt(k-1); x_k_2 = xt(k-2); y(k) = f0*yt(k) + f1*yt(k-1) + f2*yt(k-2) - g1*y(k-1) - g2*y(k-2);endplot(t,y)hold on;plot(t,sin(2*pi*10*t);title(滤波后信号与 10Hz 信号对比);xlabel(时间 t);ylabel(xt);legend(滤波后,y = sin(2*pi*10*t);N = length(t);Y = fft(y,N)/N*2; ff = 1/Dt/N*(0:1:N-1);plot(ff(1:N/2),abs(Y(1:N/2);title(滤波后频谱)xlabel(频率(Hz)ylabel(H(f);
