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

哈尔滨工业大学-试验方法与数

字信号处理大作业

HarbinInstituteofTechnology

大作业一

解程名椰：武验方法芍数字憎号处理

院系：机楸电子

施级：15so825

学号：...

姓.名：___________

哈尔澳工业大学

给出信号x(t)=sin(2?r-101t)+sin(2?r-801

t)+sin(2?r-200-t)

1.绘出信号波形。

利用matla软件，绘制出的原信号波形如图1

所示。

原信号温影

3000100200300400500600700B00901

t

图1原波形信号

2.低通滤波，分别用FIR,IIR滤波器，保留

10Hz,去除80Hz和200Hz,并画出波形，并与

10Hz信号对比。

解：原信号的最大Fmax=200Hz,取：

11

At=IO_3<——=——=0.0025

2Fmax400

此时，满足采样定理。

(1)、用FIR滤波器(附录1)

选择低通滤波的截止频率为50Hz,滤波器

项数为80,通过FIR滤波器公式，可得到滤波

后的信号。编写matlab程序，对比滤波后信号

和10Hz信号，如图2所示。

08

06

04

02

0

-02

•0.4

-06

­08

图2FIR滤波后信号与10Hz信号对比

通过图2可以发现，滤波后的信号大致反应

了10Hz信号的变化，相位一致，幅值衰减了一

部分，说明滤波后，确实去除了80Hz,200Hz的

信号。

为了进一步说明问题，绘制滤波后信号的频谱

图，如图3所示。从图3可以看出，随着N的增

大,10Hz信号幅值衰减的程度变小，会趋于至

原幅值的一半，其余信号幅值衰减的程度变大，

滤波效果更加明显。

注/盛疆用N«8

图3FIR滤波后频谱(N=8,30,80,800)

10Hz

尝试用汉宁窗口对泄漏进行修正，修正前后

的波形如图4所示。

采用汉宁窗口修正对比

《11Ittllll1

'00050101502026030950404505

付间I

图4采用汉宁窗口修正

(2)、用IIR滤波器(附录2)

选择低通滤波的截止频率为50Hz的二阶HR

滤波器，根据相关公式，可以得到HR滤波器的

滤波因子，进而可得到滤波后的信号。

编写matlab程序，对比滤波后信号和10Hz

信号，如图5所示。

图5HR滤波后信号与10Hz信号对比

通过图5可以发现，滤波后的信号大致反应

了10Hz信号的变化，相位一致，幅值衰减了一

部分，说明滤波后，确实去除了80Hz,200Hz的

信号。在滤波信号开始阶段，会出现一较大的波

动，该波动会随滤波的进行而消失。

为了便于说明问题，绘制出滤波后信号的频谱，

如图6所示。从图6可以看出，滤波后的信号幅

值基本与原幅值一样，且高频信号衰减幅度比较

大,滤波效果比FIR滤波效果好。

算率(Hz)

图6HR滤波后频谱

3、带通滤波，分别用FIR,IIR滤波器，保留

80Hz,去除10Hz和200Hz,并画出波形，并与

10Hz信号对比。

解：原信号的最大取x=200Hz,取:

211

△t=10-3<———=0.0025

^2LFmax400

满足采样定理。

(1)、用FIR滤波器(附录3)

选择带通频率为40~120Hz,即Fl=40Hz,

F2=120Hz,滤波器项数为80,根据公式，可

得相应的滤波因子，编写相应的程序，可得到滤

波后的信号，如图7所示。

图7FIR滤波后信号与80Hz信号对比

通过图7可以发现，滤波后的信号大致反应

了80Hz信号的变化，相位一致，幅值衰减了一

部分，说明滤波后，确实去除了10Hz,200Hz的

信号。

为了进一步说明问题，绘制滤波后信号的频谱

图，如图8所示。从图8可以看出，随着N的增

大,80Hz幅值衰减的程度变小，会趋于至原幅

值的一半，10Hz和200Hz信号幅值衰减程度变

大，滤波效果更加明显。

图8FIR滤波后频谱(N=8,30,80,800)

80Hz

(2)用FIR滤波(附录4)

选择带通频率为40~120Hz,即Fl=40Hz,

F2=120Hz,,根据公式，可得相应的滤波因子,

编写相应的程序，可得到滤波后的信号，如图9

所示。

图9IIR滤波后信号与80Hz信号对比

通过图9可以发现，滤波后的信号大致反应

了80Hz信号的变化，相位一致，幅值衰减了一

部分，说明滤波后，确实去除了10Hz,200Hz的

信号。在滤波信号开始阶段，会出现一较大的波

动，该波动会随滤波的进行而消失。

为了便于说明问题，绘制出滤波后信号的频谱，

如图10所示。从图10可以看出，滤波后的信号

幅值基本与原幅值一样，10Hz信号和200Hz信

号的幅值衰减较大，滤波效果比FIR滤波效果

好。

身率（H：）

图10HR滤波后频谱

(4)、原信号波形加5%的白噪声信号，进行滤

