您的当前位置：首页数字信号处理实验报告

数字信号处理实验报告

来源：爱够旅游网

数字信号处理实验报告

Part1 分别用脉冲响应不变法和双线性变换法将图中的RC 低通滤波器转换为数字滤波器，RC=0.001。

Ha(s)a1,a asRCa1eaT脉冲响应不变法: H(z)双线性变换法：

z1

a1(1z1)H2(z)Ha(s)|21z11s1az12T1za1aTaT2,a2aT2aT2

利用matlab绘图,分别取T=0.0005s、T=0.001s、T=0.002s

1. 模拟滤波器的频响曲线

aHas

as程序及结果如下： RC=0.001; a=1/RC;

f=0:10:1000;

Hajf=a./(a+j*2*pi*f); A=abs(Hajf); plot(f,A)

10.90.80.70.60.50.40.30.20.101002003004005006007008009001000

2. 分别用脉冲响应不变法和双线性变换法得到数字滤波器的系统函数H(z)；

冲激响应不变法得到函数：

H(z)a1eaTz1

aTaTz1双线性变换法得到函数： 1(2aT)(aT2)z3.当采样间隔T=0.0005s、T=0.001s、T=0.002s，分别画出两种数字滤波器的频响曲线，示例程序如下：

%模拟 RC=0.001; a=1/RC; f=0:10:1000; Hajf=a./(a+j*2*pi*f); A=abs(Hajf); f=f/max(f); figure(1) plot(f,A) hold on

%数字 a=1000; T=0.0005; B=[a]; A=[1,-exp(-a*T)]; [H,W]=freqz(B,A,100); n=0:99; n=n/max(n); H=abs(H); H=H/max(H); figure(1) plot(n,H,'r')

4.将模拟滤波器的频响曲线放置在同一坐标中(蓝色曲线)，T不同时需要模拟频率f的取值范围也不相同；曲线如下，其中红色曲线为数字滤波器，蓝色曲线为模拟滤波器。

冲激相应不变法：

10.90.8T=0.0005s

0.70.60.50.40.30.20.100.10.20.30.40.50.60.70.80.9110.90.80.70.60.50.400.10.20.30.40.50.60.70.80.91

T=0.001s

10.950.90.850.80.750.70.650.60.550.500.10.20.30.40.50.60.70.80.91T=0.002s

双线形变换法

10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.60.70.80.91T=0.0005s

10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.60.70.80.91T=0.001s

10.90.80.70.60.50.40.30.20.1000.10.20.30.40.50.60.70.80.91T=0.002s

比较两种方法设计的特点：

在通带范围内，这两种设计方法的响应曲线与模拟滤波器的响应曲线基本吻合。但在阻带范围内，这两种设计方法的响应曲线与模拟滤波器的响应曲线有较大偏差。脉冲响应不变法在高频时性能不如模拟滤波器，即阻带衰减慢；而双线性变换法在高频时性能比模拟滤波器好，阻带衰减快！

脉冲响应不变法: 频率坐标变换是线形的，由于由s到z是多值映射，造成高频混叠。在w = 附近的高频混叠。从图中可以看出，减小采样时间可以减小频率混叠。

双线性变换法: 频率坐标变换是非线性的，由s到z的映射是单值的，采用非线性频率压缩的方法，因此消除了频谱高频混叠。但同时也会造成频响曲线形状的偏离。采样时间越小，频响曲线形状偏离减小。

Part2 用巴特沃斯模拟滤波器设计一个IIR 数字低通滤波器，滤波器的技术要求为：

p0.2rad,ap1dB,s0.3rad,as15dB,T1s1、利用脉冲响应不变法设计：源程序为： clear

