数字滤波器的设计实验.docx
《数字滤波器的设计实验.docx》由会员分享，可在线阅读，更多相关《数字滤波器的设计实验.docx（25页珍藏版）》请在冰豆网上搜索。
数字滤波器的设计实验
实验二IIR数字滤波器的设计
实验内容及步骤:
数字滤波器的性能指标：
通带临界频率fp、阻带临界频率fr；通带内的最大衰减Ap；阻带内的最小衰减Ar；采样周期T；
（1）、fp=0.3KHz,Ap=0.8dB,fr=0.2KHz,Ar=20dB,T=1ms；设计一Chebyshev高通滤波器；观察其通带损耗和阻带衰减是否满足要求。
程序如下：
fp=300;fr=200;
Ap=0.8;Ar=20;
T=0.001;fs=1/T;
wp=2*pi*fp*T;
wr=2*pi*fr*T;
Wp=2/T*tan（wp/2）;
Wr=2/T*tan（wr/2）;
[N,Wn]=cheb1ord（Wp,Wr,Ap,Ar,'s'）;
[B,A]=cheby1（N,Ap,Wn,'high','s'）;
[num,den]=bilinear（B,A,1/T）;
[h,w]=freqz（num,den）;
plot（w*fs/（2*pi）,20*log10（abs（h）））;%衰减及频率都用归一化的1为单位显示
axis（[0,500,-30,0]）;
title（'Chebyshev高通滤波器'）;
xlabel（'频率'）;
ylabel（'衰减'）;
gridon;
根据下图知道通带损耗与阻带衰减满足要求
（2）、fp=0.2KHz,Ap=1dB,fr=0.3KHz,Ar=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。
比较这两种方法的优缺点。
程序如下：
fp=200;fr=300;
Ap=1;Ar=25;
T=0.001;fs=1/T;
wp=2*pi*fp*T;
wr=2*pi*fr*T;
Wp=2/T*tan（wp/2）;
Wr=2/T*tan（wr/2）;
[N,Wn]=buttord（Wp,Wr,Ap,Ar,'s'）;
[B,A]=butter（N,Wn,'s'）;
[num1,den1]=impinvar（B,A,1/T）;%脉冲响应不变法得出设计的传递函数
[num2,den2]=bilinear（B,A,1/T）;%双线性变换法得出设计的传递函数
[h1,w]=freqz（num1,den1）;
plot（w*fs/（2*pi）,20*log10（abs（h2））,w*fs/（2*pi）,20*log10（abs（h1））,'r.'）;gridon;%衰减及频率都用归一化的1为单位显示
axis（[0,500,-30,0]）;
title（'Butterworth低通滤波器（红线—脉冲响应不变法蓝线—双线性变换法）'）;
xlabel（'ƵÂÊ'）;
ylabel（'Ë¥¼õ'）;
gridon;
优缺点：
采用脉冲响应不变法
优点：
1.h（n）完全模仿模拟滤波器的单位抽样响应时域逼近良好
2线性相位模拟滤波器转变为线性相位数字滤波器
缺点：
1.对时域的采样会造成频域的“混叠效应”，故有可能使所设计数字滤波器的频率响应与原来模拟滤波器的频率响应相差很大
2不能用来设计高通和带阻滤波器。
只适用于限带的低通、带通滤波器
采用双线性变换法
优点：
1避免了频率响应的混迭失真现象
2在特定数字滤波器和特定模拟滤波器处，频率响应是严格相等的，它可以较准确地控制截止频率的位置。
3它是一种简单的代数关系，设计十分方便。
缺点：
1除了零频率附近，w与W之间严重非线性，即线性相位模拟滤波器变为非线性相位数字滤波器
2要求模拟滤波器的幅频响应为分段常数型，不然会产生畸变
3对于分段常数型AF滤波器，经双线性变换后，仍得到幅频特性为分段常数的DF.但在各个分段边缘的临界频率点产生畸变，这种频率的畸变，可通过频率预畸变加以校正
（3）、利用双线性变换法分别设计满足下列指标的Butterworth型和Chebyshev型数字低通滤波器，并作图验证设计结果。
fp=1.2kHz,Ap≤0.5dB,fr=2KHz,Ar≥40dB,fs=8KHz
程序如下：
fp=1200;fr=2000;Ap=0.5;Ar=40;
fs=8000;T=1/fs;wp=2*pi*fp*T;wr=2*pi*fr*T;
Wp=2/T*tan（abs（wp/2））;Wr=2/T*tan（abs（wr/2））;
[N,Wn]=buttord（Wp,Wr,Ap,Ar,'s'）;
[B,A]=butter（N,Wn,'s'）;
[num2,den2]=bilinear（B,A,1/T）;%双线性变换法得到设计的传递函数
[h2,w]=freqz（num2,den2）;
figure
（1）;
subplot（2,1,1）
plot（w/pi*fs/2,10*log10（abs（h2）.^2））;
axis（[0fs/2-500]）;
title（'Butterworth低通滤波器'）;
xlabel（'频率Hz'）;ylabel（'衰减dB'）;gridon;
subplot（2,1,2）
[N,Wn]=cheb1ord（Wp,Wr,Ap,Ar,'s'）;
[B,A]=cheby1（N,Ap,Wn,'s'）;
[num,den]=bilinear（B,A,1/T）;
[h1,w]=freqz（num,den）;
plot（w/pi*fs/2,10*log10（abs（h1）.^2））;
axis（[0fs/2-500]）;
title（'chebyshev低通滤波器'）;
xlabel（'频率Hz'）;ylabel（'衰减dB'）;gridon;
（4）、利用双线性变换法设计一Butterworth型数字带通滤波器，已知fs=30KHz，其等效的模拟滤波器指标为
Ap＜3dB,2KHz＜f≤3KHz,Ar≥5dB,f≥6KHz,Ar≥20dB,f≤1.5KHz
程序如下：
f1=2000;f3=3000;
fsl=1500;fsh=6000;
rp=3;rs=20;
Fs=30000;
wp1=2*pi*f1/Fs;
wp3=2*pi*f3/Fs;
wsl=2*pi*fsl/Fs;
wsh=2*pi*fsh/Fs;
wp=[wp1wp3];ws=[wslwsh];
wap=2*Fs*tan（wp./2）;
was=2*Fs*tan（ws./2）;
[n,wn]=buttord（,'s'）;
[z,p,k]=buttap（n）;
[bp,ap]=zp2tf（z,p,k）;
bw=wap
（2）-wap
（1）;
w0=sqrt（wap
（2）*was
（1））;
[bs,as]=lp2bp（bp,ap,w0,bw）;
[bz1,az1]=bilinear（bs,as,Fs）;
[h,w]=freqz（bz1,az1,256,Fs）;
plot（w,20*log10（abs（h）））;
axis（[08000-300]）;
title（'Butterworth带通滤波器'）;
xlabel（'频率Hz'）;
ylabel（'衰减dB'）;
gridon;
1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？
从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？
答：
在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相
位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。
如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。
2．能否利用公式完成脉冲响应不变法的数字滤波器设计？
为什么？
答：
不能，这样会使得H（z）中的z以对数形式出现，使H（z）的分子分母不再是z的有理多项式，这样会给系统的分析和实现带来很大的困难。
实验三、FIR数字滤波器的设计
实验内容及步骤
（1）N=15,
。
用Hanning窗设计一线性相位带通滤波器，观察它的实际3dB和20dB带宽。
N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N变化的影响；
程序如下：
N1=15;
w1=0.3;
w2=0.5;
wn=[w1,w2];
b1=fir1（N1-1,wn,hanning（N1））;%用Hanning窗作为冲击响应的窗函数
figure
（1）;
freqz（b1,1）;
title（'N=15hanning窗幅频相频特性'）;
N2=45;
b2=fir1（N2-1,wn,hanning（N2））;
figure
（2）;
freqz（b2,1）;
title（'N=45hanning窗幅频相频特性'）;
由图可知：
随着N值的增大，取样值增大，主瓣的宽度减小，但是幅频特性与相频特性曲线的波动频率会增加。
（2）分别改用矩形窗和Blackman窗，设计
（1）中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点；
程序如下：
N1=15;
w1=0.3;w2=0.5;
wn=[w1,w2];
b1=fir1（N1-1,wn,boxcar（N1））;%用boxcar作为冲击响应的窗函数
figure
（1）;
freqz（b1,1）;
title（'N=15boxcar窗幅频相频特性'）;
N2=45;
b2=fir1（N2-1,wn,boxcar（N2））;
figure
（2）;
freqz（b2,1）;
title（'N=45boxca窗幅频相频特性'）;
b3=fir1（N1-1,wn,blackman（N1））;%用blackman作为冲击响应的窗函数
figure（3）;
freqz（b3,1）;
title（'N=15blackman窗幅频相频特性'）;
N2=45;
b4=fir1（N2-1,wn,blackman（N2））;
figure（4）;
freqz（b4,1）;
title（'N=45blackman窗幅频相频特性'）;
分析：
由图可知，在三个窗函数中，Blackman窗的衰减性最好，矩形窗的衰减最差。
，而矩形窗的过渡带带宽最小，Blackman窗的过渡带带宽最大。
Hanning窗的性能比较平均。
其实，改善阻带衰减的一种方法是加宽过渡带宽，以牺牲过渡带宽换取阻带衰减的增大。
也就是以增加主瓣宽度为代价来降低旁瓣。
（3）用Kaiser窗设计一专用线性相位滤波器，N=40，
如图，当β0=4，6，10时，分别设计，比较它们的幅频和相频特性，注意β0取不同值时的影响；
程序如下：
N=40;
beta=[4,6,10];
wn=[0.2,0.4,0.6,0.8];
s=['rgb'];
fori=1:
3
window=kaiser（N,beta（i））;
b1=fir1（N-1,wn,window）;
[H,w]=freqz（b1,1）;
magH=abs（H）;angH=angle（H）;
subplot（2,1,1）
plot（w/pi,magH,[s（i）]）;holdon;gridon
xlim（[0,1]）;
ylim（[0,1]）;
xlabel（'w/（2*pi）'）;ylabel（'|H（e^jw）|'）
title（'N=40幅度部分'）;
legend（'β=4','β=6','β=10'）;
subplot（2,1,2）;
plot（w/pi,angH,[s（i）]）;holdon;gridon;
xlim（[0,1]）;
xlabel（'w/（2*pi）'）;ylabel（'Φ（w）'）;
title（'N=40相频部分'）;
end
由下图可知，随着β的增大，幅频曲线变得越平缓，但其边峰瓣值的衰减越大，相频特性变得越好。
（4）、用频率采样法设计（3）中的滤波器，过渡带分别设一个过渡点，令H（k）=0.5。
比较两种不同方法的结果；
分析：
，
，
，
。
而在通带两侧的过渡带需要插入0.5。
采用第一种相位特性，因此，可以得到以下式子
程序如下：
Hr=[zeros（1,3）,0.5,ones（1,5）,0.5,0,0.5,ones（1,5）,0.5,zeros（1,5）,-0.5,-ones（1,5）,-0.5,0,-0.5,-ones（1,5）,-0.5,0,0];
N=40;
k=0:
N-1;
subplot（2,2,1）;
stem（k,Hr,'.'）;gridon;
xlabel（'n'）,ylabel（'h（n）'）,title（'插值图'）;
wm=2*pi*k./N;
Hd=Hr.*exp（-j*（N-1）*wm*0.5）;
hh=real（ifft（Hd））;
w1=linspace（0,pi-0.1,1000）;
H1=freqz（hh,1,w1）;
amgH=abs（H1）;
dbH=20*log10（amgH）;
angH=angle（H1）;
subplot（2,2,2）;
plot（w1/pi,dbH）,gridon;
title（'幅度响应'）;
xlabel（'w/pi'）,ylabel（'|H（e^jw）|/dB'）
axis（[01-10020]）;
subplot（2,2,3）;
plot（w1/pi,angH）,gridon;
title（'相频特性'）;
xlabel（'w/pi'）,ylabel（'Φ'）
subplot（2,2,4）;
plot（w1/pi,amgH）,gridon;
title（'幅频特性'）;
xlabel（'w/pi'）,ylabel（'|H（e^jw）|'）
（5）、用雷米兹（Remez）交替算法设计（3）中的滤波器，并比较（3）、（4）、（5）三种方法的结果。
程序如下：
f=[00.380.420.780.821];
A=[001100];
weigh=[181];
b=remez（40,f,A,weigh）;
%
[h,w]=freqz（b,1,512,'whole'）;
magh=abs（h）;
dbh=20*log10（magh）;
phah=angle（h）;
subplot（2,1,1）
plot（w/（2*pi）,dbh）;gridon;
xlabel（'w/（2*pi）'）,ylabel（'|H（e^jw）|/dB'）,title（'幅频特性'）;
axis（[01-705]）
subplot（2,1,2）
plot（w/（2*pi）,phah）;gridon;
xlabel（'w/（2*pi）'）,ylabel（'Φ（w）'）,title（'相频特性'）;
四、实验思考
1.定性地说明用本实验程序设计的FIR滤波器的3dB截止频率在什么位置？
它等于理想频率响应Hd（ejω）的截止频率吗？
答：
在0.3-0.5之间，大致等于理想频率响应Hd（ejω）的截止频率。
2.如果没有给定h（n）的长度N，而是给定了通带边缘截止频率ωc和阻带临界频率ωp，以及相应的衰减，你能根据这些条件用窗函数法设计线性相位FIR低通滤波器吗？
答：
能，可以根据通带边缘截止频率ωc和阻带临界频率ωp，以及相应的衰减与所选用的窗函数来大致判断其长度N，从而设计线性相位FIR低通滤波器。
实验四、用MATLAB信号处理工具箱进行滤波器设计训练
针对一个含有5Hz、50Hz和100Hz的混和正弦波信号，设计一个FIR带通滤波器并对其进行滤波。
利用MATLAB实现的三种方法：
程序设计法、FDATool设计法和SPTool设计法。
参数要求：
采样频率fs=100Hz，通带下限截止频率fc1=10Hz，通带上限截止频率fc2=20Hz，过渡带宽6Hz，通阻带波动0.01，采用凯塞窗设计。
方法1：
程序设计法
fsamp=100;
fcuts=[4102026];
mags=[010];
devs=[0.010.010.01];
[n,wn,beta,ftype]=kaiserord（fcuts,mags,devs,fsamp）;
n=n+rem（n,2）;
hh=fir1（n,wn,ftype,kaiser（n+1,beta）,'noscale'）;
[H,f]=freqz（hh,1,1024,fsamp）;
figure
（1）
plot（f,abs（H））,gridon;
xlabel（'f/Hz'）,ylabel（'|H（e^jw）|'）,title（'kaiser数字带通滤波器'）;
用一个序列测试滤波器性能
figure
（2）
N=256;
n=0:
1:
N-1;
t=n/fsamp;
xn=sin（2*pi*5*t）+sin（2*pi*15*t）+sin（2*pi*35*t）;
subplot（2,2,1）
stem（n,xn,'.'）;gridon;
xlabel（'n'）;ylabel（'x（n）'）;title（'序列x（n）'）;
axis（[0256-22]）;
X=fft（xn,N）;
magX=abs（X）;
f=n*fsamp/N;
subplot（2,2,2）;
plot（f（1:
N/2）,magX（1:
N/2））;gridon;
%axis（[0fs01]）;
xlabel（'f/Hz'）;ylabel（'|X（e^jw）|'）;title（'序列x（n）幅频特性'）;
x=xn;
y=filter（hh,1,x）;
Y=fft（y,N）;
magY=abs（Y）;
subplot（2,2,3）;
stem（n,y,'.'）;gridon;
xlabel（'n'）;ylabel（'y（n）'）;title（'滤波后的序列y（n）'）;
axis（[0256-22]）;
subplot（2,2,4）;
plot（f（1:
N/2）,magY（1:
N/2））;gridon;
xlabel（'f/Hz'）;ylabel（'|Y（e^jw）|'）;title（'序列y（n）幅频特性'）;
方法2：
利用FDATool设计法
在MATLAB命令窗口输入FDATool后回车就会弹出FDATool界面。
带通滤波器设计:
已知滤波器的阶数n=38，beta=3.4。
（1）首先在Filter Type中选择Bandpass；
（2）在Design Method选项中选择FIR Window，
（3）接着在Window选项中选取Kaiser，Beta值为3.4；
（4）指定Filter Order项中的Specify order为38；采样频率Fs=100Hz,截止频率Fc1=10Hz,Fc2=20Hz。
（5）设置完以后点击窗口下方的Design Filter，在窗口上方就会看到所设计滤波器的幅频响应，通过菜单选项Analysis还可以看到滤波器的相频响应、组延迟、脉冲响应、阶跃响应、零极点配置等。
（6）设计完成后将结果保存为附件test0802.fda。
在simulink中进行验证，保存为附件test0802.mdl其原理图如下：
滤波前后的图如下：
方法3：
利用SPTool设计法
带通滤波器设计:
已知滤波器的阶数n=38，beta=3.4。
（1） 创建并导入信号源。
在MATLAB命令窗口输入命令：
Fs=100;n = 0:
255;t=n/Fs;
xn = sin（2*pi*t*5）+sin（2*pi*t*15）+sin（2*pi*t*30）; 此时，变量Fs、t、s将显示在workspace列表中。
在命令窗口键入Sptool，将弹出Sptool主界面，点击菜单将信号xn导入并取名为sig1。
（2）单击Filters列表下的New，按照参数要求设计出kaiser窗型带通滤波器filt1
（3）将滤波器filt1应用到sig1信号序列。
分别在Signals、Filters、Spectra列表中选择sig1、filt1、mtlbse，单击Filters列表下的Apply按钮，在弹出的Apply Filter对话框中将输出信号命名为sig2。
（4）进行频谱分析。
在Signals中选择sig1，单击Spectra下的Create按钮，在弹出
Spectra Viewer界面中选择Method为FFT，Nfft=256，单击Apply按钮生成sig1的频谱spect1。
同样的步骤可以生成信号sig2的频谱spect2。
分别选中信号sig1、sig2、spect1、spect2，单击各自列表下方的View按钮，即可观察他们的波形。
按住shift再点击两个需要观察的信号，然后点击view，即可同时观察两个信号的图形。
下图中蓝色是滤波前的信号，红色是滤波后的信号。
频域图像：
时域图像：
升级会员 