波(附录5)

解：利用matlab的awgn函数，对原信号添

加50%的白噪声，命令如下：

y=awgn(x,SNR)在信号x中加入高斯白噪

声。

信噪比SNR,本例中，SNR=

2o

加入白噪声信号之后的信号波形如图11所

Zj\O

图11添加白噪声信号之后的信号波形

采用低通HR滤波器，滤去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

.%得

到的滤波因子序列

fork=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');

holdon;

plot(t,sin(2*pi*10*t));

title('滤波后信号与10Hz信号对比');

xlabelC时间t');

ylabelCxt');

legend('滤波后'，'y=sin(2*pi*10*t,);

%%采用汉宁窗口对泄漏进行修正

holdon;

fi_hanning=0.5*fi_fir.*(1+

cos(pi*(l:N)/N));

f_hanning=[f0_firfi_hanning];

fork=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));

end

figure;

holdon

plot(t,y,'b,，t,y_hanning,Jg)

titleC采用汉宁窗口修正对比');

xlabel(，时间;

ylabel('xt');

legendC未修正'，'修正后');

%%频谱分析幅值频谱

subplot(4,1,4);

N=length(t);

Y=fft(y,N)/N*2;.

ff=1/Dt/N*(O:1:N-1);

plot(ff(l:N/20),abs(Y(l:N/20)));

title('滤波后频谱N=800')

xlabelC频率(Hz)')

ylabelCH(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);

fO=omega人2/(l+sqrt⑵*omega+omega-2);

fl=

2*omega2/(1+sqrt(2)*omega+omega"2);

f2=omega"2/(1+sqrt(2)*omega+omegaZK2);

gl=-2*(1-

omega"2)/(1+sqrt(2)*omega+omega"2);

g2=

(1-sqrt(2)*omega+omegaz'2)/(1+sqrt(2)*omeg

a+omega-2);

y(l)=0；

y(2)=xt(2);

fork=3:length(t)

%x_k=xt(k);x_k_l=xt(k-1)

x_k_2=xt(k-2);

y(k)=f0*xt(k)+fl*xt(k-1)+

f2*xt(k-2)-gl*y(k-1)-g2*y(k-2);

end

plot(t,y)

holdon;

plot(t,sin(2*pi*10*t));

title('滤波后信号与10Hz信号对比');

xlabelC时间t');

ylabel('xt');

legend('滤波后'，'y=sin(2*pi*10*t，);

%%频谱分析幅值频谱

N=length(t);

Y=fft(y,N)/N*2;

ff=1/Dt/N*(O:1:N-1)

plot(ff(l:N/2),abs(Y(l:N/2)));

titleC滤波后频谱')

xlabelC频率(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);

Fl=40;F2=120;

N=800;%滤波器项数；

fO=2*Dt*(F2-F1);

fi=

2./(pi.*(1:N)).*sin(pi*(F2-Fl).*(1:N)*Dt)

.*cos(pi*(F2+Fl).*(1:N)*Dt);

f=[fOfi];

fork=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)

holdon;

plot(t,sin(2*pi*80*t));

titleC滤波后信号与80Hz信号对比');

xlabelC时间t');

ylabel(，xt')；

legendC滤波后'，'y=sin(2*pi*80*t，);

subplot(4,1,4);

N=length(t);

Y=fft(y,N)/N*2;

ff=1/Dt/N*(O:1:N-1)

plot(ff(l:N/20),abs(Y(l:N/20)));

title('滤波后频谱N=800')

xlabelf频率(Hz)')

ylabelCH(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);

Fl=40;F2=120;

omega=tan(pi*(F2-Fl)*Dt);

beta=

cos(pi*(F2+F1)*Dt)/cos(pi*(F2-Fl)*Dt);

K=1+sqrt(2)*omega+omega"2;

fO=omega-2/K;fl=0;f2=-2*f0;

f3=0;f4=fO;

gl=-2*beta*(2+sqrt(2)*omega)/K;

g2=2*(l+2*beta"2-omega"2)/K;

g3=-2*beta*(2-sqrt(2)*omega)/K;

g4=(1-sqrt(2)*omega+omega"2)/K;

y(l)=xt(l)；

y(2)=xt(2);

y(3)=xt(3);

y(4)=xt(4);

fork=5:length(t)

y(k)=f0*xt(k)+fl*xt(k-1)+

f2*xt(k-2)+f3*xt(k-3)+f4*xt(k-4)

gl*y(k-1)-g2*y(k-2)-g3*y(k-3)

-g4*y(k-4);

end

plot(t,y)

holdon;

plot(t,sin(2*pi*80*t));

title('滤波后信号与80Hz信号对比');

xlabel('时间t');

ylabel('xt');

legendC滤波后‘，'y二sin(2*pi*80*t,);

N=length(t);

Y=fft(y,N)/N*2;

ff=1/Dt/N*(O:1:N-1);

plot(ff(l:N/2),abs(Y(l:N/2)));

titleC滤波后频谱')

xlabelC频率(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);