T=1;Wp=0.2*pi/T;Rp=1;Ws=0.3*pi/T;Rs=15; [N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');

Wc=Wp/((10^(0.1*abs(Rp))-1)^(1/(2*N))); [z,p,k]=buttap(N); [Bp,Ap]=zp2tf(z,p,k); [Bs,As]=lp2lp(Bp,Ap,Wc); [Bz,Az]=impinvar(Bs,As,T); [H,W]=freqz(Bz,Az,100); w=W/pi;h=abs(H);

subplot(211);plot(w,h); grid on;

hh=20*log10(h);

subplot(212);plot(w,hh); grid on;

所求得的H（z）为：

画出的数字滤波器的幅度响应和增益曲线如下：

10.80.60.40.2000.10.20.30.40.50.60.70.80.910-20-40-60-8000.10.20.30.40.50.60.70.80.91

2、利用双线性变换法设计：程序为： clear;

T=1;Wp=2*tan(0.1*pi)/T;Rp=1;Ws=2*tan(0.15*pi)/T;Rs=15; [N,Wc]=buttord(Wp,Ws,Rp,Rs,'s'); [z,p,k]=buttap(N); [Bp,Ap]=zp2tf(z,p,k); [Bs,As]=lp2lp(Bp,Ap,Wc); [Bz,Az]=bilinear(Bs,As,T); [H,W]=freqz(Bz,Az,100); w=W/pi;h=abs(H); subplot(211); plot(w,h); grid on;

axis([0 1 0 1]) hh=20*log10(h); subplot(212); plot(w,hh); grid on;

axis([0 1 -80 0])所求得的H（z）为：

0.00070.0044z10.0111z20.0148z30.0111z40.0044z50.0007z6H(z)1234561.00003.1836z4.6222z3.7795z1.8136z0.4800z0.04z画出的数字滤波器的幅度响应和增益曲线如下：

10.80.60.40.2000.10.20.30.40.50.60.70.80.910-20-40-60-8000.10.20.30.40.50.60.70.80.91

比较两种设计方法的特点：

在通带范围内，频响特性曲线无太大区别；在过渡带和阻带内，差别很大。脉冲响应不变法: 频率坐标变换是线形的，有高频混叠。

双线性变换法: 由于 Ω与ω之间的非线性关系，数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变。

Part3 用窗函数法设计FIR 数字低通滤波器，截止频率

c0.25源程序为： clear;

N=33; %%% 窗口长度 w=0.25*pi; %% 截止频率wc n=0:N-1;

m=n-0.5*(N-1)+eps;

hd=sin(w*m)./(pi*m); %% 理想低通滤波器单位脉冲响应 figure;

subplot(2,3,1); stem(n,hd);

axis([0 N-1 -0.1 0.3]); hold on;

xlabel('n');

title('理想低通滤波器的hd(n)'); B=boxcar(N); %%矩形窗 %B=hanning(N); %%汉宁窗 %B=hamming(N); %%海明窗 %B=blackman(N); %%布莱克曼窗 subplot(2,3,4); stem(n,B);

axis([0 N-1 0 1]); hold on;

xlabel('n');

title('窗函数w(n)'); string=['N=',num2str(N)];

h=hd.*(B)'; %% 得到FIR数字滤波器(理想滤波器乘以窗函数) subplot(2,3,2); n=0:N-1;

stem(n,h,'.');

axis([0 N-1 -0.1 0.3]); hold on;

xlabel('n');

title('实际低通滤波器的h(n)'); text((0.3*N),0.27,string);

[H,W]=freqz(h,1,1024,'whole'); %% 求其频率响应 mag=abs(H);

db=20*log10(mag/max(mag)); %% 增益(dB) mag=mag/mag(1); %% 得到归一化幅值特性

pha=angle(H); %% 得到相位 subplot(2,3,3); plot(W/pi,db);

axis([0 1 -100 0]); title('幅频特性-增益(dB)'); xlabel('w/pi'); grid;

subplot(2,3,5); plot(W/pi,pha); hold on;

title('相频特性'); xlabel('w/pi'); axis([0 1 -4 4]); subplot(2,3,6); plot(W/pi,mag); axis([0 1 0 1.5]); title('幅频特性');

xlabel('w/pi'); grid;

改变窗函数为矩形窗(Boxcar)、汉宁窗(Hanning)、海明窗(Hamming)和布莱克曼窗(Blackman)，每一种窗的长度分别取N=15,33。实验结果为：

1、矩形窗(Boxcar)，N=15,33

N=15

理想低通滤波器的hd(n)实际低通滤波器的h(n)幅频特性-增益(dB)0.30.30N=150.20.10-0.105100.20.10-0.10510-10000.5w/pi幅频特性1-50n窗函数w(n)1420.50-2005n10-40n相频特性1.510.500.5w/pi100.5w/pi1N=33

理想低通滤波器的hd(n)0.30.20.10-0.10102030实际低通滤波器的h(n)0.30.20.10-0.10102030-10000.5w/pi幅频特性1-50N=330幅频特性-增益(dB)

n窗函数w(n)1420.50-20010n2030-40n相频特性1.510.500.5w/pi100.5w/pi1

2、汉宁窗(Hanning)，N=15,33

N=15

理想低通滤波器的hd(n)实际低通滤波器的h(n)幅频特性-增益(dB)0.30.30N=330.20.10-0.1020n窗函数w(n)1420.50-200n20-400.5w/pi10.500.20.10-0.1020n相频特性1.51-10000.5w/pi幅频特性1-50

00.5w/pi13、海明窗(Hamming) , N=15,33

N=15

00.5w/pi14、布莱克曼窗(Blackman)， N=15,33

N=15

00.5w/pi1

滤波器的性能与N 的取值有关 ,N 越大滤波器的性能越好。不同窗函数滤波器性能也不一样。矩形窗有吉布斯效应，增加 N 或者改善窗函数可以减小吉布

斯效应。

1. 矩形窗：窗形状：矩形；性能最差；

N15,s23.47dB；N33,s22.71dB； 2、汉宁窗：窗函数：升余弦窗；性能稍好；

N15,s44.34dB；N33,s44.01dB； 3、海明窗：窗函数：改进的升余弦窗；性能其次；

N15,s47.06dB；N33,s52.14dB； 4、布莱克曼：窗函数：改进的海明窗；性能最好；

N15,s75.02dB；N33,s75.73dB；

Part4 IIR 和 FIR 滤波效果的比较，分别用上述实验中的IIR 和FIR 低通滤波器对输入信号xn=x1+x2+x3 进行滤波，滤除高频干扰信号x3，输出得到低频信号yn=x1+x2。

1. IIR 低通滤波器

利用上面双线性变换法设计的IIR 低通滤波器对xn进行滤波。 clear; T=1;

Wp=(2/T)*tan(0.5*0.2*pi/T); Rp=1;

Ws=(2/T)*tan(0.5*0.3*pi/T); Rs=15;

[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s'); [Bs,As]=butter(N,Wc,'s'); [Bz,Az]=bilinear(Bs,As,T); N=127; n=0:N-1; fs=;

x1=sin(pi*n*6.4/fs); x2=sin(pi*n*9.6/fs); x3=sin(pi*n*28.8/fs); xn=x1+x2+x3;

yn=filter(Bz,Az,xn); subplot(3,1,1) plot(n,xn);

axis([0 N-1 -3 3]); subplot(3,1,2) plot(n,x1+x2);

axis([0 N-1 -2 2]); subplot(3,1,3) plot(n,yn);

axis([0 N-1 -2 2]); 实验结果：

20-2020-220-220406080100120020406080100120020406080100120

2. FIR 低通滤波器

利用上面N=33 的汉宁Hanning 窗设计的FIR 滤波器对xn 滤波。 clear; N=33;

w=0.25*pi; n=0:1:N-1;

m=n-0.5*(N-1)+eps; hd=sin(w*m)./(pi*m); B=hanning(N); h=hd.*(B)'; N=127;

n=0:1:N-1; fs=;

x1=sin(pi*n*6.4/fs); %% 0.1*pi x2=sin(pi*n*9.6/fs); %% 0.15*pi x3=sin(pi*n*28.8/fs); %% 0.45*pi xn=x1+x2+x3;

yn=filter(h,1,xn); subplot(3,1,1)

plot(n,xn);

axis([0 N-1 -3 3]); subplot(3,1,2) plot(n,x1+x2);

axis([0 N-1 -2 2]); subplot(3,1,3) plot(n,yn);

axis([0 N-1 -2 2]); 实验结果如下：

20-2020-220-220406080100120020406080100120020406080100120

2. 两种滤波器输出的细微差别：

FIR_LPF相对IIR_LPF而言有一定的延时。 3. FIR 滤波器延迟的大小：

16 * 1/fs = 0.25 s

Part5 滤波器的实现 1. 利用差分方程

y(n)bix(ni)aiy(ni)

i0i1MN对上面双线性变换法设计的IIR滤波器，自己编写滤波程序（只允许包含加减乘除运算），对Part4 中的信号xn 滤波，结果应与filter(Bz, Az, xn)函数结果基本一致，要注意滤波器系数ai, bi 的下标；

解：

利用双线性变换法设计的IIR滤波器：

0.00070.0044z10.0111z20.0148z30.0111z40.0044z50.0007z6H(z)1.00003.1836z14.6222z23.7795z31.8136z40.4800z50.04z6差分方程如下：

y(n)3.1816y(n1)4.6222y(n2)3.7795y(n3)1.8136y(n4)0.4800y(n5)0.04y(n6)0.0007x(n)0.0044x(n1)0.0111x(n2)0.0148x(n3)0.0111x(n4)0.0044x(n5)0.0007x(n6)

源程序如下： clear N=127; n=0:N-1; fs=;

x1=sin(pi*n*6.4/fs); x2=sin(pi*n*9.6/fs); x3=sin(pi*n*28.8/fs); xn=x1+x2+x3; y(1:6)=0; for i=7:127

y(i)=0.0007*xn(i)+0.0044*xn(i-1)+0.0111*xn(i-2)+0.0148*xn(i-3)+0.0111*xn(i-4)+0.0044*xn(i-5)+0.0007*xn(i-6)...

+3.1836*y(i-1)-4.6222*y(i-2)+3.7795*y(i-3)-1.8136*y(i-4)+0.48*y(i-5)-0.04*y(i-6); end

N=length(y); n=1:1:N; plot(n,y);

axis([0 N-1 -2 2]);

21.510.50-0.5-1-1.5-2020406080100120

2. 在Matlab 中调用fft 函数，用FFT 实现FIR 滤波，FIR 低通滤波器为Part3 中的汉宁窗N=33，对Part4 中的信号xn 滤波；源程序如下： clear; N=33;

w=0.25*pi; n=0:N-1;

m=n-0.5*(N-1)+eps; hd=sin(w*m)./(pi*m); B=hanning(N); h=hd.*(B)'; N=127; n=0:N-1; fs=;

x1=sin(pi*n*6.4/fs); %% 0.1*pi x2=sin(pi*n*9.6/fs); %% 0.15*pi x3=sin(pi*n*28.8/fs); %% 0.45*pi xn=x1+x2+x3;

L=length(xn)+length(h);

hn=[h,zeros(1,L-length(xn))]; xn=[xn,zeros(1,L-length(h))]; Hk=fft(hn,L); Xk=fft(xn,L); Yk=Hk.*Xk;

yn=ifft(Yk,L); ynws=yn(1:127); N=length(ynws);

n=1:1:N;

plot(n,ynws);

axis([0 N-1 -2 2]);

21.510.50-0.5-1-1.5-2020406080100120

对比实验结果可以看出,实验结果与 part4 完全相同.

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文