基于频率采样法FIR带通滤波器设计

**大学电子信息工程学院DSP课程设计报告题目：基于频率采样法的FIR带通滤波器的设计专业班级：通信工程组员：***指导老师：***2011年6月13日目录一、设计目的………………………………………………3二、设计要求………………………………………………3三、设计原理………………………………………………4四、程序代码………………………………………………6五、调试分析………………………………………………12六、设计总结………………………………………………13七、设计心得………………………………………………13八、参考文献………………………………………………一、设计目的：1、掌握用频率取样法设计FIRDF的方法，并掌握该方法的MATLAB编程。2、熟悉频率取样理论，熟悉内插函数及其应用。3、了解FIRDF的频率特性和相位特性，观察过渡带取样点对滤波器幅频特性的影响。二、设计内容:基于MATLAB结合FFT和IFFT，利用频率采样法设计FIR数字带通滤波器，然后用自己设计的滤波器对采集的加噪后的语音信号进行滤波，并将滤波前后的信号进行比较，回放语音信号。三、设计原理：FIR滤波器的单位脉冲响应nh是有限长的（10Nn），其z变换为z1的（N-1）阶多项式：znNnnhzXzYzH10可得FIR滤波器的系统差分方程：nxnbmnxmbNnxNbnxbnxbnyNm1011110因此，FIR滤波器又称为卷积滤波器。FIR滤波器的频率响应表达式为：eejwnNnjwnhH10信号通过FIR滤波器不失真条件是在通带内具有恒定的幅频特性和相位线性特性。理论上可以证明：当FIR滤波器的系数满足下列中心对称条件：nNhnhnNhnh11或者时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个N阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频域内信号通过滤波器后仍保持原有波形形状而无相位失真。FIR滤波器设计的方法很多，如窗函数法、频率采样法以及其他的各种优化设计方法，本次课程设计使用频率采样法设计FIR带通滤波器。频率采样法是从频域出发。因为有限长序列h(n)又可用其离散傅立叶变换H(k)来唯一表示，H(k)与所要求的FIR滤波器系统函数Hd(z)之间存在着频率取样关系。即Hd(z)在Z平面单位圆上按角度等分的取样值等于Hd(k)的各相应值，就以此Hd(k)值作为实际FIR数字滤波器频率特性的取样值H(k),或者说H(k)正是所要求的频率响应H(ejw)的N各等间隔的取样值。频率采样法就是根据频域采样理论，由滤波特性指标构造希望逼近的滤波器频响函数Hd(ejω)，对其在[0,2π]上采样得到:1102NkeHkHkNjdd,,,然后，就可求出单位脉冲响应h(n)，或是系统函数H(z)。这样，h(n)或是H(z)就是滤波器的设计结果。频率取样法设计的基本思想：把给出的理想频率响应进行取样，通过IDFT从频谱样点直接求得有限脉冲响应。其设计过程如下频率取样法的关键是正确确定数字频域系统函数H(k)在Ω∈[0，2π]内的N个样点，其约束条件为H(k)=H(N-k)ϕ（m）=-ϕ(N-m)0≤k≤N-1频率采样法的优点是可以在频域直接设计，并且适合最优化设计；缺点是采样频率只能等于N/2的整数倍，因而不能确保截止频率wc的自由取值，要想实现自由地选择截止频率，必须增加采样点数N，但是这又使计算量加大。内插公式1,,1,0111,,1,0101NkzWkHNzzHNnkHIDFTnhNkkNN)()()()(H2jDTFTIDTFTNkjdeHnhkHek四、程序代码:语音程序：filename='111';[s,fs,nbits]=wavread(filename);sound(s,fs,nbits);%回放语音信号n=length(s);%求出语音信号的长度t=0:1/fs:(n-1)/fs;Y=fft(s);%傅里叶变换figure(1)subplot(2,1,1);plot(s);title('原始信号波形');gridonsubplot(2,1,2);plot(abs(Y));title('原始信号频谱')sound(s,fs);gridon仿真结果：图一带通滤波器程序：N=40;alfa=(40-1)/2;k=0:N-1;w1=(2*pi/N)*k;T1=0.109021;T2=0.59417456;hrs=[zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4)];hdr=[0,0,1,1,0,0];wd1=[0,0.2,0.35,0.65,0.8,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alfa*(2*pi)/N*k1,alfa*(2*pi/N*(N-k2))];H=hrs.*exp(j*angH);h=real(ifft(H));[db,mag,pha,grd,w]=freqz_m(h,1);[Hr,ww,a,L]=Hr_Type2(h);subplot(2,2,1)plot(w1(1:21)/pi,hrs(1:21),'o',wd1,hdr)axis([0,1,-0.1,1.1]);title('带通：N=40，T1=0.109021，T2=0.59417456')ylabel('Hr(k)');set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1])set(gca,'YTickMode','manual','YTick',[0,0.059,0.109,1]);grid%绘制带网格的图像subplot(2,2,2);stem(k,h);axis([-1,N,-0.4,0.4])title('脉冲响应');ylabel('h(n)');text(N+1,-0.4,'n')subplot(2,2,3);plot(ww/pi,Hr,w1(1:21)/pi,hrs(1:21),'o');axis([0,1,-0.1,1.1]);title('振幅响应')xlabel('频率（单位：pi）');ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1]);set(gca,'YTickMode','manual','YTick',[0,0.059,0.109,1]);gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-100,10]);gridtitle('幅度响应');xlabel('频率（单位：pi）');ylabel('分贝')set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1])set(gca,'YTickMode','manual','YTick',[-60;0]);set(gca,'YTickLabelMode','manual','YTickLabels',[60;0]);仿真结果：图二总程序：%设计一个通带为0.35π到0.65π，阻带下截止为0.2π，上截止为0.8πClear;N=40;alfa=(40-1)/2;k=0:N-1;w1=(2*pi/N)*k;T1=0.109021;T2=0.59417456;hrs=[zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4)];hdr=[0,0,1,1,0,0];wd1=[0,0.2,0.35,0.65,0.8,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alfa*(2*pi)/N*k1,alfa*(2*pi/N*(N-k2))];H=hrs.*exp(j*angH);h=real(ifft(H));[db,mag,pha,grd,w]=freqz_m(h,1);[Hr,ww,a,L]=Hr_Type2(h);figure(1)subplot(2,2,1)plot(w1(1:21)/pi,hrs(1:21),'o',wd1,hdr)axis([0,1,-0.1,1.1]);title('带通：N=40，T1=0.109021，T2=0.59417456')ylabel('Hr(k)');set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1])set(gca,'YTickMode','manual','YTick',[0,0.059,0.109,1]);grid%绘制带网格的图像subplot(2,2,2);stem(k,h);axis([-1,N,-0.4,0.4])title('脉冲响应');ylabel('h(n)');text(N+1,-0.4,'n')subplot(2,2,3);plot(ww/pi,Hr,w1(1:21)/pi,hrs(1:21),'o');axis([0,1,-0.1,1.1]);title('振幅响应')xlabel('频率（单位：pi）');ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1]);set(gca,'YTickMode','manual','YTick',[0,0.059,0.109,1]);gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-100,10]);gridtitle('幅度响应');xlabel('频率（单位：pi）');ylabel('分贝')set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1])set(gca,'YTickMode','manual','YTick',[-60;0]);set(gca,'YTickLabelMode','manual','YTickLabels',[60;0]);[s,fs,nbits]=wavread('sj.wav');%信号de取样频率为44100HZx=s(:,1);sound(x,fs);L=length(x);f=fs*(0:L-1)/L;t=0:1/fs:(L-1)/fs;%将所加噪声信号的点数调整到与原始信号相同d=0.03*abs(max(x))*[cos(2*pi*22000*t)]';%噪声为500和3300Hz的余弦信号dz=[cos(2*pi*11025*t)]';xd=x.*dz;xz=xd+d;sound(xz,fs);%播放加噪声后的语音信号X=fft(x);%求信号的频谱XD=fft(xd);%信号调制后的频谱XZ=fft(xz);figure(2)subplot(3,1,1);plot(t,x)title('未加噪的信号');xlabel('times');ylabel('幅度');subplot(3,1,2);plot(t,xd)title('调制后的信号');xlabel('times');ylabel('幅度');subplot(3,1,3);plot(t,xz)tit