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 信号的变化，相位
4、一致，幅值衰减 了一部分，说明滤波后，确实去除了 80Hz，200Hz 的信号。在滤波信号开始阶段，会出现 一较大的波动，该波动会随滤波的进行而消失。 为了便于说明问题，绘制出滤波后信号的频谱，如图 6 所示。从图 6 可以看出，滤波后的 信号幅值基本与原幅值一样，且高频信号衰减幅度比较大，滤波效果比 FIR 滤波效果好。 图 6 IIR 滤波后频谱 3、带通滤波，分别用 FIR，IIR 滤波器，保留 80Hz，去除 10Hz 和 200Hz，并画出波形，并 与 10Hz 信号对比。 解：原信号的最大 Fmax = 200Hz， 取： =103 12= 1400=0.0025 满足采样定理。
5、（1 ） 、用 FIR 滤波器（附录 3） 选择带通频率为 40120Hz， 即 F1 = 40Hz， F2 = 120Hz，滤波器项数为 80，根据公式， 可得相应的滤波因子，编写相应的程序，可得到滤波后的信号，如图 7 所示。 图 7 FIR 滤波后信号与 80Hz 信号对比 通过图 7 可以发现，滤波后的信号大致反应了 80Hz 信号的变化，相位一致，幅值衰减 了一部分，说明滤波后，确实去除了 10Hz，200Hz 的信号。 为了进一步说明问题，绘制滤波后信号的频谱图，如图 8 所示。从图 8 可以看出，随着 N 的增大，80Hz 幅值衰减的程度变小，会趋于至原幅值的一半，10Hz 和
6、200Hz 信号幅值衰 减程度变大，滤波效果更加明显。 图 8 FIR 滤波后频谱（N = 8，30， 80， 800） 80Hz (2)用 FIR 滤波（附录 4） 选择带通频率为 40120Hz， 即 F1 = 40Hz， F2 = 120Hz， ，根据公式，可得相应的滤波 因子，编写相应的程序，可得到滤波后的信号，如图 9 所示。 图 9 IIR 滤波后信号与 80Hz 信号对比 通过图 9 可以发现，滤波后的信号大致反应了 80Hz 信号的变化，相位一致，幅值衰减 了一部分，说明滤波后，确实去除了 10Hz，200Hz 的信号。在滤波信号开始阶段，会出现 一较大的波动，该波动会随滤波的
7、进行而消失。 为了便于说明问题，绘制出滤波后信号的频谱，如图 10 所示。从图 10 可以看出，滤波后 的信号幅值基本与原幅值一样，10Hz 信号和 200Hz 信号的幅值衰减较大，滤波效果比 FIR 滤波效果好。 图 10 IIR 滤波后频谱 （4 ） 、原信号波形加 5%的白噪声信号，进行滤波（附录 5） 解：利用 matlab 的 awgn 函数，对原信号添加 50%的白噪声，命令如下： y = awgn(x,SNR) 在信号 x 中加入高斯白噪声。 信噪比 SNR，本例中，SNR = 2。 加入白噪声信号之后的信号波形如图 11 所示。 图 11 添加白噪声信号之后的信号波形 采用低通
8、 IIR 滤波器，滤去 80Hz,200Hz 信号，保留 10Hz 信号，滤波后信号如图 12 所 示。 图 12 加白噪声之后滤波信号与 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
9、*80*t) + sin(2*pi*200*t); F = 50; %低通滤波的频率； N = 80; %滤波器项数； fi_fir = sin(2*pi*F*(1: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); end figure; plot(t,y,r); hold on; plot
10、(t,sin(2*pi*10*t); title(滤波后信号与 10Hz 信号对比); xlabel(时间 t); ylabel(xt); legend(滤波后,y = sin(2*pi*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
11、) = w(length(f_hanning); end figure; hold on plot(t,y,b-,t,y_hanning,g) title(采用汉宁窗口修正对比); xlabel(时间 t); ylabel(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); 附录
12、 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 = 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)
13、/(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); 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); end plot(t,y) hold on; plot(t,sin(2*pi*10*t); title(滤波后信号与
14、 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); 附录 3 %fir 滤波器 带通 clc;clear; Dt = 0.0001; t = 0:Dt:0.1; xt = (t)sin(2*pi*10*t) + sin(2*pi*80*
15、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 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); end plot(t,y) hold on; plot(t,sin(2*pi*80*t
16、); 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:N/20); title(滤波后频谱 N = 800) xlabel(频率(Hz) ylabel(H(f); 附录 4 %iir 滤波器 带通 clc;clear; %绘制信号波形 Dt = 1/1000; t = 0:Dt:0.2; x
17、t = 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 = 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*b
18、eta*(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) -. g1*y(k-1) - g2*y(k-2) - g3*y(k-3) - g4*y(k-4); end plot(t,y) hold on; plot(t,sin(2*pi*80*t); ti
19、tle(滤波后信号与 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 %添加高斯白噪声，信号比为 2 %采用二阶 IIR 滤波器 clc;clear; %绘制信号波形 Dt = 1/1000; t = 0:Dt:0.2; xt = sin(2*pi*
20、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- omega2)/(1+sqrt(2)*omega+omega2); g2 = (1-sqrt(2)*omega+omega2)/(1
21、+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); end plot(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);
