数字信号处理实验报告.docx

数字信号处理 信息工程 05111003班数字信号处理实验报告姓名： 学号：专业：信息工程班级：05111003日期：2012年12月目录实验一 利用dft分析信号频谱3实验二 利用fft计算线性卷积20实验三 脉冲响应不变法设计iir数字滤波器31实验四 频率取样法设计fir数字滤波器44实验一 利用dft分析信号频谱（对应实验教材中实验二）一、 实验目的1.加深对dft原理的理解。2.应用dft分析信号频谱。3.深刻理解利用dft分析信号频谱的原理，分析现实过程现象及解决办法。二、 实验设备与环境计算机、matalab软件环境。三、 实验基础理论1.dft和dtft的关系有限长序列的离散时间傅里叶变换在频率区间的n个等分点上的n个取样值可以由下式表示： xej|=2k/n=k=0n-1xne-j2nkn=x(k) 0kn-1 (2-1)由上式可知，序列的n点dft，实际上就是序列的dtft在n个等间隔频率点上样本。2.利用dft求dtft 方法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. 已知，完成如下要求：（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(ejomega); title(dtft); axis tight; subplot(212); plot(w,angle(x1)/pi); % 绘制相频特性曲线xlabel(omega); ylabel(pi) title(phase); axis tight; 实验结果：（2）计算四点dft，把结果显示在（1）所画的图形中。实验代码：n=0:3; x=2,-1,1,1; w=0:0.01*pi:2*pi; x1=x*exp(-j*n*w); % 利用定义计算dtft plot(w, x1); % 绘制dtft xlabel(omega); ylabel(x(k); title(dft of 4 points) axis tight; hold on; 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(dft of 64 points); axis tight; 实验结果：（4）根据实验结果，分析是否可以由dft计算dtft，如果可以，如何实现？答：由定义可知，dft是在数字域以2n为间隔，对dtft的等间隔取样得到的n点序列，因此当序列x(n)的有效值一定时，可以通过补零来增加dft的点数，也就是提高在dtft上的取样频率，从而克服“栅栏效应”。对比(2)和(3)的结果，可以发现(3)的波形已经非常接近(2)了。也就是说对x(n)补零至64点后进行dft，可以很好地克服栅栏效应，模拟x(n)的dtft。因此，我们可以通过对x(n)补零之后做dft来近似计算dtft。用plot函数绘制的补零后的dft波形如下，可以看出几乎与dtft没有区别。2. 考察序列（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(dft of 11 points); axis tight; 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(100 point); axis tight; 实验结果：（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(dft of 100 points); axis tight; 实验结果：分析：这样的频谱符合预期：在=0.48和=0.52附近各有一个冲激。它们刚好对应着两个余弦函数的角频率。（3）根据实验结果，分析怎样提高频谱的分辨率。答：由(1)和(2)中的结果对比可知，当我们增加对x(n)的取样点数（在这里相当于增加了x(n)的记录时间）之后，其频谱分辨率明显提高。但单纯的补零只能避免栅栏效应，而无法提高频谱分辨率。因此，提高频谱分辨率的方法是增长记录时间tp。3. 已知信号，xt=0.15sin2f1t+sin2f2t-0.1sin2f3t其中看，从可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用dft作出频谱分析，确定适合的参数，使得频谱的分辨率符合要求。实验代码：fs=input(sampling freqency= ); % 输入采样频率，保存在fs中tp=input(record time=); % 输入记录时间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); axis tight; 实验结果：首先应当明确的是，什么叫做“频谱分辨率符合要求“。在本题中，一共有三个正弦序列组成了函数x，由信号与系统知识可以知道，正弦函数的频谱结构应当是在虚轴与实轴上各有一个冲激。因此本题中如果能够在数字频率域中，清晰地看到关于=对称的两组谱线，且每组谱线由三个明显的冲激构成，那么就算是”频谱分辨率符合要求“了。本题中，为了分辨出不同频率分量的谱，就要求频率分辨率至少应为不同频率分量的频差的一般，即f=1hz，所以频率分辨率至少为f=0.5hz，因此记录时间tp1f=2。又由nyquist采样定理知，采样频率至少应该是最高频率分量的2倍，因此fs23=6。 基于以上的分析，我们首先取tp=3，fs=10l sampling freqency= 10；record time=3分析：此时已经能够大致看出频谱中有三部分主要频率，但仍不够清晰，下面继续修改参数。l sampling freqency= 10；record time=6分析：当采样频率等于3fh时，从频谱中可以明显地分辨出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；再来观察sampling freqency= 9时的实验图形，发现左半边的三组谱线的中心位置对应的数字频率恰好是0.11、0.22、0.33，说明这三组谱线确实是三个正弦函数所对应的谱线。综上所述，当我们采用大于3fh(即9hz左右)左右的采样频率，记录时间长于6s时，就可以得到符合“频谱分辨率要求“的图像。4. 利用dft近似分析连续时间信号的频谱（幅度谱）。分析采用不同的采样间隔和截取长度计算得结果，并且最终确定合适的参数。 实验代码：tp=input(record time= ); % 输入记录时间fs=input(sampling frequncy= ); % 输入采样频率n=0:1/fs:tp; %以采样频率的倒数作为采样周期，以记录时95f4作为截断时间 x=exp(-0.1*n).*heaviside(n); y=fft(x); subplot(211); % 绘制时域波形stem(n,x); title(time domain) ; axis tight; subplot(212); % 绘制频域波形stem(n/tp*2,abs(y),fill);xlabel(omega(pi);ylabel(x(k);title(dft); axis tight; 实验结果：固定采样频率fs=10hz时，改变记录时间：l record time= 5；sampling frequncy=10l record time= 10；sampling frequncy=10l record time= 100；sampling frequncy=10分析：当我们固定采样频率为10hz时，记录时间从1改变到100。发现当记录时间为100以上时才能得到我们想要的波形。因此选取tp=100s。固定记录时间tp=100s，改变采样频率：l record time= 100；sampling frequncy= 1l record time= 100；sampling frequncy= 5l record time= 100；sampling frequncy= 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=xnhn=m=0n-1x(m)hn-mnrn(n)圆周卷积与线性卷积相等而不产生交叠的必要条件为nl+m-1圆周卷积定理：根据dft性质，x(n)和h(n)的n点圆周卷积的dft等于它们的dft的乘积：dftxnhn=xkh(k)2.快速卷积快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用fft算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：(1)必须选择nl+m-1；为了能使用基-2算法，要求n=2。采用补零的办法使得x(n)和h(n)的长度均为n。(2)计算x(n)和h(n)的n点fft。xn fft x(k)hn fft h(k)(3)组成乘积yk=xkh(k)(3)利用ifft计算y(k)的idft，得到线性卷积y(n)yk i fft y(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=xm ilmil+n-1，i0,0nn-1(3)计算每一段与h(n)的n点圆周卷积，利用fft计算圆周卷积xin fft xikhn fft hk yik fft xikh(k)yik ifft y(n)(4)社区每一段卷积结果的前m-1个样本，连接剩下的样本得到卷积结果y(n)。重叠相加法：设h(n)长度为m，将信号x(n)分解成长为l的子段。以xin表示没断信号，则：xn=i=0xinxin=xn+il, x(n)0nl-10 其他xn*hn=i=0xin*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=0yin-il四、实验内容假设要计算序列x(n)=un-un-l，0nl-1和h(n)=cos0.2n，0nm的线性卷积，完成以下实验内容：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(线性卷积耗时：);ticnewx,newy=meshgrid(x,h); % 仿照上学期信号与系统中求线性卷积方法，先把两序列a=newx.*newy; % 各个元素分别相乘，得到矩阵a。for n=1:lfor m=1:my(n+m-1)=y(n+m-1)+a(m,n);endendtocz=conv(x,h); % 利用conv函数验证上述计算结果是否和正确值吻合。minus=y-z;if(minus1.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(minus1.0e-10) disp(线性卷积=快速卷积);end%-实验结果：l=8m=8线性卷积耗时：elapsed time is 0.000166 seconds.线性卷积算法正确l=8m=8快速卷积耗时：elapsed time is 0.000033 seconds.线性卷积=快速卷积l=16m=16线性卷积耗时：elapsed time is 0.000191 seconds.线性卷积算法正确l=16m=16快速卷积耗时：elapsed time is 0.000038 seconds.线性卷积=快速卷积l=32m=32线性卷积耗时：elapsed time is 0.000315 seconds.线性卷积算法正确l=32m=32快速卷积耗时：elapsed time is 0.000064 seconds.线性卷积=快速卷积l=64m=64线性卷积耗时：elapsed time is 0.000479 seconds.线性卷积算法正确l=64m=64快速卷积耗时：elapsed time is 0.000118 seconds.线性卷积=快速卷积l=256m=256线性卷积耗时：elapsed time is 0.004545 seconds.线性卷积算法正确l=256m=256快速卷积耗时：elapsed time is 0.000095 seconds.线性卷积=快速卷积l=512m=512线性卷积耗时：elapsed time is 0.018785 seconds.线性卷积算法正确l=512m=512快速卷积耗时：elapsed time is 0.000128 seconds.线性卷积=快速卷积l=1024m=1024线性卷积耗时：elapsed time is 0.090633 seconds.线性卷积算法正确l=1024m=1024快速卷积耗时：elapsed time is 0.000204 seconds.线性卷积=快速卷积分析： 根据以上结果可以发现，快速卷积耗时明显小于线性卷积。并且点数越多，这个现象越明显。2.当l=2048，m=256时，比较直接计算线性卷积和快速卷积所需的时间，进一步考察当l=4096，m=256时两种算法所需要的时间。实验代码：代码与(1)相同。实验结果：l=2048m=256线性卷积耗时：elapsed time is 0.046976 seconds.线性卷积算法正确l=2048m=256快速卷积耗时：elapsed time is 0.000402 seconds.线性卷积=快速卷积l=4096m=256线性卷积耗时：elapsed time is 0.103451 seconds.线性卷积算法正确l=4096m=256快速卷积耗时：elapsed time is 0.000717 seconds.线性卷积=快速卷积分析：快速卷积耗费的时间依然小于线性卷积，并且卷积的点数越多，快速卷积的优越性就越明显。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(重叠相加法耗时：);ticfor i=0:num-1 xx=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) yy zeros(1,(l-(i+1)*m); % 将yy补零， y=y+yy; % 再与最终结果y相加。endtoc yyy=conv(x,h); %利用conv函数验证上述计算结果是否和正确值吻合。minus=y-yyy;if(minus1.0e-10) disp(重叠相加法卷积=线性卷积);end实验结果：l=2048m=256重叠相加法耗时：elapsed time is 0.000697 seconds.重叠相加法卷积=线性卷积分析：重叠相加法耗费时间多于快速卷积法(0.000402 seconds)，但远少于线性卷积(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) x zeros(1,n); %创建新序列，这个序列是对x(n)前后都补零得到的。y=zeros(1,l+m-1);p=pow2(nextpow2(n+m-1);disp(重叠保留法耗时：);ticfor i=0:num-1 xx=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) yy zeros(1,l+m-1-(1+i)*(n-m+1); yy=yy(1:l+m-1); y=y+yy;endtoc yyy=conv(x,h); %利用conv函数验证上述计算结果是否和正确值吻合。minus=y-yyy;if(minus1.0e-10) disp(重叠保留法卷积=线性卷积);end实验结果：l=2048m=256重叠保留法耗时：elapsed time is 0.001359 seconds.重叠保留法卷积=线性卷积分析：重叠保留法耗费时间多于快速卷积法(0.000402 seconds)，但少于线性卷积(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);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(h)/max(abs(h);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|(db);subplot(223);plot(w/pi,angle(h)/pi);grid on;xlabel(omega(pi);ylabel(phase of h(ejomega)(pi);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi);ylabel(group delay);实验结果：分析：如上图所示，该三阶巴特沃斯滤波器略微不满足设计指标，即当数字频率=0.5时(对应的模拟角频率为=t=0.54000=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);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(h)/max(abs(h);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|(db);subplot(223);plot(w/pi,angle(h)/pi);grid on;xlabel(omega(pi);ylabel(phase of h(ejomega)(pi);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi);ylabel(group delay);实验结果： 此时，n=6分析：下面我们来检验这个巴特沃斯数字滤波器是否达到了设计指标的要求。根据数字频率和模拟频率的变换公式=t，我们得到数字频率域的指标为：通带截止频率：p=0.2，通带波动：rp=1db阻带截止频率：st=0.3，阻带衰减：as=15db实验结果显示，当=0.2时，幅度的衰减约为衰减到0.9，因此通带内的波动为-10lg0.920.9db15db，满足阻带衰减指标。因此这个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(a2-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);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(h)/max(abs(h);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|(db);subplot(223);plot(w/pi,angle(h)/pi);grid on;xlabel(omega(pi);ylabel(phase of h(ejomega)(pi);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi);ylabel(group delay);实验结果： 此时 n=4分析：下面我们来检验这个切比雪夫i型数字滤波器是否达到了设计指标的要求。实验结果显示，当=0.2时，幅度衰减至0.9，所以通带内的波动为-10lg0.920.9db15db，不仅满足阻带衰减指标，而且还有一定的裕度。因此这个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(a2-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);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(h)/max(abs(h);grid on;xlabel(omega(pi);ylabel(|h(ejomega)|(db);subplot(223);plot(w/pi,angle(h)/pi);grid on;xlabel(omega(pi);ylabel(phase of h(ejomega)(pi);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi);ylabel(group delay);实验结果： 此时 n=4图表 4-1 切比雪夫ii型用脉冲响应不变法数字化的结果分析：下面我们来检验这个切比雪夫ii型数字滤波器是否达到了设计指标的要求。实验结果显示，在00.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/(a2-1)0.5);n=ceil(ellipke(k)*ellipke(1-k12)0.5)/(ell