《数字信号处理》课程研究性学习报告
IIR和FIR滤波器设计专题研讨
【目的】
(1)掌握根据滤波器指标设计IIR和FIR数字滤波器的原理和方法。
(2)熟悉通过IIR和FIR数字滤波器进行实际系统设计的方法。
(3)培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。
【研讨题目】
1．设计一个数字滤波器，在频率低于的范围内，通带衰减不低于0.75dB。在频率和之间，阻带衰减至少为40dB。
(1)试求满足这些条件的最低阶Butterworth滤波器。
%butterworth单独计算wc双线性
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
wp=0.325;ws=0.726;
N=buttord(wp,ws,Ap,As,''s'');
wc=wp/(10^(0.1Ap)-1)^(1/2/N);
[numa,dena]=butter(N,wc,''s'');
[numd,dend]=bilinear(numa,dena,0.5);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.7476为什么不是0.75wc=wp/(10^(0.1Ap)-1)^(1/2/N);
As=41.6663
N=7.0000
%butterworth单独计算wc脉冲响应不变
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
Fs=1;
wp=WpFs;ws=WsFs;
N=buttord(wp,ws,Ap,As,''s'');
wc=wp/(10^(0.1Ap)-1)^(1/2/N);
[numa,dena]=butter(N,wc,''s'');
[numd,dend]=impinvar(numa,dena,Fs);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.7500
As=40.9184
N=8.0000
%butterworth公式计算wc双线性
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
wp=0.325;ws=0.726;
[N,wc]=buttord(wp,ws,Ap,As,''s'');
[num,den]=butter(N,wc,''s'');
[numd,dend]=bilinear(num,den,0.5);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.5282
As=40.0454
N=7.0000
%butterworth公式计算wc脉冲响应不变
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
Fs=1;
wp=WpFs;ws=WsFs;
[N,wc]=buttord(wp,ws,Ap,As,''s'');
[num,den]=butter(N,wc,''s'');
[numd,dend]=impinvar(num,den,Fs);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.6167
As=40.0001
N=8.0000
通过比较选用第一个，butterworth单独计算wc双线性
(2)试求满足这些条件的最低阶ChebyshevI滤波器。
%切比雪夫脉冲响应
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
Fs=1;
wp=WpFs;ws=WsFs;
[N,wc]=cheb1ord(wp,ws,Ap,As,''s'');
[num,den]=cheby1(N,Ap,wc,''s'');
[numd,dend]=impinvar(num,den,Fs);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.7500
As=43.9315
N=5.0000
%切比雪夫双线性
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
wp=0.325;ws=0.726;
[N,wc]=cheb1ord(wp,ws,Ap,As,''s'');
[num,den]=cheby1(N,Ap,wc,''s'');
[numd,dend]=bilinear(num,den,0.5);
w=linspace(0,pi,512);
h=freqz(numd,dend,w);
norm=max(abs(h));
numd=numd/norm;
plot(w/pi,20log10(abs(h/norm)));
xlabel(''Normalizedfrequency'');
ylabel(''Gain,dB'');
w=[WpWs];
h=freqz(numd,dend,w);
fprintf(''Ap=%.4f\n'',-20log10(abs(h(1))));
fprintf(''As=%.4f\n'',-20log10(abs(h(2))));
fprintf(''N=%.4f\n'',N);
Ap=0.7415
As=49.4169
N=5.0000
通过比较，选用第二个切比雪夫双线性
(3)自主选择一段带限信号，通过所设计的（1）、（2）两种滤波器，比较各自的输入和输出信号。讨论两种滤波器在结构和性能上的差异。
%3-1
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
wp=0.325;ws=0.726;
N=buttord(wp,ws,Ap,As,''s'');
wc=wp/(10^(0.1Ap)-1)^(1/2/N);
[numa,dena]=butter(N,wc,''s'');
[numd,dend]=bilinear(numa,dena,0.5);
k=0:99;
x=cos(0.1pik)+cos(0.9pik);
a=numd;b=dend;
y=filter(b,a,x);
figure(1)
subplot(2,1,1);
stem(k,x);
subplot(2,1,2);
stem(k,y);
L=512;
h1=fftshift(fft(x,L));
h2=fftshift(fft(y,L));
w=linspace(-10pi,10pi,512);
figure(2)
subplot(2,1,1);
plot(w/pi,abs(h1));
subplot(2,1,2);
plot(w/pi,abs(h2));
%3-2
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
wp=0.325;ws=0.726;
[N,wc]=cheb1ord(wp,ws,Ap,As,''s'');
[num,den]=cheby1(N,Ap,wc,''s'');
[numd,dend]=bilinear(num,den,0.5);
k=0:99;
x=cos(0.1pik)+cos(0.9pik);
a=numd;b=dend;
y=filter(b,a,x);
figure(1)
subplot(2,1,1);
stem(k,x);
subplot(2,1,2);
stem(k,y);
L=512;
h1=fftshift(fft(x,L));
h2=fftshift(fft(y,L));
w=linspace(-10pi,10pi,512);
figure(2)
subplot(2,1,1);
plot(w/pi,abs(h1));
subplot(2,1,2);
plot(w/pi,abs(h2));
2．分别用Hamming窗，Blackman窗和Kaiser窗设计，满足下列指标的FIR低通滤波器：
(1)是否可以用这三种窗函数逼近该指标？
(2)如果可以，试设计滤波器，并对这三种滤波器进行比较。
【题目分析】
本题讨论窗函数法设计数字FIR滤波器。结合课本分析不同窗函数法的设计结果。
窗函数法设计FIR滤波器的基本思想是在时域逼近理想滤波器的单位脉冲响应，首先根据待逼近的理想滤波器的频率响应求出滤波器的单位脉冲响应，再将无限长的hd[k]加窗截短得有限长的序列h[k]。通过分析题目中的滤波器指标并对比Hamming窗，Blackman窗和Kaiser窗这三种窗的近似过渡带宽度、通带阻带衰减等参数可知，三种窗均可设计满足题目要求的低通滤波器。
【FIR模拟滤波器设计的基本方法】
根据所需设计的滤波器，确定线形相位滤波器的类型：可选Ⅰ型和Ⅱ型，be本题中我们选用Ⅰ型；
计算理想低通滤波器的单位脉冲响应hd[k];
用不同窗函数截短以上无限长序列为可实现的有限长序列：h[k]=hd[k]w[k]。
【仿真结果】
【结果分析】
通过分析图形可知，利用三种窗函数截短都达到了题中的滤波器指标，其中blackman窗的阻带衰减最大。
【阅读文献】
数字信号处理
【仿真程序】
%利用Hamming窗设计FIR低通滤波器
Wp=0.4pi;Ws=0.6pi;Ap=0.5;As=45;
N=ceil(7pi/(Ws-Wp));
N=mod(N+1,2)+N;
M=N-1;
w=hamming(N)'';
Wc=(Wp+Ws)/2;
k=0:M;
hd=(Wc/pi)sinc(Wc(k-0.5M)/pi);
h=hd.w;
omega=linspace(0,pi,512);
mag1=freqz(h,[1],omega);
magdb1=20log10(abs(mag1));
%利用Blackman窗设计FIR低通滤波器
Wp=0.4pi;Ws=0.6pi;Ap=0.5;As=45;
N=ceil(11.4pi/(Ws-Wp));
N=mod(N+1,2)+N;
M=N-1;
w=blackman(N)'';
Wc=(Wp+Ws)/2;
k=0:M;
hd=(Wc/pi)sinc(Wc(k-0.5M)/pi);
h=hd.w;
omega=linspace(0,pi,512);
mag2=freqz(h,[1],omega);
magdb2=20log10(abs(mag2));
%利用Kaiser窗设计FIR低通滤波器
Ap=0.5;As=45;
Rp=1-10.^(-0.05Ap);Rs=10.^(-0.05As);
f=[0.4,0.6];a=[1,0];dev=[Rp,Rs];
[M,Wc,beta,ftype]=kaiserord(f,a,dev);
M=mod(M,2)+M;
h=fir1(M,Wc,ftype,kaiser(M+1,beta))
omega=linspace(0,pi,512);
mag3=freqz(h,[1],omega);
magdb3=20log10(abs(mag3));
plot(omega/pi,magdb1,''r'',omega/pi,magdb2,''gr'',omega/pi,magdb3,''b'');grid;
legend(''Hamming'',''Blakman'',''Kaiser'');
3．附件给出了一段含有噪声的音频信号。
（1）分析该信号的频谱特点；
（2）利用所学的IIR设计方法设计一个数字低通滤波器对其进行处理，得到有用信息，自主确定各项指标。
（3）利用所学的FIR设计方法设计一个数字低通滤波器对其进行处理，得到有用信息，自主确定各项指标。
（4）试定量比较上述两种滤波器的各项性能，画出能说明性能差异的相关图形，对比并解释。
（5）请尝试采用其它的音频信号，混入不同的噪声，利用所学的滤波方法进行分析，会得到什么样的效果？
【题目分析】
本题讨论用IIR和FIR数字滤波器进行实际系统设计的方法。
【仿真结果】
（1）
（2）
（3）
（4）
数字滤波在数字信号处理中，占有重要的地位。数字滤波包括FIR和IIR两种滤波方式，其中FIR滤波具有很多优点，可以在幅度特性随意设计的同时，保证精确、严格的线性相位，滤波稳定，不会出现递归型结构中的极限振荡等不稳定现象，且误差较小，可采用FFT算法实现，因此运算效率高。设计FIR滤波器常用的方法有窗函数法与频率抽样法，但是这两种方法均不易精确控制通带与阻带的边界频率，所以在实际应用中有一定的局限性。文中用Matlab语言实现了最佳等波纹FIR滤波器的设计，通过比较显示了它在等波纹方脉冲响应方面的优化特性。
(5)
【结果分析】
对几种滤波器得到的音频信号进行对比。
【问题探究】
在IIR滤波器设计过程中，由于利用指标数据的不同，造成裕量的出现。讨论利用不同指标出现裕量对滤波器性能的影响，以及如何有效地利用它？
在FIR滤波器设计过程中，由于所选用窗函数的不同，导致对信号滤波的效果不同。本题的语音信号加入的是双频噪声，讨论如果对信号加入其它形式的噪声，采用何种滤波器更合适。
【仿真程序】
（1）
fs=40000;
x1=wavread(''nihao_final'');
sound(x1,40000);
y1=fft(x1);
subplot(2,1,1);
plot(x1)
title(''滤波前的时域波形'');
subplot(2,1,2);
plot(abs(y1));
title(''滤波前的频谱'');
（2）
y1=wavread(''nihao_final'');
sound(y1,40000);
Wp=0.2pi;Ws=0.4pi;Ap=0.75;As=40;
Fs=1;
wp=WpFs;ws=WsFs;
[N,wn]=buttord(wp,ws,Ap,As,''s'');
[Z,P,K]=buttap(N);
[Bap,Aap]=zp2tf(Z,P,K);
[b,a]=lp2lp(Bap,Aap,wn);
[bz,az]=bilinear(b,a,Fs);
[H,W]=freqz(bz,az);
y2=filter(bz,az,x1);
sound(y2,40000);
z1=fft(y1);
z2=fft(y2);
subplot(2,2,1);plot(y1);
title(''滤波前的时域波形'');
subplot(2,2,2);plot(y2);
title(''滤波后的时域波形'');
subplot(2,2,3);plot(abs(z1));
title(''滤波前的频谱'')
subplot(2,2,4);plot(abs(z2));
title(''滤波后的频谱'')
（3）
fs=40000;
x1=wavread(''nihao_final'');
sound(x1,40000);
wp=0.25pi;
ws=0.3pi;
wdelta=ws-wp;
N=ceil(6.6pi/wdelta);
wc=(0.25+0.3)pi/2;
b=fir1(N,wc/pi,hamming(N+1));
f2=filter(bz,az,x1);
sound(f2,40000);
m1=fft(x1);
m2=fft(f2)
subplot(2,2,1);plot(x1);
title(''滤波前的时域波形'');
subplot(2,2,2);plot(f2);
title(''滤波后的时域波形'');
subplot(2,2,3);plot(abs(m1));
title(''滤波前的频谱'')
subplot(2,2,4);plot(abs(m2));
title(''滤波后的频谱'')
(5)
z1=wavread(''【司马懿】难道真是天命难违'');
sound(z1,30000);
t=0:1/30000;
d=[1cos(2pi5000t)];
z2=z1+d;
sound(z2,30000);
fs=30000;
wp=0.4pi;
ws=0.5pi;
wdelta=ws-wp;
N=ceil(6.6pi/wdelta);
wc=(0.25+0.3)pi/2;
b=fir1(N,wc/pi,hamming(N+1));
freqz(b,1,512)
z3=filter(bz,az,z2);
sound(z3,30000);
subplot(2,2,1)
plot(z1)
title(''原始信号波形'');
n1=fft(z1)
subplot(2,2,2)
plot(abs(n1))
title(''原始信号频谱'');
subplot(2,2,3)
n2=fft(z2);
plot(abs(n2))
title(''加噪后信号频谱'');
subplot(2,2,4)
n3=fft(z3);
plot(abs(n3))
title(''加噪又去噪后的信号频谱'');
|
|
