数字信号处理实验报告 完整版

./实验1利用DFT分析信号频谱一、实验目的1.加深对DFT原理的理解。2.应用DFT分析信号的频谱。3.深刻理解利用DFT分析信号频谱的原理,分析实现过程中出现的现象及解决方法。二、实验设备与环境计算机、MATLAB软件环境三、实验基础理论1.DFT与DTFT的关系有限长序列的离散时间傅里叶变换在频率区间的N个等间隔分布的点上的N个取样值可以由下式表示：由上式可知,序列的N点DFT,实际上就是序列的DTFT在N个等间隔频率点上样本QUOTE。2.利用DFT求DTFT方法1：由QUOTE恢复出QUOTE的方法如下：IDFTDTFTIDFTDTFT𝑋<𝑒𝑗𝜔>QUOTE由图2.1所示流程可知：由上式可以得到：其中QUOTE为插函数方法2：实际在MATLAB计算中,上述插值运算不见得是最好的办法。由于DFT是DTFT的取样值,其相邻两个频率样本点的间距为2π/N,所以如果我们增加数据的长度N,使得到的DFT谱线就更加精细,其包络就越接近DTFT的结果,这样就可以利用DFT计算DTFT。如果没有更多的数据,可以通过补零来增加数据长度。3.利用DFT分析连续信号的频谱采用计算机分析连续时间信号的频谱,第一步就是把连续信号离散化,这里需要进行两个操作：一是采样,二是截断。对于连续时间非周期信号QUOTE,按采样间隔T进行采样,阶段长度M,那么：对QUOTE进行N点频域采样,得到因此,可以将利用DFT分析连续非周期信号频谱的步骤归纳如下：〔1确定时域采样间隔T,得到离散序列〔2确定截取长度M,得到M点离散序列,这里为窗函数。〔3确定频域采样点数N,要求N≥M。〔4利用FFT计算离散序列的N点DFT,得到.〔5根据上式由计算采样点QUOTE的近似值。采用上述方法计算信号QUOTE的频谱需要注意如下三个问题：〔1频谱混叠。如果不满足采样定理的条件,频谱会出现混叠误差。对于频谱无限宽的信号,应考虑覆盖大部分主要频率分量的围。〔2栅栏效应和频谱分辨率。使用DFT计算频谱,得到的结果只是N个频谱样本值,样本值之间的频谱是未知的,像通过一个栅栏观察频谱,称为"栅栏效应"。频谱分辨率与记录长度成反比,要提高频谱分辨率,就要增加记录时间。〔3频谱泄露。对信号截断会把窗函数的频谱引入信号频谱,造成频谱泄露。解决这个问题的主要办法是采用旁瓣小的窗函数,频谱泄露和窗函数均会引起误差。因此,要合理选取采样间隔和截取长度,必要时还需考虑加适当的窗。对于连续时间周期信号,我们在采用计算机进行计算时,也总是要进行截断,序列总是有限长的,仍然可以采用上述方法近似计算。4.可能用到的MATLAB函数与代码实验中DFT运算可采用MATLAB中提供的函数fft来实现。DTFT可采用MATLAB矩阵运算的方法进行计算,如下式所示：四、实验容1、已知,完成如下要求：〔1计算其DTFT,并画出区间的波形。〔2计算4点DFT,并把结果显示在〔1所画的图形中。〔3对x<n>补零,计算64点DFT,并显示结果。〔4根据实验结果,分析是否可以由DFT计算DTFT,如果可以,如何实现。解：〔1计算其DTFT,并画出区间的波形。>>n=0:3;>>x=[2-111];>>w=-pi:0.01*pi:pi;>>X=x*exp<-j*n'*w>;>>subplot<211>;>>plot<w,abs<X>>;>>xlabel<'\Omega/\pi'>;>>title<'Magnitude'>;>>axistight;>>subplot<212>;>>plot<w,angle<X>/pi>;>>xlabel<'\Omega/\pi'>;>>title<'Phase'>;>>axistight;〔2计算4点DFT,并把结果显示在〔1所画的图形中。>>Xk=fft<x>;>>subplot<211>;>>holdon;>>stem<n,abs<Xk>,'filled'>;>>plot<w,abs<X>>;>>axistight;>>xlabel<'\Omega/\pi'>;>>title<'Magnitude'>;>>subplot<212>;>>holdon;>>plot<w,angle<X>/pi>;>>stem<n,angle<Xk>,'filled'>;>>axistight;>>xlabel<'\Omega/\pi'>;>>title<'Phase'>;运行结果如下：〔3对x<n>补零,计算64点DFT,并显示结果。>>x=[2-111zeros<1,60>];>>n=0:63;>>Xk=fft<x>;>>subplot<211>;>>stem<n,abs<Xk>,'filled'>;>>subplot<212>;>>stem<n,angle<Xk>,'filled'>;〔4根据实验结果,分析是否可以由DFT计算DTFT,如果可以,如何实现。 可以由DFT计算DTFT,序列补零后,长度越长,DFT点越多,其DFT越逼近DTFT连续波形。所以令序列补零至足够长时其DFT序列的波形与DTFT的波形在一定的分辨率下已经相同。2、考察序列x<n>=cos<0.48πn>+cos<0.52πn>〔10<=n<=10时,用DFT估计x<n>的频谱；将x<n>补零加长到长度为100点序列用DFT估计x<n>的频谱,要求画出相应波形。〔20<=n<=100时,用DFT估计x<n>的频谱。并画出波形。〔3根据实验结果,分析怎样提高频谱分辨率解：〔10<=n<=10时,用DFT估计x<n>的频谱；将x<n>补零加长到长度为100点序列用DFT估计x<n>的频谱,要求画出相应波形。>>n=0:10;>>x=cos<0.48*pi.*n>+cos<0.52*pi.*n>;>>Xk=fft<x>;>>subplot<211>;>>stem<n,Xk,'filled'>;>>x=[cos<0.48*pi.*n>+cos<0.52*pi.*n>zeros<1,89>];>>Xk=fft<x>;>>subplot<212>;>>stem<Xk,'filled'>;.〔20<=n<=100时,用DFT估计x<n>的频谱。并画出波形。程序代码如下：>>n=0:100;>>x=cos<0.48*pi.*n>+cos<0.52*pi.*n>;>>Xk=fft<x>;>>stem<Xk,'filled'>;>>axistight;〔3根据实验结果,分析怎样提高频谱分辨率可以通过如下三种方式来增加分辨率。a、增加时域信号采样时间b、提高采样频率c、补零3、已知信号x<t>=0.15sin<2πf1t>+sin<2πf2t>-0.1sin<2πf3t>,其f1=1Hz,f2=2Hz,f3=3Hz。从x<t>的表达式可以看出,它包含三个频率的正弦波,但是,从其时域波形来看,似乎是一个正弦信号,利用DFT做频谱分析,确定适合的参数,使得到的频谱的频率分辨率符合需要。n=0:10;x=0.15*sin<2*pi.*n>+sin<4*pi.*n>-0.1*sin<6*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;n=0:100;x=0.15*sin<0.2*pi.*n>+sin<0.4*pi.*n>-0.1*sin<0.6*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;n=0:200;x=0.15*sin<0.1*pi.*n>+sin<0.2*pi.*n>-0.1*sin<0.3*pi.*n>;Xk=fft<x>;stem<abs<Xk>,'filled'>;结果分析：上图为x<t>信号截取过后的连续时间信号的傅里叶变换幅频特性曲线,截取周期为100s〔即为采样时间M,根据的特点,知1hz处的幅值为,,与图像相符。4、利用DFT近似分析连续时间信号x<t>=e-0.1u<t>的频谱〔幅度值。分析采用不同的采样间隔和截取长度进行计算的结果,并最终确定合适的参数。n=0:10;x=exp<-0.1.*n>;Xk=fft<x>;stem<Xk,'filled'>;n=0:20;x=exp<-0.05.*n>;Xk=fft<x>;stem<Xk,'filled'>;n=0:40;x=exp<-0.025.*n>;Xk=fft<x>;stem<Xk,'filled'>;五、心得体会通过本次实验,加深了对DFT原理的理解。学会了应用DFT分析信号的频谱。深刻理解到利用DFT分析信号频谱的原理,能够分析实现过程中出现的现象及解决方法。实验二利用FFT计算线性卷积一、实验目的1.掌握利用FFT计算线性卷积的原理及具体实现方法。2.加深理解重叠相加法和重叠保留法。3.考察利用FFT计算线性卷积各种方法的适用围。二、实验基础理论1.线性卷积与圆周卷积设x<n>为L点序列,h<n>为M点序列,x<n>和h<n>的线性卷积为的长度为L+M-1。x<n>和h<n>的圆周卷积为圆周卷积与线性卷积相等而不产生交叠的必要条件为N圆周卷积定理：根据DFT性质,x<n>和h<n>的N点圆周卷积的DFT等于它们的DFT的乘积：2.快速卷积快速卷积发运用圆周卷积实现线性卷积,根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：<1>必须选择；为了能使用基-2算法,要求。采用补零的办法使得x<n>和h<n>的长度均为N。<2>计算x<n>和h<n>的N点FFT。<3>组成乘积<3>利用IFFT计算Y<k>的IDFT,得到线性卷积y<n>3.分段卷积我们考察单位取样响应为h<n>的线性系统,输入为x<n>,输出为y<n>,则如果x<n>极长时,如果要等x<n>全部集齐时再开始进行卷积,会使输出有较大延时；如果序列太长,需要大量存储单元。为此,我们把x<n>分段,为别求出每段的卷积,合在一起得到最后的总输出。这称为分段卷积。分段卷积可以细分为重叠保留法和重叠相加法。重叠保留法：设x<n>的长度为,h<n>的长度为M。把序列x<n>分成多段N点序列,每段雨前一段重写M-1个样本。并在第一个输入段前面补M-1个零。计算每一段与h<n>的圆周卷积,其结果中前M-1个不等与线性卷积,应当舍去,只保留后面N-M+1个正确的输出样本,把它们合起来得到总的输出。利用FFT实现重叠保留法的步骤如下：<1>在x<n>前面填充M-1个零,扩大以后的序列为<2>将x<n>分为若干段N点子段,设L=N-M+1为每一段的有效长度,则第i段的数据为：<3>计算每一段与h<n>的N点圆周卷积,利用FFT计算圆周卷积<4>设每一段卷积结果的前M-1个样本,连接剩下的样本得到卷积结果y<n>。重叠相加法：设h<n>长度为M,将信号x<n>分解成长为L的子段。以表示没断信号,则：每一段卷积的长度为L+M-1,所以在做求和时,相邻两段序列由M-1个样本重叠,即前一段的最后M-1个样本和下一段前M-1个样本序列重叠,这个重叠部分相加,再与不重叠的部分共同组成y<n>。利用FFT实现重叠保留法的步骤如下：<1>将x<n>分为若干L点子段。<2>计算每一段与h<n>的卷积,根据快速卷积法利用FFT计算卷积。<3>将各段相加,得到输出y<n>。三、实验容假设要计算序列和的线性卷积,完成以下实验容：设L=M,根据线性卷积的表达式和快速卷积的原理,分别编程实现计算两个序列线性卷积的方法,比较当序列长度分别为8,16,32,64,256,512,1024时,两种计算方法计算线性卷积所需时间。>>fori=1:7L=input<'L:'>;n=0:L;x=heaviside<n>-heaviside<n-L>;h=cos<0.2*pi.*n>;ticy=conv<x,h>;tocendL:8时间已过0.000050秒。L:16时间已过0.000037秒。L:32时间已过0.000056秒。L:64时间已过0.000065秒。L:256时间已过0.000101秒。L:512时间已过0.000188秒。L:1024时间已过0.000254秒。>>fori=1:7L=input<'L:'>;n=0:L;x=heaviside<n>-heaviside<n-L>;h=cos<0.2*pi.*n>;ticXk=fft<x,L+1>;Hk=fft<h,L+1>;Yk=Xk.*Hk;y=ifft<Yk>;tocendL:8时间已过0.000048秒。L:16时间已过0.000044秒。L:32时间已过0.000042秒。L:64时间已过0.000055秒。L:256时间已过0.000134秒。L:512时间已过0.000147秒。L:1024时间已过0.000176秒。当L=2048且M=256时,比较计算线性卷积和快速卷积所需的时间,进一步考察当L=4096且M=256时两种算法所需时间。>>fori=1:2L=input<'L:'>;M=input<'M:'>;n1=0:L;n2=0:M;x=heaviside<n1>-heaviside<n1-L>;h=cos<0.2*pi.*n2>;ticXk=fft<x,L+1>;Hk=fft<h,L+1>;Yk=Xk.*Hk;y=ifft<Yk>;tocendL:2048M:256时间已过0.001139秒。L:4096M:256时间已过0.001503秒。>>fori=1:2L=input<'L:'>;M=input<'M:'>;n1=0:L;n2=0:M;x=heaviside<n1>-heaviside<n1-L>;h=cos<0.2*pi.*n2>;ticy=conv<x,h>;tocendL:2048M:256时间已过0.000263秒。L:4096M:256时间已过0.000368秒。3.>>L=input<'L:'>;M=input<'M:'>;n2=0:M;k=L/M;h=cos<0.2*pi.*n2>;Y=0;ticfori=1:kn1=<i-1>*M:i*M-1;x=heaviside<n1>-heaviside<n1-L>;Xk=fft<x,M>;Hk=fft<h,M>;Yk=Xk.*Hk;y=ifft<Yk,M>;Y=Y+y;endtocL:2048M:256时间已过0.025870秒。>>L=input<'L:'>;M=input<'M:'>;n2=0:M;k=L/M;h=cos<0.2*pi.*n2>;Y=0;ticfori=1:kn1=<i-1>*M:i*M-1;x=heaviside<n1>-heaviside<n1-L>;Xk=fft<x,M>;Hk=fft<h,M>;Yk=Xk.*Hk;y=ifft<Yk,M>;Y=Y+y;endtocL:4096M:256时间已过0.022300秒。L=4096;n1=0:L;x=heaviside<n1>-heaviside<n1-L>;n2=0:256;h=cos<0.2*pi.*n2>;ticy1=conv<x,h>;tocticXk=fft<x>;Hk=fft<h>;Yk=Xk.*Hk;y2=ifft<Yk>;toc编程实现利用重叠相加法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与第二题结果进行比较。clearticN=512;m=0:256;h=cos<0.2*pi*m>;n=0:2048;x=heaviside<n>-heaviside<n-2048>;Lenx=length<x>;M=length<h>;M1=M-1;L=N-M1;h=fft<h,N>;K=ceil<Lenx/L>;fori=Lenx:K*L-1x<i+1>=0;endY=zeros<K,N>;YY=zeros<1,<K-1>*L+N>;fork=0:K-1xk=[x<k*L+1:k*L+L>,zeros<1,M1>];Y<k+1,:>=<ifft<fft<xk>.*h>>;YY<k*L+1:k*L+N>=YY<k*L+1:k*L+N>+Y<k+1,:>;endtoc时间已过0.028816秒。编程实现利用重叠保留法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与第二题结果进行比较。clcclearticN=512;m=0:256;h=cos<0.2*pi*m>;n=0:2048;x=heaviside<n>-heaviside<n-2048>;Lenx=length<x>;M=length<h>;M1=M-1;L=N-M1;h=fft<h,N>;K=floor<<Lenx+M1-1>/L>+1;p=<K>*L-Lenx;x1=[zeros<1,M1>,x,zeros<1,p>];Y=zeros<K,N>;fork=0:K-1xk=fft<x1<k*L+1:k*L+N>>;Y<k+1,:>=<ifft<xk.*h>>;endZ=reshape<Y<:,M:N>',1,[]>;toc时间已过0.044424秒。四、实验心得通过本次实验,掌握了利用FFT计算线性卷积的原理及具体实现方法,加深了理解重叠相加法和重叠保留法。实验3IIR数字滤波器设计一、实验目的掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。加深理解数字滤波器和模拟滤波器之间的技术指标转化。掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用围。二、实验容设采样频率为10kHz,设计数字低通滤波器,满足如下指标 通带截止频率：1kHz,通带波动：1dB： 阻带截止频率：1.5kHz,阻带衰减：15dB：要求分别采用巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和椭圆模拟原型滤波器,并分别结合脉冲响应不变法和双线性变换法进行设计。结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求,分析脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用围。1巴特沃斯脉冲响应不变法：fp=1000fs=1500f=10000Wp=2*pi*fpWs=2*pi*fswp=2*f*tan<2*pi*fp/<2*f>>ws=2*f*tan<2*pi*fs/<2*f>>Rp=1As=15N1=ceil<<log10<<10^<Rp/10>-1>/<10^<As/10>-1>>>/<2*log10<Wp/Ws>>>N2=ceil<<log10<<10^<Rp/10>-1>/<10^<As/10>-1>>>/<2*log10<wp/ws>>>Omegac1=Wp/<<10^<Rp/10>-1>^<1/<2*N1>>>Omegac2=wp/<<10^<Rp/10>-1>^<1/<2*N1>>>[z,p,k]=buttap<N1>b1=k*Omegac1^N1a1=poly<p*Omegac1>[H,w]=freqs<b1,a1>subplot<221>plot<w/pi,abs<H>>gridonsubplot<223>plot<w/pi,angle<H>>gridon[b2,a2]=butter<N2,Omegac2,'S'>[H,w]=freqs<b2,a2>subplot<222>plot<w/pi,abs<H>>axistightgridonsubplot<224>plot<w/pi,angle<H>>gridon双线性法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<1/2>*<log10<<10^<As/10>-1>/<10^<Rp/10>-1>>>/<log10<ws1/wp1>>>wc=wp1/<10^<Rp/10>-1>^<1/2/N>[b1,a1]=butter<N,wc,'low','s'>;[b,a]=bilinear<b1,a1,F>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;2切比雪夫I型脉冲响应不变法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=wp/T;ws1=ws/T;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=wp/pi[b,a]=cheby1<N,Rp,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;双线性变换法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=wp/pi[b,a]=cheby1<N,Rp,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;3切比雪夫II型脉冲响应不变法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=wp/T;ws1=ws/T;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=ws/pi[b,a]=cheby2<N,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;双线性变换法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;N=ceil<<acosh<<A^2-1>^<1/2>/e>>/<acosh<ws1/wp1>>>wc=ws/pi[b,a]=cheby2<N,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;4椭圆原型模拟滤波器脉冲响应不变法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=wp/T;ws1=ws/T;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;k=wp1/ws1;k1=e/<A^2-1>^<1/2>;K1=ellipke<k>;K2=ellipke<<1-k1^2>^<1/2>>K3=ellipke<k1>K4=ellipke<<1-k^2>^<1/2>>;N=ceil<K1*K2/<K3*K4>>wc=wp/pi[b,a]=ellip<N,Rp,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;双线性法：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;wp1=2/T*tan<wp/2>;ws1=2/T*tan<ws/2>;e=<10^<Rp/10>-1>^<1/2>;A=10^<As/20>;k=wp1/ws1;k1=e/<A^2-1>^<1/2>;K1=ellipke<k>;K2=ellipke<<1-k1^2>^<1/2>>K3=ellipke<k1>K4=ellipke<<1-k^2>^<1/2>>;N=ceil<K1*K2/<K3*K4>>wc=wp/pi[b,a]=ellip<N,Rp,As,wc>;w=[0:500]*pi/500;[H,w]=freqz<b,a>;subplot<221>;plot<w/pi,abs<H>>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'|H<e^j^\Omega>|'>;subplot<222>;plot<w/pi,20*log10<<abs<H>>/max<abs<H>>>>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'|H<e^j^\Omega>|<dB>'>;subplot<223>;plot<w/pi,angle<H>/pi>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'PhaseofH<e^j^Omega><\pi>'>;grd=grpdelay<b,a,w>;subplot<224>;plot<w/pi,grd>;gridon;xlabel<'\Omega<\pi>'>;ylabel<'Groupdelay'>;三、结果分析：〔1对于设计要求的符合情况 根据以dB为单位的幅频特性图可以看出,能满足设计要求的滤波器有：脉冲响应不变法的切比雪夫Ⅰ型,以及双线性变换法的全部四种滤波器。而脉冲响应不变法的巴特沃斯、切比雪夫Ⅱ型和椭圆数字滤波器不满足设计要求。〔2脉冲响应不变法的优缺点及适用围优点：幅频特性与模拟滤波器线性对应。脉冲响应不变法的一个重要特点是频率坐标的变换是线性的,ω=ΩΤ,ω与Ω是线性关系。因此如果模拟滤波的频响带限于折叠频率以的话,通过变换后滤波器的频响可不失真地反映原响应与频率的关系。稳定性好。如果Ha<s>是稳定的,即其极点在S左半平面,映射到H<Z>也是稳定的。缺点：设计结果不固定。在设计时,往往不能直接根据要求的参数得到合适的滤波器,需要设计者多次修改设计参数。可能产生混叠失真。因频谱周期延拓效应,脉冲响应不变法只能用于带限的频响特性,如衰减特性很好的低通或带通。而高频衰减越大,频响的混淆效应越小,至于高通和带限滤波器,由于它们在高频部分不衰减,因此将完全混淆在低频响应中。所以用脉冲响应不变法实现高通和带限滤波器时,应增加一保护滤波器,滤掉高于折叠频率以上的频带,然后再用脉冲响应不变法转换为数字滤波器,这会增加设计的复杂性和滤波器的阶数。适用围：根据幅频特性与模拟滤波器相对应的特点,在要求时域脉冲响应能模仿模拟滤波器的场合,一般使用脉冲响应不变法。〔3双线性变换法的优缺点及适用围优点：消除了脉冲响应不变法所固有的混叠误差。靠频率的严重非线性关系得到S平面与Z平面的单值一一对应关系,整个jΩ轴单值对应于单位圆一周,其中ω和Ω为非线性关系。设计结果具有一般性。双线性变换法根据设计要求和计算流程计算得到的滤波器,一般能够满足原来的要求,得到的频率特性较好。双线性变换比脉冲响应法的设计计算更直接和简单。由于s与z之间的简单代数关系,所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。缺点：频率响应有畸变。Ω与ω的非线性关系,导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变,<使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变>。例如,一个模拟微分器,它的幅度与频率是线性关系,但通过双线性变换后,就不可能得到数字微分器。一个线性相位的模拟滤波器经双线性变换后,滤波器就不再有线性相位特性。适用围：双线性变换目前仍是使用得最普遍、最有成效的一种设计工具。这是因为大多数滤波器都具有分段常数的频响特性,如低通、高通、带通和带阻等,它们在通带要求逼近一个衰减为零的常数特性,在阻带部分要求逼近一个衰减为∞的常数特性,这种特性的滤波器通过双线性变换后,虽然频率发生了非线性变化,但其幅频特性仍保持分段常数的特性。四、实验心得通过本次实验,学会了掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。加深了理解数字滤波器和模拟滤波器之间的技术指标转化。掌握到脉冲响应不变法设计IIR数字滤波器的优缺点及适用围实验4FIR数字滤波器设计一、实验目的掌握窗函数法设计FIR数字滤波器的原理及具体方法。二、实验设备与环境计算机、MATLAB软件环境。三、实验容1.设计一个数字低通滤波器FIR滤波器,其技术指标如下,,分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。〔1矩形窗程序代码如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%录入给定的参数tr_width=wst-wp;N=ceil<1.8*pi/tr_width>+1%根据性能指标确定滤波器长度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%确定理想滤波器的单位脉冲响应hd<n>w_boxcar=boxcar<N>';h=hd.*w_boxcar;%产生矩形窗函数,并对hd<n>加窗得到h<n>subplot<221>;%画出时域和频域图形,以确定设计出来的滤波器是否符合给定的指标stem<n,hd,'filled'>;[Hr,wl]=zerophase<h>;subplot<222>;plot<wl/pi,Hr>;axistight;xlabel<'\omega/\pi'>;ylabel<'H<\omega>'>;subplot<223>;stem<n,h,'filled'>;axistight;xlabel<'n'>;ylabel<'h<n>'>;[H,w]=freqz<h,1>subplot<224>;plot<w/pi,20*log10<abs<H>/max<abs<H>>>>;xlabel<'\omega/\pi'>;ylabel<'dB'>;gridon;%set<gca,'xtick',<0:0.2:1>>;line<[0.3,0.3],[0-100],'linestyle',':'>;line<[0,1],[-0.25-0.25],'linestyle',':'>;line<[0,1],[-50-50],'linestyle',':'>;%axis<[01-20]>;运行结果如下：N=19〔2汉宁窗关键代码如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%录入给定的参数tr_width=wst-wp;N=ceil<6.2*pi/tr_width>+1%根据汉宁窗的性能指标确定滤波器的长度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%确定理想滤波器的单位脉冲响应hd<n>w_hanning=hanning<N>';h=hd.*w_hanning;%产生矩形窗函数,并对hd<n>加窗得到h<n>运行结果如下：N=63〔3海明窗关键代码如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%录入给定的参数tr_width=wst-wp;N=ceil<6.6*pi/tr_width>+1%根据海明窗的性能指标确定滤波器的长度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%确定理想滤波器的单位脉冲响应hd<n>w_hamming=hamming<N>';h=hd.*w_hamming;%产生矩形窗函数,并对hd<n>加窗得到h<n>实验结果如下：N=67〔4布莱克曼窗关键代码如下：wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%录入给定的参数tr_width=wst-wp;N=ceil<11*pi/tr_width>+1%根据布莱克曼的性能指标确定滤波器的长度Nn=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%确定理想滤波器的单位脉冲响应hd<n>w_blackman=blackman<N>';h=hd.*w_blackman;%产生矩形窗函数,并对hd<n>加窗得到h<n>运行结果如下：N=111〔5凯瑟窗wp=0.2*pi;Rp=0.25;wst=0.3*pi;As=50;%录入给定的参数tr_width=wst-wp;N=ceil<<As-7.95>/<2.285*tr_width>>+1%根据凯瑟窗的性能指标确定滤波器的长度Nbeta=0.1102*<As-8.7>;n=0:<N-1>;wc=<wp+wst>/2;alpha=<N-1>/2;hd=<wc/pi>*sinc<<wc/pi>*<n-alpha>>;%确定理想滤波器的单位脉冲响应hd<n>w_kaiser=kaiser<N,beta>';h=hd.*w_kaiser;%产生矩形窗函数,并对hd<n>加窗得到h<n>运行结果如下：N=60上面5种设计方案中,矩形窗和汉宁窗不满足设计要求,而海明窗、布莱克曼窗和凯瑟窗满足设计要求。2.设计一个数字带通FIR滤波器,其技术指标如下下阻带边缘：,下阻带边缘：,上通带边缘：,上阻带边缘：,用布莱克曼窗设计FIR数字滤波器程序代码如下：wp1=0.2*pi;Rp=1;wst1=0.35*pi;As=60;%录入给定的参数wp2=0.65*pi;wst2=0.8*pi;tr_width1=wst1-wp1;tr_width2=wst2-wp2;N1=ceil<11*pi/tr_width1>+1;%根据性能指标确定滤波器长度NN2=ceil<11*pi/tr_width2>+1;N=max<N1,N2>%带通滤波器长度N应为N1、N2中的较大值n=0:<N-1>;wc1=<wp1+wst1>/2;wc2=<wp2+wst2>/2;alpha=<N-1>/2;hd=<wc2/pi>*sinc<<wc2/pi>*<n-alpha>>-<wc1/pi>*sinc<<wc1/pi>*<n-alpha>>w_blackman=blackman<N>';h=hd.*w_blackman;%产生矩形窗函数,并对hd<n>加窗得到h<n>subplot<221>;stem<n,hd,'filled'>;[Hr,wl]=zerophase<h>;subplot<222>;plot<wl/pi,Hr>;axistight;xlabel<'\omega/\pi'>;ylabel<'H<\omega>'>;subplot<223>;stem<n,h,'filled'>;axistight;xlabel<'n'>;ylabel<'h<n>'>;[H,w]=freqz<h,1>;subplot<224>;plot<w/pi,20*log10<abs<H>/max<abs<H>>>>;xlabel<'\omega/\pi'>;ylabel<'dB'>;gridon;set<gca,'xtick',<0:0.2:1>>;line<[0.35,0.35],[0-100],'linestyle',':'>;line<[0.65,0.65],[0-100],'linestyle',':'>;line<[0,1],[-1-1],'linestyle',':'>;line<[0,1],[-60-60],'linestyle',':'>;axis<[0.30.7-20]>;运行结果如下：N=753.采用频率取样设计法设计FIR数字低通滤波器,满足以下指标〔1取N=20,过渡带没有样本。〔2取N=40,过渡带有一个样本,T=0.39。〔3取N=60,过渡带有两个样本,T1=0.5925,T2=0.1099。〔4分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。解：〔1取N=20,过渡带没有样本。关键代码如下：N=20;alpha=<N-1>/2;l=0:N-1;wl=<2*pi/N>*l;Hrs=[1,1,1,zeros<1,15>,1,1];Hdr=[1100];wdl=[00.250.251];%用于绘制理想滤波器幅度函数曲线k1=0:floor<N-1>/2;k2=floor<<N-1>/2>+1:N-1;angH=[<-alpha*<2*pi>/N*k1>,alpha*<2*pi>/N*<N-k2>];H=Hrs.*exp<angH*j>;h=ifft<H,N>;w=[0:500]*pi/500;H=freqz<h,1,w>;[Hrwr]=zerophase<h>;运行结果如下：〔2取N=40,过渡带有一个样本,T=0.39。关键代码如下：N=40;alpha=<N-1>/2;l=0:N-1;wl=<2*pi/N>*l;Hrs=[ones<1,5>,0.39,zeros<1,29>,0.39,ones<1,4>];Hdr=[1100];wdl=[00.250.251];%用于绘制理想滤波器幅度函数曲线k1=0:floor<N-1>/2;k2=floor<<N-1>/2>+1:N-1;angH=[<-alpha*<2*pi>/N*k1>,alpha*<2*pi>/N*<N-k2>];H=Hrs.*exp<angH*j>;h=ifft<H,