哈尔滨工业大学-试验方法与数字信号处理大作业

Harbin Institute of Technology大作业一课程名称： 试验方法与数字信号处理 院 系： 机械电子 班 级： 15S0825 学 号： 姓 名： 哈尔滨工业大学给出信号xt=sin210t+sin280t+ sin2200t 1. 绘出信号波形。利用matla软件，绘制出的原信号波形如图1所示。图1 原波形信号2. 低通滤波，分别用FIR，IIR滤波器，保留10Hz，去除80Hz和200Hz，并画出波形，并与10Hz信号对比。解：原信号的最大Fmax = 200Hz， 取：t=10-3 12Fmax= 1400=0.0025 此时，满足采样定理。（1）、用FIR滤波器（附录1） 选择低通滤波的截止频率为50Hz，滤波器项数为80，通过FIR滤波器公式，可得到滤波后的信号。编写matlab程序，对比滤波后信号和10Hz信号，如图2所示。图2 FIR滤波后信号与10Hz信号对比 通过图2可以发现，滤波后的信号大致反应了10Hz信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了80Hz，200Hz的信号。为了进一步说明问题，绘制滤波后信号的频谱图，如图3所示。从图3可以看出，随着N的增大，10Hz信号幅值衰减的程度变小，会趋于至原幅值的一半，其余信号幅值衰减的程度变大，滤波效果更加明显。图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，200Hz的信号。在滤波信号开始阶段，会出现一较大的波动，该波动会随滤波的进行而消失。为了便于说明问题，绘制出滤波后信号的频谱，如图6所示。从图6可以看出，滤波后的信号幅值基本与原幅值一样，且高频信号衰减幅度比较大，滤波效果比FIR滤波效果好。图6 IIR滤波后频谱3、带通滤波，分别用FIR，IIR滤波器，保留80Hz，去除10Hz和200Hz，并画出波形，并与10Hz信号对比。解：原信号的最大Fmax = 200Hz， 取：t=10-3 12Fmax= 1400=0.0025 满足采样定理。（1）、用FIR滤波器（附录3） 选择带通频率为40120Hz， 即F1 = 40Hz， F2 = 120Hz，滤波器项数为80，根据公式，可得相应的滤波因子，编写相应的程序，可得到滤波后的信号，如图7所示。图7 FIR 滤波后信号与80Hz信号对比 通过图7可以发现，滤波后的信号大致反应了80Hz信号的变化，相位一致，幅值衰减了一部分，说明滤波后，确实去除了10Hz，200Hz的信号。为了进一步说明问题，绘制滤波后信号的频谱图，如图8所示。从图8可以看出，随着N的增大，80Hz幅值衰减的程度变小，会趋于至原幅值的一半，10Hz和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的信号。在滤波信号开始阶段，会出现一较大的波动，该波动会随滤波的进行而消失。为了便于说明问题，绘制出滤波后信号的频谱，如图10所示。从图10可以看出，滤波后的信号幅值基本与原幅值一样，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 加白噪声之后滤波信号与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: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*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(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 = 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); 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/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 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: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 = 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) -. 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%添加高斯白噪声，信号比为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- omega2)/(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);