数字信号课程设计.docx
《数字信号课程设计.docx》由会员分享，可在线阅读，更多相关《数字信号课程设计.docx（14页珍藏版）》请在冰点文库上搜索。
数字信号课程设计
2010秋季学期数字信号处理课程设计
IIR低/高通、带通/阻数字滤波器的MATLAB实现
输入信号：
频率为59Hz的周期方波x（t）的抽样序列x[k],其频谱为X（m）。
数字滤波器的技术指标：
低通截止频率为0.1π，高通截止频率为0.5π，带通通带截止频率为0.1π和0.5π，带阻阻带截止频率为0.2π和0.3π，通带衰减最大1dB，阻带衰减最小50dB。
数字滤波器的单位冲击响应：
h[k]，其幅度谱为H（Ω）。
输出信号：
y[k]=x[k]*h[k]，则其频谱为Y（m）。
一、MATLAB实现程序
I脉冲响应不变法设计BW型低通数字滤波器，指标：
Wp=0.1*pi;Ws=0.55*pi;Ap=1;As=50;
程序如下：
%方波抽样序列及其频谱
f=59;N=80;T=1/f;dt=T/N;n=0:
N-1;tn=n*dt;
x=square（2*pi*f*tn,50）;
L=256;
X=fftshift（fft（x,L））;
w=linspace（-pi,pi,L）/pi;
%脉冲响应不变法设计BW型LPDF
%数字滤波器指标
Wp=0.1*pi;Ws=0.55*pi;Ap=1;As=50;Fs=1;
%确定原型低通滤波器指标
wp=Wp*Fs;ws=Ws*Fs;
%设计滤波器
[N,wc]=buttord（wp,ws,Ap,As,'s'）
[num,den]=butter（N,wc,'s'）;
[numd,dend]=impinvar（num,den,Fs）;
%幅度频率响应曲线
omega1=linspace（0,Wp,500）;
omega2=linspace（Wp,Ws,200）;
omega3=linspace（Ws,pi,500）;
H1=20*log10（abs（freqz（numd,dend,omega1）））;
H2=20*log10（abs（freqz（numd,dend,omega2）））;
H3=20*log10（abs（freqz（numd,dend,omega3）））;
%与给定的指标对比
fprintf（'Ap=%.4f\n',max（-H1））;
fprintf（'As=%.4f\n',max（-H3））;
%获取冲击响应h（k）
[h,k]=impz（numd,dend）;
%获取输出信号及其频谱
y=conv（x,h）;
Y=fftshift（fft（y,L））;
%画图
subplot（3,2,1）,stem（n,x,'filled'）;
xlabel（'k'）;ylabel（'x[k]'）;
title（'输入一个周期方波序列'）;
subplot（3,2,2）,plot（w,abs（X））;
xlabel（'w（pirad）'）;ylabel（'X（m）'）;
title（'输入一个周期方波序列频谱'）;
subplot（3,2,3）,plot（k,h）;
xlabel（'k'）;ylabel（'h[k]'）;
title（'数字滤波器的单位冲击响应序列'）;
subplot（3,2,4）,plot（[omega1omega2omega3]/pi,[H1H2H3]）;
ylabel（'幅度dB'）;
xlabel（'w（pirad）'）;
title（'数字滤波器幅度谱'）;axis（[0,1,-60,10]）;
gridon
subplot（3,2,5）,plot（y）;
xlabel（'k'）;ylabel（'y[k]'）;
title（'输出序列'）;
subplot（3,2,6）,plot（w,abs（Y））;
xlabel（'w（pirad）'）;ylabel（'Y（m）'）;
title（'输出序列频谱'）;
II双线性变换法设计BW型高通数字滤波器，指标：
Wp=0.5*pi;Ws=0.12*pi;Ap=1;As=50;
程序如下：
%方波抽样序列及其频谱
f=59;N=80;T=1/f;dt=T/N;n=0:
N-1;tn=n*dt;
x=square（2*pi*f*tn,50）;
L=256;
X=fftshift（fft（x,L））;
w=linspace（-pi,pi,L）/pi;
%双线性变换法设计BW型HPDF
%数字滤波器指标
Wp=0.5*pi;Ws=0.12*pi;Ap=1;As=50;Ts=2;Fs=1/Ts;
%转换成模拟滤波器指标
wph=tan（Wp/2）;wsh=tan（Ws/2）;
%确定原型低通滤波器指标
wp=1/wph;ws=1/wsh;
%设计滤波器
[N,wc]=buttord（wp,ws,Ap,As,'s'）
[num,den]=butter（N,wc,'s'）;
[numt,dent]=lp2hp（num,den,1）;
[numd,dend]=bilinear（numt,dent,Fs）;
%幅度频率响应曲线
omega1=linspace（0,Ws,500）;
omega2=linspace（Ws,Wp,200）;
omega3=linspace（Wp,pi,500）;
H1=20*log10（abs（freqz（numd,dend,omega1）））;
H2=20*log10（abs（freqz（numd,dend,omega2）））;
H3=20*log10（abs（freqz（numd,dend,omega3）））;
%与给定的指标对比
fprintf（'As=%.4f\n',min（-H1））;
fprintf（'Ap=%.4f\n',min（-H2））;
%获取冲击响应h（k）
[h,k]=impz（numd,dend）;
%获取输出信号及其频谱
y=conv（x,h）;
Y=fftshift（fft（y,L））;
%画图
subplot（3,2,1）,stem（n,x,'filled'）;
xlabel（'k'）;ylabel（'x[k]'）;
title（'输入一个周期方波序列'）;
subplot（3,2,2）,plot（w,abs（X））;
xlabel（'w（pirad）'）;ylabel（'X（m）'）;
title（'输入一个周期方波序列频谱'）;
subplot（3,2,3）,plot（k,h）;
xlabel（'k'）;ylabel（'h[k]'）;
title（'数字滤波器的单位冲击响应序列'）;
subplot（3,2,4）,plot（[omega1omega2omega3]/pi,[H1H2H3]）;
ylabel（'幅度dB'）;
xlabel（'w（pirad）'）;
title（'数字滤波器幅度谱'）;axis（[0,1,-60,10]）;
gridon
subplot（3,2,5）,plot（y）;
xlabel（'k'）;ylabel（'y[k]'）;
title（'输出序列'）;
subplot（3,2,6）,plot（w,abs（Y））;
xlabel（'w（pirad）'）;ylabel（'Y（m）'）;
title（'输出序列频谱'）;
III双线性变换法设计BW型带通数字滤波器，指标：
Wp1=0.1*pi;Wp2=0.5*pi;Ws1=0.04*pi;Ws2=0.75*pi;Ap=1;As=50;
程序如下：
%方波抽样序列及其频谱
f=59;N=80;T=1/f;dt=T/N;n=0:
N-1;tn=n*dt;
x=square（2*pi*f*tn,50）;
L=256;
X=fftshift（fft（x,L））;
w=linspace（-pi,pi,L）/pi;
%双线性变换法设计BW型BPDF
%数字滤波器指标
Wp1=0.1*pi;Wp2=0.5*pi;Ws1=0.04*pi;Ws2=0.75*pi;Ap=1;As=50;
%转换成模拟滤波器指标
Ts=2;Fs=1/Ts;
wp1=tan（Wp1/2）;wp2=tan（Wp2/2）;ws1=tan（Ws1/2）;ws2=tan（Ws2/2）;
B=wp2-wp1;w0=（wp1*wp2）^0.5;
%确定原型低通滤波器指标
bpws1=（ws1^2-w0^2）/（B*ws1）;bpws2=（ws2^2-w0^2）/（B*ws2）;
wp=1;ws=min（abs（bpws1）,abs（bpws2））;
%设计滤波器
[N,wc]=buttord（wp,ws,Ap,As,'s'）
[num,den]=butter（N,wc,'s'）;
[numt,dent]=lp2bp（num,den,w0,B）;
[numd,dend]=bilinear（numt,dent,Fs）;
%幅度频率响应曲线
omega1=linspace（0,Ws1,500）;
omega2=linspace（Ws1,Wp1,200）;
omega3=linspace（Wp1,Wp2,500）;
omega4=linspace（Wp2,Ws2,500）;
omega5=linspace（Ws2,pi,200）;
H1=20*log10（abs（freqz（numd,dend,omega1）））;
H2=20*log10（abs（freqz（numd,dend,omega2）））;
H3=20*log10（abs（freqz（numd,dend,omega3）））;
H4=20*log10（abs（freqz（numd,dend,omega4）））;
H5=20*log10（abs（freqz（numd,dend,omega5）））;
%与给定的指标对比
fprintf（'As1=%.4f\n',min（-H1））;
fprintf（'Ap1=%.4f\n',min（-H2））;
fprintf（'Ap2=%.4f\n',min（-H4））;
fprintf（'As2=%.4f\n',max（-H4））;
%获取冲击响应h（k）
[h,k]=impz（numd,dend）;
%获取输出信号及其频谱
y=conv（x,h）;
Y=fftshift（fft（y,L））;
%画图
subplot（3,2,1）,stem（n,x,'filled'）;
xlabel（'k'）;ylabel（'x[k]'）;
title（'输入一个周期方波序列'）;
subplot（3,2,2）,plot（w,abs（X））;
xlabel（'w（pirad）'）;ylabel（'X（m）'）;
title（'输入一个周期方波序列频谱'）;
subplot（3,2,3）,plot（k,h）;
xlabel（'k'）;ylabel（'h[k]'）;
title（'数字滤波器的单位冲击响应序列'）;
subplot（3,2,4）,plot（[omega1omega2omega3omega4omega5]/pi,[H1H2H3H4H5]）;
ylabel（'幅度dB'）;
xlabel（'w（pirad）'）;
title（'数字滤波器幅度谱'）;axis（[0,1,-60,10]）;
gridon
subplot（3,2,5）,plot（y）;
xlabel（'k'）;ylabel（'y[k]'）;
title（'输出序列'）;
subplot（3,2,6）,plot（w,abs（Y））;
xlabel（'w（pirad）'）;ylabel（'Y（m）'）;
title（'输出序列频谱'）;
IV双线性变换法设计BW型带阻数字滤波器，指标：
Wp1=0.12*pi;Wp2=0.48*pi;Ws1=0.2*pi;Ws2=0.3*pi;Ap=1;As=50;
程序如下：
%方波抽样序列及其频谱
f=59;N=80;T=1/f;dt=T/N;n=0:
N-1;tn=n*dt;
x=square（2*pi*f*tn,50）;
L=256;
X=fftshift（fft（x,L））;
w=linspace（-pi,pi,L）/pi;
%双线性变换法设计BW型BSDF
%数字滤波器指标
Wp1=0.12*pi;Wp2=0.48*pi;Ws1=0.2*pi;Ws2=0.3*pi;Ap=1;As=50;
%转换成模拟滤波器指标
Ts=2;Fs=1/Ts;
wp1=tan（Wp1/2）;wp2=tan（Wp2/2）;ws1=tan（Ws1/2）;ws2=tan（Ws2/2）;
B=ws2-ws1;w0=（ws1*ws2）^0.5;
%确定原型低通滤波器指标
bpwp1=（B*wp1）/（-wp1^2+w0^2）;bpwp2=（B*wp2）/（-wp2^2+w0^2）;
ws=1;wp=min（abs（bpwp1）,abs（bpwp2））;
%设计滤波器
[N,wc]=buttord（wp,ws,Ap,As,'s'）
[num,den]=butter（N,wc,'s'）;
[numt,dent]=lp2bs（num,den,w0,B）;
[numd,dend]=bilinear（numt,dent,Fs）;
%幅度频率响应曲线
omega1=linspace（0,Wp1,200）;
omega2=linspace（Wp1,Ws1,500）;
omega3=linspace（Ws1,Ws2,500）;
omega4=linspace（Ws2,Wp2,500）;
omega5=linspace（Wp2,pi,200）;
H1=20*log10（abs（freqz（numd,dend,omega1）））;
H2=20*log10（abs（freqz（numd,dend,omega2）））;
H3=20*log10（abs（freqz（numd,dend,omega3）））;
H4=20*log10（abs（freqz（numd,dend,omega4）））;
H5=20*log10（abs（freqz（numd,dend,omega5）））;
%与给定的指标对比
fprintf（'Ap1=%.4f\n',max（-H1））;
fprintf（'As1=%.4f\n',max（-H2））;
fprintf（'As2=%.4f\n',max（-H4））;
fprintf（'Ap2=%.4f\n',min（-H4））;
%获取冲击响应h（k）
[h,k]=impz（numd,dend）;
%获取输出信号及其频谱
y=conv（x,h）;
Y=fftshift（fft（y,L））;
%画图
subplot（3,2,1）,stem（n,x,'filled'）;
xlabel（'k'）;ylabel（'x[k]'）;
title（'输入一个周期方波序列'）;
subplot（3,2,2）,plot（w,abs（X））;
xlabel（'w（pirad）'）;ylabel（'X（m）'）;
title（'输入一个周期方波序列频谱'）;
subplot（3,2,3）,plot（k,h）;
xlabel（'k'）;ylabel（'h[k]'）;
title（'数字滤波器的单位冲击响应序列'）;
subplot（3,2,4）,plot（[omega1omega2omega3omega4omega5]/pi,[H1H2H3H4H5]）;
ylabel（'幅度dB'）;
xlabel（'w（pirad）'）;
title（'数字滤波器幅度谱'）;axis（[0,1,-60,10]）;
gridon
subplot（3,2,5）,plot（y）;
xlabel（'k'）;ylabel（'y[k]'）;
title（'输出序列'）;
subplot（3,2,6）,plot（w,abs（Y））;
xlabel（'w（pirad）'）;ylabel（'Y（m）'）;
title（'输出序列频谱'）;
二、仿真结果图形
I低通数字滤波器
II高通数字滤波器
III带通数字滤波器
IV带阻数字滤波器
三、结果分析及总结
分析：
题目中没有给出低通、高通、带通滤波器的阻带截频以及带阻滤波器的通带截频，根据题目要求我们可以有很多种选择，而根据选择不同，所得波形的清晰度也不太相同。
考虑到实现滤波器的成本问题，我选择了能使滤波器在达到要求的设计指标时的阶数尽可能低的参数。
另外，我在设计低通滤波器时采用了脉冲响应不变法，这样得到的波形更稳定，又发现设计题中所给定指标的高通，带通，带阻滤波器用脉冲响应不变法设计时，所得到的滤波器不稳定，故采用了双线性变换法设计。
我通过设计滤波器了解了MATLAB的运行环境,从图形和程序两方面了解了四种滤波器的不同,让我在理论知识的了解方面更轻松.
总结：
以下是各种滤波器的阶数和所设计滤波器的通带衰减Ap和阻带衰减As
I低通
N=4
wc=0.4097
Ap=0.4902
As=65.1339
可以看出通带衰减和阻带衰减均符合设计要求。
II高通
N=4
wc=1.2431
As=50.0000
Ap=0.7017
可以看出通带衰减和阻带衰减均符合设计要求。
III带通
N=7
wc=1.2262
As1=52.6822
Ap1=0.2431
Ap2=0.2431
As2=50.0000
可以看出通带衰减和阻带衰减均符合设计要求。
IV带阻
N=5
wc=
0.3162
Ap1=0.8881
As1=50.0000
As2=50.0000
Ap2=0.2896
可以看出通带衰减和阻带衰减均符合设计要求。
升级会员 