数字信号处理实验报告

数字信号处理信息工程05111003班PAGE\*MERGEFORMAT41北京理工大学信息与电子学院PAGE\*MERGEFORMAT1数字信号处理实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 姓名：学号：专业：信息工程班级：05111003日期：2012年12月TOC\o"1-3"\h\uHYPERLINK\l"_Toc1943"目录HYPERLINK\l"_Toc11861"实验一利用DFT分析信号频谱3HYPERLINK\l"_Toc10670"实验二利用FFT计算线性卷积20HYPERLINK\l"_Toc32600"实验三脉冲响应不变法设计IIR数字滤波器31HYPERLINK\l"_Toc9188"实验四频率取样法设计FIR数字滤波器44实验一利用DFT分析信号频谱（对应实验教材中实验二）实验目的1.加深对DFT原理的理解。2.应用DFT分析信号频谱。3.深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及解决 办法 鲁班奖评选办法下载鲁班奖评选办法下载鲁班奖评选办法下载企业年金办法下载企业年金办法下载 。实验设备与环境计算机、MATALAB软件环境。实验基础理论1.DFT和DTFT的关系有限长序列的离散时间傅里叶变换在频率区间的N个等分点上的N个取样值可以由下式表示：Xejω|ω=2πk/N=k=0N-1xne-j2πNkn=X(k)0≤k≤N-1(2-1)由上式可知，序列的N点DFT，实际上就是序列的DTFT在N个等间隔频率点上样本。2.利用DFT求DTFT 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 1：由恢复出的方法如图2.1所示：图2.1由图2.1所示流程图可知：由式2-2可以得到其中为内插函数方法2：然而在实际MATLAB计算中，上述插值公式不见得是最好的方法。由于DFT是DTFT的取样值，其相邻的两个频率样本点的间距为,所以如果我们增加数据的长度N，使得得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样可以利用DFT来近似计算DTFT。如果没有更多的数据，可以通过补零来增加数据长度。3.利用DFT分析连续时间信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散化，这里需要进行连个操作：一是采样，二是截断。对于连续非周期信号，按采样间隔T进行采样，截取长度为M，那么对进行N点的频率采样，得到因此，可以将利用DFT分析连续非周期信号频谱的步骤归纳如下：（1）确定时域采样间隔T，得到离散序列；（2）确定截取长度M，得到M点离散序列，这里的为窗函数。（3）确定频域采样点数N，要求。（4）利用FFT计算离散序列的N点DFT，得到。（5）根据式（2-6）由计算采样点的近似值。采用上述方法计算的频谱，需要注意如下三点问题：（1）频谱混叠。如果不满足采样定理的条件，频谱会很出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率的范围。（2）栅栏效应和频谱分辨率。使用DFT计算频谱，得到的结果只是N个频谱样本值，样本值之间的频谱是未知的，就像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成正比，提高频谱分辨率，就要增加记录时间。（3）频谱泄露。对于信号截断会把窗函数的频谱会引入到信号频谱中，造成频谱泄露。解决这问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。因此，要合理选取采样间隔和截取长度，必要时还需考虑适当的窗。对于连续周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上诉方法近似计算。4.可能用到MATLAB函数与代码实验中的DFT运算可以采用MATLAB中提供的FFT来实现。实验 内容 财务内部控制制度的内容财务内部控制制度的内容人员招聘与配置的内容项目成本控制的内容消防安全演练内容 已知，完成如下要求：（1）计算器DTFT，并画出区间的波形。实验代码：n=0:3;x=[2,-1,1,1];w=0:0.01*pi:pi;X1=x*exp(-j*n'*w);%利用定义计算DTFTsubplot(211);plot(w,abs(X1));%绘制幅频特性曲线xlabel('\Omega');ylabel('X(e^j^\Omega)');title('DTFT');axistight;subplot(212);plot(w,angle(X1)/pi);%绘制相频特性曲线xlabel('\omega');ylabel('\pi')title('phase');axistight;实验结果：（2）计算四点DFT，把结果显示在（1）所画的图形中。实验代码：n=0:3;x=[2,-1,1,1];w=0:0.01*pi:2*pi;X1=x*exp(-j*n'*w);%利用定义计算DTFTplot(w,X1);%绘制DTFTxlabel('\Omega');ylabel('X(k)');title('DFTof4points')axistight;holdon;X2=fft(x);n=pi/2.*n;%绘制4点DFTstem(n,X2,'filled');实验结果：（3）对补零，计算其64点DFT，并显示结果。实验代码：n=0:63;x=[2,-1,1,1];z=zeros(1,60);%建立一个零矩阵对x(n)进行补零x=[x,z];%对x(n)进行补零X2=fft(x);%求64点DFTstem(n,X2,'filled');xlabel('k');ylabel('X(k)');title('DFTof64points');axistight;实验结果：（4）根据实验结果，分析是否可以由DFT计算DTFT，如果可以，如何实现？答：由定义可知，DFT是在数字域以2πN为间隔，对DTFT的等间隔取样得到的N点序列，因此当序列x(n)的有效值一定时，可以通过补零来增加DFT的点数，也就是提高在DTFT上的取样频率，从而克服“栅栏效应”。对比(2)和(3)的结果，可以发现(3)的波形已经非常接近(2)了。也就是说对x(n)补零至64点后进行DFT，可以很好地克服栅栏效应，模拟x(n)的DTFT。因此，我们可以通过对x(n)补零之后做DFT来近似计算DTFT。用plot函数绘制的补零后的DFT波形如下，可以看出几乎与DTFT没有区别。考察序列（1）时，用DFT估计的频谱；将补零加长到长度为100点序列用DFT估计的频谱，要求画出相应的波形。实验代码：n=0:10;%取11个样X1=cos(0.48*n*pi)+cos(0.52*n*pi);X2=fft(X1);%做11点DFTsubplot(211);plot(n/5,X2);xlabel('\omega(\pi)');ylabel('X(k)');title('DFTof11points');axistight;subplot(212);a=zeros(1,89);%把11点的采样值序列补零至100点X1=[X1,a];n=0:99;X2=fft(X1);%求补零后的序列的DFTplot(n/50,X2);xlabel('\omega(\pi)');ylabel('X(k)');title('DFT(100point)');axistight;实验结果：（2）时，用DFT估计的频谱，并画出波形实验代码：n=0:99;%题目中要求n上限为100，但是发现n上限为99时图形更准确。X1=cos(0.48*n*pi)+cos(0.52*n*pi);X2=fft(X1);%计算100点DFTplot(n/50,X2);xlabel('\omega(\pi)');ylabel('X(k)');title('DFTof100points');axistight;实验结果：分析：这样的频谱符合预期：在ω=0.48π和ω=0.52π附近各有一个冲激。它们刚好对应着两个余弦函数的角频率。（3）根据实验结果，分析怎样提高频谱的分辨率。答：由(1)和(2)中的结果对比可知，当我们增加对x(n)的取样点数（在这里相当于增加了x(n)的记录时间）之后，其频谱分辨率明显提高。但单纯的补零只能避免栅栏效应，而无法提高频谱分辨率。因此，提高频谱分辨率的方法是增长记录时间tp。已知信号，xt=0.15sin2πf1t+sin2πf2t-0.1sin2πf3t其中看，从可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT作出频谱分析，确定适合的参数，使得频谱的分辨率符合要求。实验代码：fs=input('samplingfreqency=');%输入采样频率，保存在fs中tp=input('recordtime=');%输入记录时间n=0:1/fs:tp;%取采样频率的倒数作为采样间隔，即T=1/fsx=0.15*sin(2*pi*n)+sin(4*pi*n)-0.1*sin(6*pi*n);y=fft(x);plot(n/tp,abs(y));title('DFT');ylabel('X(k)');xlabel('\omega(\pi)');axistight;实验结果：首先应当明确的是，什么叫做“频谱分辨率符合要求“。在本题中，一共有三个正弦序列组成了函数x，由信号与系统知识可以知道，正弦函数的频谱结构应当是在虚轴与实轴上各有一个冲激。因此本题中如果能够在数字频率域中，清晰地看到关于ω=π对称的两组谱线，且每组谱线由三个明显的冲激构成，那么就算是”频谱分辨率符合要求“了。本题中，为了分辨出不同频率分量的谱，就要求频率分辨率至少应为不同频率分量的频差的一般，即∆f=1Hz，所以频率分辨率至少为F=0.5Hz，因此记录时间tp≥1F=2。又由Nyquist采样定理知，采样频率至少应该是最高频率分量的2倍，因此fs≥2×3=6。基于以上的分析，我们首先取tp=3，fs=10samplingfreqency=10；recordtime=3分析：此时已经能够大致看出频谱中有三部分主要频率，但仍不够清晰，下面继续修改参数。samplingfreqency=10；recordtime=6分析：当采样频率等于3∙fh时，从频谱中可以明显地分辨出f1、f2、f3所对应的频谱。下面我们用模拟频率与数字频率转换的公式ω=Ω∙T，来验证这些谱线是否属于三个正弦函数。三个正弦函数的模拟频率分别为:Ω1=1*2π，Ω2=2*2π，Ω3=3*2π因此其数字频率分别为:ω1=1*2π*118=0.11π；ω2=2*2π*118=0.22π；ω3=3*2π*118=0.33π；再来观察samplingfreqency=9时的实验图形，发现左半边的三组谱线的中心位置对应的数字频率ω恰好是0.11π、0.22π、0.33π，说明这三组谱线确实是三个正弦函数所对应的谱线。综上所述，当我们采用大于3∙fh(即9Hz左右)左右的采样频率，记录时间长于6s时，就可以得到符合“频谱分辨率要求“的图像。利用DFT近似分析连续时间信号的频谱（幅度谱）。分析采用不同的采样间隔和截取长度计算得结果，并且最终确定合适的参数。实验代码：tp=input('recordtime=');%输入记录时间fs=input('samplingfrequncy=');%输入采样频率n=0:1/fs:tp;%以采样频率的倒数作为采样周期，以记录时95F4作为截断时间x=exp(-0.1*n).*heaviside(n);y=fft(x);subplot(211);%绘制时域波形stem(n,x);title('TimeDomain');axistight;subplot(212);%绘制频域波形stem(n/tp*2,abs(y),'fill');xlabel('\omega(\pi)');ylabel('X(k)');title('DFT');axistight;实验结果：固定采样频率fs=10Hz时，改变记录时间：recordtime=5；samplingfrequncy=10recordtime=10；samplingfrequncy=10recordtime=100；samplingfrequncy=10分析：当我们固定采样频率为10Hz时，记录时间从1改变到100。发现当记录时间为100以上时才能得到我们想要的波形。因此选取tp=100s。固定记录时间tp=100s，改变采样频率：recordtime=100；samplingfrequncy=1recordtime=100；samplingfrequncy=5recordtime=100；samplingfrequncy=10分析：由上述三种采样频率下得到的图像可知，取fs=10Hz就可得到比较准确合理的波形。可以发现，当记录时间不足时，即使采样频率足够，也无法得到准确的图形，因此足够的；当记录时长足够而采样频率不足时同样无法得到准确波形。这是因为记录时间决定了图形的频谱分辨率，而采样频率则决定着波形是否发生混叠。综上所述，选取采样频率fs=10Hz，记录时间tp=100s。实验心得与总结通过这次实验，我加深了对于DFT的理解，并且学会了用DFT分析信号的频谱。同时，我更加深刻地理解了记录时间tp和采样频率fs与频谱分辨率的关系，掌握了提高频谱分辨率的方法。另外，这次实验中第一次尝试用MATLAB的script脚本模式调试程序，体会了其方便之处。实验二利用FFT计算线性卷积（对应实验教材中实验三）实验目的1.掌握利用FFT计算线性卷积的原理及具体实现方法。2.加深理解重叠相加法和重叠保留法。3.考察利用FFT计算线性卷积各种方法的适用范围。实验设备与环境计算机、MATLAB软件环境。实验基础理论1.线性卷积与圆周卷积设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为yln=xn*hn=m=-∞∞x(m)hn-myln的长度为L+M-1。x(n)和h(n)的圆周卷积为yn=xn⊙hn=m=0N-1x(m)hn-mNRN(n)圆周卷积与线性卷积相等而不产生交叠的必要条件为N≥L+M-1圆周卷积定理：根据DFT性质，x(n)和h(n)的N点圆周卷积的DFT等于它们的DFT的乘积：DFTxn⊙hn=XkH(k)2.快速卷积快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：(1)必须选择N≥L+M-1；为了能使用基-2算法，要求N=2γ。采用补零的办法使得x(n)和h(n)的长度均为N。(2)计算x(n)和h(n)的N点FFT。xnFFTX(k)hnFFTH(k)(3)组成乘积Yk=XkH(k)(3)利用IFFT计算Y(k)的IDFT，得到线性卷积y(n)YkIFFTy(n)3.分段卷积我们考察单位取样响应为h(n)的线性系统，输入为x(n)，输出为y(n)，则yn=xn*hn如果x(n)极长时，如果要等x(n)全部集齐时再开始进行卷积，会使输出有较大延时；如果序列太长，需要大量存储单元。为此，我们把x(n)分段，为别求出每段的卷积，合在一起得到最后的总输出。这称为分段卷积。分段卷积可以细分为重叠保留法和重叠相加法。重叠保留法：设x(n)的长度为Nx，h(n)的长度为M。把序列x(n)分成多段N点序列xi(n)，每段雨前一段重写M-1个样本。并在第一个输入段前面补M-1个零。计算每一段与h(n)的圆周卷积，其结果中前M-1个不等与线性卷积，应当舍去，只保留后面N-M+1个正确的输出样本，把它们合起来得到总的输出。利用FFT实现重叠保留法的步骤如下：(1)在x(n)前面填充M-1个零，扩大以后的序列为xn={0,0,⋯0,x(n)}(2)将x(n)分为若干段N点子段，设L=N-M+1为每一段的有效长度，则第i段的数据为：xin=xmiL≤m≤iL+N-1，i≥0,0≤n≤N-1(3)计算每一段与h(n)的N点圆周卷积，利用FFT计算圆周卷积xinFFTXikhnFFTHkYikFFTXikH(k)YikIFFTy(n)(4)社区每一段卷积结果的前M-1个样本，连接剩下的样本得到卷积结果y(n)。重叠相加法：设h(n)长度为M，将信号x(n)分解成长为L的子段。以xin表示没断信号，则：xn=i=0∞xinxin=xn+iL,x(n)0≤n≤L-10其他xn*hn=i=0∞xin*hn每一段卷积yin的长度为L+M-1，所以在做求和时，相邻两段序列由M-1个样本重叠，即前一段的最后M-1个样本和下一段前M-1个样本序列重叠，这个重叠部分相加，再与不重叠的部分共同组成y(n)。利用FFT实现重叠保留法的步骤如下：(1)将x(n)分为若干L点子段xin。(2)计算每一段与h(n)的卷积，根据快速卷积法利用FFT计算卷积。(3)将各段相加，得到输出y(n)。y(n)=i=0∞yin-iL四、实验内容假设要计算序列x(n)=un-un-L，0≤n≤L-1和h(n)=cos0.2πn，0≤n≤M的线性卷积，完成以下实验内容：1.设L=M，根据线性卷积的表达式和快速卷积的原理，分别编程实现计算两个序列线性卷积的方法，比较序列长度分别为8,16,32,64,256,512,1024时，两种方法计算线性卷积所需的时间。实验代码：%----------------用自己编制的程序，根据原理计算线性卷积-----------L=input('L=');M=input('M=');x=ones(1,L);k=0:1:M-1;h=cos(0.2*pi*k);y=zeros(1,M+L-1);%先用0对y占位，以避免在循环中反复添加新元素造成的时间浪费。disp('线性卷积耗时：');tic[newx,newy]=meshgrid(x,h);%仿照上学期信号与系统中求线性卷积方法，先把两序列A=newx.*newy;%各个元素分别相乘，得到矩阵A。forn=1:Lform=1:My(n+m-1)=y(n+m-1)+A(m,n);endendtocz=conv(x,h);%利用conv函数验证上述计算结果是否和正确值吻合。minus=y-z;if(minus<1.0e-10)%如果两者吻合，显示“线性卷积算法正确”。disp('线性卷积算法正确');end%------------------------计算快速卷积-------------------------L=input('L=');M=input('M=');x=ones(1,L);k=0:1:M-1;h=cos(0.2*pi*k);N=pow2(nextpow2(L+M-1));%利用pow2函数，确定FFT的点数。disp('快速卷积耗时：');ticX=fft(x,N);H=fft(h,N);Y=X.*H;y=ifft(Y,N);y=y(1:L+M-1);tocyy=conv(x,h);%利用conv函数验证上述计算结果是否和正确值吻合。minus=y-yy;if(minus<1.0e-10)disp('线性卷积=快速卷积');end%--------------------------------------------------------------------实验结果：L=8M=8线性卷积耗时：Elapsedtimeis0.000166seconds.线性卷积算法正确L=8M=8快速卷积耗时：Elapsedtimeis0.000033seconds.线性卷积=快速卷积L=16M=16线性卷积耗时：Elapsedtimeis0.000191seconds.线性卷积算法正确L=16M=16快速卷积耗时：Elapsedtimeis0.000038seconds.线性卷积=快速卷积L=32M=32线性卷积耗时：Elapsedtimeis0.000315seconds.线性卷积算法正确L=32M=32快速卷积耗时：Elapsedtimeis0.000064seconds.线性卷积=快速卷积L=64M=64线性卷积耗时：Elapsedtimeis0.000479seconds.线性卷积算法正确L=64M=64快速卷积耗时：Elapsedtimeis0.000118seconds.线性卷积=快速卷积L=256M=256线性卷积耗时：Elapsedtimeis0.004545seconds.线性卷积算法正确L=256M=256快速卷积耗时：Elapsedtimeis0.000095seconds.线性卷积=快速卷积L=512M=512线性卷积耗时：Elapsedtimeis0.018785seconds.线性卷积算法正确L=512M=512快速卷积耗时：Elapsedtimeis0.000128seconds.线性卷积=快速卷积L=1024M=1024线性卷积耗时：Elapsedtimeis0.090633seconds.线性卷积算法正确L=1024M=1024快速卷积耗时：Elapsedtimeis0.000204seconds.线性卷积=快速卷积分析：根据以上结果可以发现，快速卷积耗时明显小于线性卷积。并且点数越多，这个现象越明显。2.当L=2048，M=256时，比较直接计算线性卷积和快速卷积所需的时间，进一步考察当L=4096，M=256时两种算法所需要的时间。实验代码：代码与(1)相同。实验结果：L=2048M=256线性卷积耗时：Elapsedtimeis0.046976seconds.线性卷积算法正确L=2048M=256快速卷积耗时：Elapsedtimeis0.000402seconds.线性卷积=快速卷积L=4096M=256线性卷积耗时：Elapsedtimeis0.103451seconds.线性卷积算法正确L=4096M=256快速卷积耗时：Elapsedtimeis0.000717seconds.线性卷积=快速卷积分析：快速卷积耗费的时间依然小于线性卷积，并且卷积的点数越多，快速卷积的优越性就越明显。3.编程实现利用重叠相加法计算两个序列的线性卷积，考察L=2048且M=256使计算线性卷积所需的时间，与第二题的结果进行比较。实验代码：L=input('L=');M=input('M=');num=L/M;k=0:1:M-1;h=cos(0.2*pi*k);x=ones(1,L);y=zeros(1,L+M-1);%先用0对y占位，避免在循环中反复添加新元素造成的时间浪费。N=pow2(nextpow2(M+M-1));%用pow2函数，确定FFT的点数。disp('重叠相加法耗时：');ticfori=0:num-1xx=x((i*M+1):(i+1)*M);%对x(n)分段取值。X=fft(xx,N);H=fft(h,N);Y=X.*H;yy=ifft(Y,N);yy=yy(1:(2*M-1));yy=[zeros(1,i*M)yyzeros(1,(L-(i+1)*M))];%将yy补零，y=y+yy;%再与最终结果y相加。endtocyyy=conv(x,h);%利用conv函数验证上述计算结果是否和正确值吻合。minus=y-yyy;if(minus<1.0e-10)disp('重叠相加法卷积=线性卷积');end实验结果：L=2048M=256重叠相加法耗时：Elapsedtimeis0.000697seconds.重叠相加法卷积=线性卷积分析：重叠相加法耗费时间多于快速卷积法(0.000402seconds)，但远少于线性卷积(0.046976)。这说明重叠保留法的计算速度强于线性卷积，但略逊于快速卷积。4.编程实现利用重叠保留法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积所需的时间，与第二题的结果进行比较。实验代码：L=input('L=');M=input('M=');N=2*M;num=ceil((L+M-1)/(N-M+1));%利用ceil函数取整，得到的是分段的段数。k=0:1:M-1;h=cos(0.2*pi*k);x=ones(1,L);x_1=[zeros(1,M-1)xzeros(1,N)];%创建新序列，这个序列是对x(n)前后都补零得到的。y=zeros(1,L+M-1);P=pow2(nextpow2(N+M-1));disp('重叠保留法耗时：');ticfori=0:num-1xx=x_1((i*(N-M+1)+1):i*(N-M+1)+N);%对补零后的序列分段取值。X=fft(xx,P);H=fft(h,P);Y=X.*H;yy=ifft(Y,P);yy=yy(M:N);%yy=[zeros(1,i*(N-M+1))yyzeros(1,L+M-1-(1+i)*(N-M+1))];yy=yy(1:L+M-1);y=y+yy;endtocyyy=conv(x,h);%利用conv函数验证上述计算结果是否和正确值吻合。minus=y-yyy;if(minus<1.0e-10)disp('重叠保留法卷积=线性卷积');end实验结果：L=2048M=256重叠保留法耗时：Elapsedtimeis0.001359seconds.重叠保留法卷积=线性卷积分析：重叠保留法耗费时间多于快速卷积法(0.000402seconds)，但少于线性卷积(0.046976)。这说明重叠保留法的计算速度强于线性卷积，但逊于快速卷积。五、实验心得通过这次实验，我掌握了利用FFT算法计算线性卷积的原理和方法。同时，通过对比利用原理直接计算线性卷积和利用快速卷积算法计算线性卷积，我更加深刻地体会到了在计算多点数的线性卷积时FFT算法的优越性。另外，通过重叠保留法和重叠相加法与上述两种算法的对比，我了解到了它们彼此之间运算速度的快慢关系，从而加深了对重叠相加法和重叠保留法的理解。实验三脉冲响应不变法设计IIR数字滤波器（对应实验教材中实验五）一、实验目的1.掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。二、实验设备与环境计算机、MATLAB软件环境。三、实验基础理论1.基本原理从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应，h(n)等于的取样值。2.变换方法思路：(1)将进行部分分式展开(2)对进行拉式变换(3)对时域采样得到h(n)(4)对h(n)进行z变换3．设计步骤(1)确定数字滤波器的性能指标。(2)将数字滤波器频率指标转换成响应的模拟滤波器频率指标(3)根据指标，，和设计模拟滤波器。(4)将展成部分分式形式。(5)把模拟极点转换成数字极点，得到数字滤波器。可见至H(z)间的变换关系为在matlab中提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换。[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系数向量bz和az，fs为采样频率（默认为1）。四、实验内容1.设采样频率为fs=4kHz，采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=1kHz。实验代码：N=3;OmegaC=1000*2*pi;%利用线性变换公式Ω=ω/T=ω∙fs求出模拟滤波器的指标fs=4000;[b,a]=butter(N,OmegaC,'s');%求模拟巴特沃斯滤波器[bz,az]=impinvar(b,a,fs);%数字化[H,w]=freqz(bz,az);%求幅频响应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(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('GroupDelay');实验结果：分析：如上图所示，该三阶巴特沃斯滤波器略微不满足设计指标，即当数字频率ω=0.5π时(对应的模拟角频率为Ω=ωT=0.5π×4000=2000πrad/s=1kHz)，幅度衰减为2.5dB左右，而没有达到3dB。要想使滤波器满足设计指标，可以提高阶数N。2.设采样频率为fs=10kHz，涉及数字低通滤波器满足如下指标：通带截止频率：fp=1kHz，在此处键入公式。通带波动：Rp=1dB阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型、和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合试验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定的指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。(1)巴特沃斯低通滤波器实验代码：wp=1000*2*pi;%通带截止频率ws=1500*2*pi;%阻带截止频率Rp=1;%通带波动As=15;%阻带波动N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)));%阶数OmegaC=wp/((10^(Rp/10)-1)^(1/(2*N)));%截止频率[b,a]=butter(N,OmegaC,'s');%求巴特沃斯滤波器fs=10000;%采样频率[bz,az]=impinvar(b,a,fs);%将模拟滤波器数字化[H,w]=freqz(bz,az);%求幅频响应ε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(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('GroupDelay');实验结果：此时，N=6分析：下面我们来检验这个巴特沃斯数字滤波器是否达到了设计指标的要求。根据数字频率和模拟频率的变换公式ω=Ω∙T，我们得到数字频率域的指标为：通带截止频率：ωp=0.2π，通带波动：Rp=1dB阻带截止频率：ωst=0.3π，阻带衰减：As=15dB实验结果显示，当ω=0.2π时，幅度的衰减约为衰减到0.9，因此通带内的波动为-10lg0.92≈0.9dB<1dB，满足通带指标；当ω=0.3π时，幅度大约衰减到0.17，因此阻带衰减为As=-10lg0.172≈15.4dB>15dB，满足阻带衰减指标。因此这个6阶巴特沃斯滤波器完全满足指标要求。(2)切比雪夫I型低通滤波器实验代码：Rp=1;%通带波动As=15;%阻带波动wp=1000*2*pi;%通带截止频率ws=1500*2*pi;%阻带截止频率E=((10^(Rp/10))-1)^0.5;%通带波动系数A=10^(As/20);%阻带衰减参数N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)));%阶数[b,a]=cheby1(N,Rp,wp,'s');%求切比雪夫I型滤波器fs=10000;[bz,az]=impinvar(b,a,fs);%将模拟滤波器数字化[H,w]=freqz(bz,az);%求频率响应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(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('GroupDelay');实验结果：此时N=4分析：下面我们来检验这个切比雪夫I型数字滤波器是否达到了设计指标的要求。实验结果显示，当ω=0.2π时，幅度衰减至0.9，所以通带内的波动为-10lg0.92≈0.9dB<1dB，满足通带指标；当ω=0.3π时，幅度大约衰减到略小于0.1，因此阻带衰减为As=-10lg0.12≈20dB>15dB，不仅满足阻带衰减指标，而且还有一定的裕度。因此这个4阶切比雪夫I型滤波器完全满足指标要求。(3)切比雪夫II型低通滤波器实验代码：Rp=1;%通带波动As=15;%阻带波动wp=1000*2*pi;%通带截止频率ws=1500*2*pi;%阻带截止频率E=((10^(Rp/10))-1)^0.5;%通带波动系数A=10^(As/20);%阻带衰减参数N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)));%阶数fs=10000;[b,a]=cheby2(N,As,ws,'s');%求切比雪夫II型滤波器[bz,az]=impinvar(b,a,fs);%将模拟滤波器数字化[H,w]=freqz(bz,az);%求频率响应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(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('GroupDelay');实验结果：此时N=4图表4-1切比雪夫II型用脉冲响应不变法数字化的结果分析：下面我们来检验这个切比雪夫II型数字滤波器是否达到了设计指标的要求。实验结果显示，在0<ω<0.2π的通带范围内，带内波动8dB，这显然不满足设计要求；当ω=0.3π时，阻带衰减约为6dB，这也不满足要求。另外我们可以明显看到切比雪夫II型滤波器在通带内几乎无波动，但在阻带内有明显的波动，充分利用了阻带的裕度。因此这个4阶切比雪夫II型滤波器不满足设计指标为什么会出现这样的情况？原因在于用脉冲响应不变法的最大缺点就是在对频谱的周期延拓过程中，会发生混叠。由于切比雪夫II型滤波器在阻带内有很大的波动，所以就造成了上图中所示的混叠的结果。如果在此处使用双线性变换法，就会避免混叠的问题，得到正确结果。下面我们通过实验验证：只须将数字化的代码修改为：[bz,az]=bilinear(b,a,fs);用双线性变换法改进后的结果：图表4-2切比雪夫II型用双线性变换法数字化的结果可以看出，对于同一个切比雪夫II型模拟滤波器，用双线性变换法数字化后的结果符合各项设计指标，并且波形也完全符合切比雪夫II型的特点（通带无波动，阻带有波动）。(4)椭圆模拟原型滤波器实验代码：Rp=1;%通带波动As=15;%阻带波动wp=1000*2*pi;%通带截止频率ws=1500*2*pi;%阻带截止频率E=((10^(Rp/10))-1)^0.5;%通带波动系数A=10^(As/20);%阻带衰减参数k=wp/ws;%增益k1=E/((A^2-1)^0.5);N=ceil(ellipke(k)*ellipke((1-k1^2)^0.5)/(ellipke(k1)*ellipke((1-k^2)^0.5)));%阶数fs=10000;[b,a]=ellip(N,Rp,As,wp,'s');[bz,az]=impinvar(b,a,fs);%将模拟滤波器数字化[H,w]=freqz(bz,az);%求频率响应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(bz,az,w);subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('GroupDelay');实验结果：此时，N=3；图表4-3椭圆滤波器用脉冲响应变换法数字化的结果分析：下面我们来检验这个椭圆数字滤波器是否达到了设计指标的要求。实验结果显示，在0<ω<0.2π的通带范围内，带内波动约为4dB，这明显不满足设计要求的小于1dB；当ω=0.3π时，阻带衰减13dB，这也不满足要求的15dB。因此这个3阶椭圆滤波器不符合指标要求。这样的结果依然归咎于脉冲响应不变法会引起频谱的混叠，而恰好椭圆滤波器在阻带内有波动。下面再次用双线性变换法改进。只需将数字化的代码修改为：[bz,az]=bilinear(b,a,fs);%将模拟滤波器数字化用双线性变化法改进后的结果：改进后的滤波器明显符合各项指标要求。另外，我们可以看通带和阻带都由明显的带内波动，说明它充分利用的通带和阻带的带内裕度，这也是相同指标下椭圆滤波器阶数最低的原因。总结以上各项中的分析：首先，我们发现对于相同的指标，三种类型的滤波器所需的阶数不同。具体而言，巴特沃斯滤波器所需阶数最高(本实验中N=6)，切比雪夫次之(本实验中N=4)，椭圆滤波器最低(本实验中N=3)。造成这种现象的原因在于巴特沃斯滤波器严格单减，因此不论在通带还是阻带都造成了浪费；切比雪夫滤波器在通带或阻带允许有波动(ripple)；而椭圆滤波器在通带和阻带都允许有波动。其次，考察各类滤波器是否实现指标时我们发现，只有巴特沃斯和切比雪夫I型达到了技术指标。Rp通带波动As阻带衰减是否达标巴特沃斯0.9dB15.4dB是切比雪夫I0.9dB20dB是切比雪夫II8dB6dB否椭圆4dB13dB否切比雪夫II型和椭圆滤波器不满足指标的原因是其阻带都有波动，因此在用脉冲响应不变法数字化的过程中，频谱发生了混叠，造成了严重的失真。改进的方法就是使用双线性变换法，对此我也在前文中进行了验证。最后，来分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。脉冲响应不变法最大的优点是频率坐标的变换是线性的ω=Ω∙T。因此，如果模拟滤波器的频响是限带于折叠频率以内，那么通过变换后得到的数字滤波器的频响可以不失真地反映原响应与频率的关系。它最大的缺点是频谱的周期延拓效应，因此，只能用于限带的品相特性，如衰减特性很好的低通，而且高频衰减越大，频响的混叠效应就越小。对于高通、带阻滤波器，由于他们在高频部分不衰减，当一定要追求频率线性关系而采用脉冲响应不变法时，必须先对模拟高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频率，然后再转变为数字滤波器以避免混叠失真。否则就会出现本次实验中切比雪夫II型和椭圆滤波器频谱严重失真的结果。五、实验心得1.通过这次实验，我主要练习了用脉冲响应不变法设计数字滤波器的方法，熟悉了整个设计和实现的流程。2.在分析实验结果的过程中，我曾经一度怀疑切比雪夫II型滤波器和椭圆滤波器的设计有问题，因为它们都没有达到设计指标。但经过仔细考虑我发现是脉冲响应不变法的频谱周期延拓效应造成了这样的结果，这让我对于脉冲响应不变法的缺点有了更加深刻的印象。3.另外，我还练习了用双线性变换法来改进切比雪夫II型和椭圆数字滤波器，得到了不错的结果。实验四频率取样法设计FIR数字滤波器（对应实验教材中实验八）一、实验目的掌握频率取样法设计FIR数字滤波器的原理及具体方法二、实验设备与环境计算机、MATLAB软件环境三、实验基础理论1.基本原理频率取样法从频域出发，把理想的滤波器等间隔采样得到,将作为实际设计滤波器的：得到以后可以由来确定唯一确定滤波器的单位脉冲响应,可以由求得：其中为内插函数：有求得的频率响应将逼近。如果我们设计的是线性相位FIR滤波器，则的幅度和相位满足线性相位滤波器的约束条件。我们将表示为如下形式当为实数，则由此得到即为中心偶对称。在利用线性相位条件可知，对于1型和2型线性相位滤波器：对于3型和4型线性相位滤波器2.设计步骤（1）由给定的理想滤波器给出和。（2）由求得（3）根据求得或四、实验内容及实验结果1.采用频率取样法设计FIR数字低通滤波器，满足以下指标(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=[1,1,0,0];%理想滤波器幅频响应的分段函数值wdL=[0,0.25,0.25,1];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(j*angH);%理想滤波器频率响应h=ifft(H,N);%FIR滤波器的单位脉冲响应w=[0:500]*pi/500;[Hr,wr]=zerophase(h);%计算FIR滤波器的幅度函数subplot(221);%绘制理想滤波器幅频响应及抽样值plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);%绘制单位脉冲响应stem(L,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);%绘制FIR滤波器幅频响应plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1),'o');axis([0,1,-0.2,1.2]);gridon;xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);%绘制FIR滤波器幅频响应的dB图plot(wr/pi,20*log10((abs(Hr)/max(abs(Hr)))));axis([0,1,-50,5]);gridon;xlabel('\omega(\pi)');ylabel('dB')实验结果：分析：由图可见，N=20时，滤波器对通带波动和阻带衰减两个指标都不满足(2)取N=40，过渡带有一个样本，T=0.39。实验代码：只须对(1)的代码相应位置做下列改动：N=40;Hrs=[1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1];Hdr=[1,1,0.39,0,0];wdL=[0,0.2,0.25,0.3,1];实验结果：分析：由图可见，N=40时，滤波器满足阻带衰减要求，但不满足通带波动要求。(3)取N=60，过渡带有两个样本，T1=0.5925，T2=0.1099。实验代码：只须对(1)的代码的相应位置做下列改动：N=60;Hrs=[1,1,1,1,1,1,1,0.5925,0.1099,zeros(1,43),0.1099,0.5925,1,1,1,1,1,1];Hdr=[1,1,0.5925,0.1099,0,0];wdL=[0,0.2,7/30,8/30,0.3,1];实验结果：分析：由图可见，N=60时，设计的滤波器完全满足要求。(4)分别讨论采用上述方法设计的数字滤波器是否都能满足给定的指标要求。将各阶数的滤波器是否满足各项指标的情况展示在下表中：通带波动0.25dB阻带衰减50dB是否满足指标N=20不满足不满足否N=40不满足满足否N=60满足满足是由此可见，滤波器的阶数越高，其性能也越好。2.采用频率取样技术设计下面的高通滤波器对于高通滤波器，N必须为奇数（或1型滤波器）。选择N=33，过渡带有两个样本，过渡带的最优值为T1=0.1095，T2=0.598.实验代码：N=33;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L;Hrs=[0,0,0,0,0,0,0,0,0,0,0,0.1095,0.598,1,1,1,1,1,1,1,0.598,0.1095,0,0,0,0,0,0,0,0,0,0,0];%理想滤波器幅频响应取样值Hdr=[0,0,0.1095,0.598,1,1];wdL=[0,0.6,2/3,24/33,0.8,1];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(j*angH);h=ifft(H,N);w=[0:500]*pi/500;Hr=h*cos((alpha-L)'*w);%计算幅度函数subplot(221);%绘制理想滤波器幅频响应及抽样值plot(wdL,Hdr,wL(1:17)/pi,Hrs(1:17),'o');axis([0,1,-0.1,1.1]);xlabel('\omega(\pi)');ylabel('Hr(k)');subplot(222);%绘制单位脉冲响应stem(L,h,'filled');axis([0,N-1,-0.1,0.3]);xlabel('n');ylabel('h(n)');subplot(223);%绘制幅频响应plot(w/pi,Hr,wL(1:17)/pi,Hrs(1:17),'o');axis([0,1,-0.2,1.2]);;gridon;xlabel('\omega(\pi)');ylabel('Hr(\omega)');subplot(224);%绘制幅频响应的dB图plot(w/pi,20*log10((abs(Hr)/max(abs(Hr)))));axis([0,1,-50,5]);gridon;xlabel('\omega(\pi)');ylabel('dB')