数字信号处理实验报告

数字信号处理实验报告目录实验1利用DFT分析信号频谱一、实验目的二、实验原理三、实验内容四、实验总结实验2利用FFT计算线性卷积一、实验目的二、实验原理三、实验内容四、实验总结实验3IIR数字滤波器设计一、实验目的二、实验原理三、实验内容四、实验总结实验4FIR数字滤波器设计一、实验目的二、实验原理三、实验内容四、实验总结实验1利用DFT分析信号频谱一、实验目的1.加深对DFT原理的理解。2.应用DFT分析信号频谱。3.深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及解决办法。二、实验原理1、DFT和DTFT的关系有限长序列的离散时间傅里叶变换在频率区间的N个等分点上的N个取样值可以由下式表示：由上式可知，序列的N点DFT，实际上就是序列的DTFT在N个等间隔频率点上样本。2、利用DFT求DTFT方法1：由恢复出的方法如图2.1所示：图2.1.由N点DFT恢复频谱DTFT的流程由图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来实现。DTFT可以利用MATLAB矩阵运算的方法进行计算。三、实验内容1.，完成如下要求：（1）计算其DTFT，并画出区间的波形。（2）计算4点DFT，并把结果显示在（1）所画的图形中。（3）对补零，计算64点DFT，并显示结果。（4）是否可以由DFT计算DTFT，如果可以，请编程实现。代码及图形：（1）计算其DTFT，并画出区间的波形。n=0:3;x=[2-111];w=-pi:0.001:pi;X=x*exp(-j*n'*w);subplot(211);plot(w,abs(X));xlabel('w');title('Magnitude');%画出幅度谱>>axistightsubplot(212);plot(w,angle(X)/pi);>>xlabel('w');title('Phase');%画出相位谱axistight（2）计算4点DFT，并把结果显示在（1）所画的图形中。%只需在(1)的基础上加上部分代码，并稍作更改n=0:3;x=[2-111];w=-pi:0.001:2*pi;X=x*exp(-j*n'*w);Y=fft(x);>>subplot(211);plot(w,abs(X));>>holdon>>stem([0:pi/2:3*pi/2],abs(Y));%画出4点DFT>>xlabel('w');title('Magnitude');>>axistight>>subplot(212);plot(w,angle(X)/pi);holdon>>stem([0:pi/2:3*pi/2],angle(Y));%画出4点DFT>>xlabel('w');title('Phase');axistight（3）对补零，计算64点DFT，并显示结果。n=0:3;x=[2-111];Y=fft(x,64);%为了计算64点DFTk=0:63;subplot(211);stem(k,abs(Y));xlabel('k');title('|X(k)|');axistightsubplot(212);stem(k,angle(Y));xlabel('k');title('φ(X(k))');axistight是否可以由DFT计算DTFT，如果可以，请编程实现。由（1）和（2）（3）可知，DFT相当于在DTFT的频谱上取样。当DFT的点数比较少时，无法根据根据DFT的频谱得到理想的DTFT频谱；如果通过采取在序列后面补零的方法，则可以提高采样密度，使DFT的点数更多，从而得到精细的DTFT频谱，且补零越多，频谱越精细。2.考察序列（1）时，用DFT估计的频谱；将补零加长到长度为100点序列用DFT估计的频谱。要求画出相应波形。（2）时，用DFT估计的频谱，并画出波形。（3）根据实验结果，分析怎样提高频谱分辨率。代码及图形：（1）时，用DFT估计的频谱；将补零加长到长度为100点序列用DFT估计的频谱。要求画出相应波形。步骤一：时，用DFT估计的频谱n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);>>Y=fft(x);>>subplot(211);plot(abs(Y));>>xlabel('0<=n<=10,DFT估计x(n)频谱');title('Magnitude');>>axistight>>subplot(212);plot(angle(Y));>>xlabel('0<=n<=10,DFT估计x(n)频谱');title('Phase');axistight步骤二：将补零加长到长度为100点序列用DFT估计的频谱n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);>>Y=fft(x,100);>>subplot(211);plot(abs(Y));>>xlabel('0<=n<=10,补零至100点，DFT估计x(n)频谱');title('Magnitude');>>axistight>>subplot(212);plot(angle(Y));>>xlabel('0<=n<=10,补零至100点，DFT估计x(n)频谱');title('Phase');axistight（2）时，用DFT估计的频谱，并画出波形。n=0:100;x=cos(0.48*pi*n)+cos(0.52*pi*n);>>Y=fft(x);>>subplot(211);plot(abs(Y));>>xlabel('0<=n<=100，DFT估计x(n)频谱');title('Magnitude');>>axistight>>subplot(212);plot(angle(Y));>>xlabel('0<=n<=100，DFT估计x(n)频谱');title('Phase');axistight根据实验结果，分析怎样提高频谱分辨率。由以上实验结果可知可以通过增大截取长度和增加补零的个数来使得序列列长变长，从而提高频谱分辨率。已知信号，其中，，。从的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。代码及图形：t=0:0.1:15;x=0.15*sin(2*pi*t)+sin(2*pi*2*t)-0.1*sin(2*pi*3*t);>>Y=fft(x);>>subplot(211);stem(abs(Y));>>xlabel('DFT估计x(n)频谱');title('Magnitude');>>axistight>>subplot(212);stem(angle(Y));>>xlabel('DFT估计x(n)频谱');title('Phase');axistight利用DFT近似分析连续时间信号的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。代码及图形：n=0:2:120;x=exp(-0.1*n);X=fft(x);k=0:60;stem(k,abs(X));xlabel('k');title('Magnitude');%采样区间可以为0到120，采样间隔为2。实验总结本实验研究怎样用DFT分析信号频谱的问题。在理论学习中，已经了解到DFT由DFS取主值序列得到，而DFS可以从DTFT推导而来，这说明DFT与DTFT之间有着密切关系。实际上，DTFT的频谱抽样后能够得到DFT的频谱；即DTF的点数越多，越能通过拟合得到精细的DTFT频谱。在实验过程中，我对提高频谱分辨率的方法有了更加深刻的认识，通过序列补零和选取合适的截取长度都可以达到目的。当然在实验中也遇到不少困难，由于编程能力有限，我对一些函数的用法不甚了解，在查找资料后才恍然大悟。我相信这是一个不断提高能力的过程。实验2利用FFT计算线性卷积一、实验目的1.掌握利用FFT计算线性卷积的原理及具体实现方法。2.加深理解重叠相加法和重叠保留法。3.考察利用FFT计算线性卷积各种方法的适用范围。二、实验原理1.线性卷积与圆周卷积设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为（3-1）的长度为L+M-1x(n)和h(n)的圆周卷积为（3-2）圆周卷积与线性卷积相等而不产生交叠的必要条件为(3-3)圆周卷积定理：根据DFT性质，x(n)和h(n)的N点圆周卷积的DFT等于它们的DFT的乘积：(3-4)2.快速卷积快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：(1)必须选择；为了能使用基-2算法，要求。采用补零的办法使得x(n)和h(n)的长度均为N。(2)计算x(n)和h(n)的N点FFT。(3)组成乘积(4)利用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)。4、可能得到的MATLAB函数实验中FFT运算可采用MATLAB中提供的函数fft来实现。三、实验内容假设要计算序列x(n)=u(n)-u(n-L),0≤n≤L和h(n)=cos(0.2πn),0≤n≤M的线性卷积完成以下实验内容。1.设L=M,根据线性卷积的表达式和快速卷积的原理分别编程实现计算两个序列线性卷积的方法,比较当序列长度分别为8,16,32,64,256,512,1024时两种方法计算线性卷积所需时间。2．当L=2048且M=256时比较直接计算线性卷积和快速卷积所需的时间,进一步考察当L=4096且M=256时两种算法所需的时间。3.编程实现利用重叠相加法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。4.编程实现利用重叠保留法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。题目一：快速卷积法：L=input('L=');M=input('M=');x=1;h=1;fori=1:L-1x=[1x];h=[hcos(0.2*pi*i)];endNum=pow2(nextpow2(M+L-1));ticxs=fft(x,Num);hs=fft(h,Num);ys=xs.*hs;y=ifft(ys,Num);toc结果：L=8M=8时间已过0.004177秒。L=16M=16时间已过0.007234秒。L=32M=32时间已过0.012298秒。L=64M=64时间已过0.016842秒。L=256M=256时间已过0.018824秒。L=512M=512时间已过0.020002秒。L=1024M=1024时间已过0.024016秒。线性卷积法：L=input('L=');M=input('M=');n=0:L;x=heaviside(n)-heaviside(n-L);m=0:M;h=cos(0.2*pi*m);ticy=conv(x,h);toc结果：L=8;M=8;时间已过0.002205秒。L=16;M=16;时间已过0.008691秒。L=32;M=32;时间已过0.014703秒。L=64;M=64;时间已过0.021067秒。L=256;M=256;时间已过0.036894秒。L=512;M=512;时间已过0.058499秒。L=1024;M=1024;时间已过0.080003秒。结论：在L和M相等的情况下，使用快速卷积比使用线性卷积运算更快。题目二：快速卷积法：L=input('L=');M=input('M=');x=1;h=1;fori=1:L-1x=[1x];h=[hcos(0.2*pi*i)];endNum=pow2(nextpow2(M+L-1));ticxs=fft(x,Num);hs=fft(h,Num);ys=xs.*hs;y=ifft(ys,Num);toc结果：L=2048M=256时间已过0.007841秒。L=4096M=256时间已过0.009291秒。线性卷积法：L=2048;M=256;n=0:L;x=heaviside(n)-heaviside(n-L);m=0:M;h=cos(0.2*pi*m);ticy=conv(x,h);toc结果：时间已过0.050958秒。L=4096;M=256;n=0:L;x=heaviside(n)-heaviside(n-L);m=0:M;h=cos(0.2*pi*m);ticy=conv(x,h);toc时间已过0.080958秒。结论：在L=2048，M=256；L=4096，M=256;两种情况下仍然是快速卷积运算更快。题目三：L=input('L=');M=input('M=');N=M+1;x=1;h=1;forj=1:L-1x=[1x];h=[hcos(0.2*pi*j)];endx=[x,zeros(1,N-1)];t=zeros(1,M-1);H=zeros(1,L+M-1);a=floor(L/N);ticfork=0:aA=x(k*N+1:k*N+N);H1=fft(A,L+M-1);H2=fft(h,L+M-1);HH=H1.*H2;q=ifft(HH,L+M-1);H(k*N+1:k*N+M-1)=q(1:M-1)+t(1:M-1);H(k*N+M:k*N+N)=q(M:N);t(1:M-1)=q(N+1:N+M-1);endtoc结果：L=2048M=256时间已过0.003830秒。题目四：L=input('L=');M=input('M=');N=M+1;x=1;y=1;fori=1:L-1x=[1x];y=[ycos(0.2*pi*i)];endLL=N+M-1;h=[h,zeros(1,N-M)];x=[zeros(1,M-1),x,zeros(1,N-1)];a=floor((L+M-2)/(LL))+1;Y=zeros(1,N);ticfork=0:a-1xk=x(k*LL+1:k*LL+N);b=fft(xk,N);C=fft(h,N);Z=b.*C;Y(k+1,:)=ifft(Z,N);endtoc结果：L=2048M=256时间已过0.007472秒。结论：L=2048，M=256时使用重叠相加法和重叠保留法比使用快速卷积的方法运算更快一点。实验总结本实验研究了利用FFT计算线性卷积的方法。通过比较线性卷积和快速卷积两种方法，我深深地体会到了快速卷积的优越性，它可以大大提高计算效率。快速卷积的核心思想是“分而治之”，通过利用旋转因子的对称性和周期性合并一些运算项，并且可以使长序列的DFT分更小点的DFT.尽管实验的结果是对课堂已有结论的证明，但我仍在实践中找到探索发现的乐趣。实验3IIR数字滤波器设计一、实验目的1.掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。二、实验原理1、基本原理（1）从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应，h(n)等于的取样值。（2）从频率响应出发，直接使数字滤波器的频率响应逼近模拟滤波器的频率响应，进而求得数字滤波器的系统函数H（z）。2、变换方法脉冲响应不变法思路：(1)将进行部分分式展开(2)对进行拉式变换(3)对时域采样得到h(n)(4)对h(n)进行z变换双线性变换法思路：通过正切变换Ω=tan(Ω1T/2)将s平面的整个虚轴jΩ压缩到s1平面jΩ1轴上的-jπ/T到jπ/T段上。通过z=es1T变换将s1平面映射到z平面上，得到s平面和z平面的映射关系为：s=(1-z-1)/(1+z-1)一般来说为了调节模拟频带与数字频带之间的关系，可引入常数c,使映射关系为s=c(1-z-1)/(1+z-1)。设计时可采用模拟频率特性与数字频率特性在低频出有确切对应关系的常数c=2/T。设计步骤(1)确定数字滤波器的性能指标。(2)将数字滤波器频率指标转换成响应的模拟滤波器频率指标(3)根据指标，，和设计模拟滤波器。(4)将展成部分分式形式。(5)把模拟极点转换成数字极点，得到数字滤波器。可见至H(z)间的变换关系为方法1：利用residue函数和residuez函数实现脉冲响应不变变换法，实用方法如下：[r,p,k]=residue(b,a)[b,a]=residue(r,p,k)实现多项式形式和部分分式形式之间的装换[r,p,k]=residuez(b,a)[b,a]=residuez(r,p,k)实现多项式形式和部分分式形式之间的转换方法2：（1）matlab中提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用如下：[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系数向量bz和az，fs为采样频率（默认为1）。[bz,az]=impinvar(b,a)采样频率默认为1的情况下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。（2）matlab中提供了bilinear函数采用双线性变换法实现模拟滤波器到数字滤波器的变换，其使用如下：[bz,az]=bilinear(b,a,fs)采用双线性变换法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系数向量bz和az，fs为采样频率。实验内容设采样设采样频率为fs=10kHz，设计数字低通滤波器，满足如下指标通带截止频率：fp=1kHz，通带波动：Rp=1dB阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器，并分别结合脉冲响应不变法和双线性变换法进行设计。结合实验结果，分别讨论采用上述设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。代码及图形：脉冲响应不变法：巴特沃斯型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=buttord(Wp,Ws,Rp,As,'s');[ba]=butter(N,OmegaC,'s');[bzaz]=impinvar(b,a,10000);w=[0:500]*pi/500;[Hw]=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');切比雪夫Ⅰ型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=cheb1ord(Wp,Ws,Rp,As,'s');[b,a]=cheby1(N,Rp,OmegaC,'s');[bzaz]=impinvar(b,a,10000);w=[0:500]*pi/500;[Hw]=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');切比雪夫Ⅱ型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=cheb2ord(Wp,Ws,Rp,As,'s');[b,a]=cheby2(N,As,OmegaC,'s');[bzaz]=impinvar(b,a,10000);w=[0:500]*pi/500;[Hw]=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');椭圆型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=ellipord(Wp,Ws,Rp,As,'s');[b,a]=ellip(N,Rp,As,OmegaC,'s');[bzaz]=impinvar(b,a,10000);w=[0:500]*pi/500;[Hw]=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');双线性变换法：巴特沃斯型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=buttord(Wp,Ws,Rp,As,'s');[b,a]=butter(N,OmegaC,'s');[bzaz]=bilinear(b,a,10000);w=[0:500]*pi/500;[Hw]=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');切比雪夫Ⅰ型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=cheb1ord(Wp,Ws,Rp,As,'s');[b,a]=cheby1(N,Rp,OmegaC,'s');[bzaz]=bilinear(b,a,10000);w=[0:500]*pi/500;[Hw]=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');切比雪夫Ⅱ型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=cheb2ord(Wp,Ws,Rp,As,'s');[b,a]=cheby2(N,As,OmegaC,'s');[bzaz]=bilinear(b,a,10000);w=[0:500]*pi/500;[Hw]=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');椭圆型：Wp=1000*2*pi;Ws=1500*2*pi;Rp=1;As=15;[NOmegaC]=ellipord(Wp,Ws,Rp,As,'s');[b,a]=ellip(N,Rp,As,OmegaC,'s');[bzaz]=bilinear(b,a,10000);w=[0:500]*pi/500;[Hw]=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到1kHz时满足通带波动1dB；在1.5kHz满足阻带衰减15dB。切比雪夫I型低通滤波器满足题目要求：在0到1kHz时满足通带波动1dB；在1.5kHz满足阻带衰减15dB。切比雪夫II型低通滤波器不满足题目要求，通带波动明显大于1dB，阻带衰减明显小于15dB。椭圆型低通滤波器不满足题目要求，通带波动明显大于1dB，阻带衰减明显小于15dB。优点：频率坐标的变换是线性的：。因此，如果模拟滤波器的频响是限带于折叠频率以内的话，则通过变换后所得的数字滤波器的频响可以不失真的反映原响应与频率的关系。缺点：脉冲响应不变法的最大缺点是有频谱的周期性延拓效应，因此只能用于限带的频响特性。对于高通、带阻滤波器，由于它们在高频部分不衰减，当一定要追求频率线性关系而采用脉冲响应不变法时，必须先对模拟高通和带阻滤波器加一层保护滤波器，滤除高于折叠频率以上的频率，然后再转变为数字滤波器以避免混叠失真。双线性变换法：巴特沃斯型低通滤波器满足题目要求：在0到1kHz时满足通带波动1dB；在1.5kHz满足阻带衰减15dB。切比雪夫I型低通滤波器满足题目要求：在0到1kHz时满足通带波动1dB；在1.5kHz满足阻带衰减15dB。切比雪夫II型低通滤波器不满足题目要求，椭圆型低通滤波器不满足题目要求。优点：消除了脉冲响应不变法所固有的混叠误差。缺点：频率Ω与w间的非线性。这种非线性关系要求被变换的连续时间系统的幅度响应必须是分段常数型的，不然所映射出的数字频率响应相对于原来的模拟频率响应会产生变形。解决双线性变换中的频率非线性关系的方法有：在低频段，模拟与数字两滤波器的频率关系处于近似的线性范围内，故可忽略非线性影响；或者采用补偿的方法，也可称为预畸的方法。四、实验总结本实验研究了用脉冲响应不变法和双线性变换法设计IIR滤波器的问题。在这次实验中我对两种方法的优缺点有了更加熟练的掌握，并且比较分析了巴特沃斯低通滤波器、切比雪夫I型低通滤波器、切比雪夫II型低通滤波器、椭圆低通滤波器的不同之处。理论知识的学习在实践中得到巩固。在分析实验中的问题时，（比如分析脉冲响应不变法的优缺点），我才体会到老师上课所反复讲的这方面的知识在实际设计时产生的结果。因为有些东西只靠做题时无法深刻理解的，只有具体操作过才能明白其中含义。这次实验让我明白了DSP的学习不能生搬硬套，而应该自己多动手，多时间，方可达到把这门课学精、学通的境界。实验4FIR数字滤波器设计一、实验目的 掌握窗函数法和频率取样法设计FIR数字滤波器的原理及具体方法实验原理窗函数法： 1.基本原理 窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器，然后用窗函数截取它的单位脉冲响应，得到线性相位和因果的FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器，使设计的滤波器的单位脉冲响应逼近理想滤波器的单位脉冲响应。 2.设计步骤 （1）给定理想滤波器的频率响应，在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽为的低通滤波器由下式给定：其中为采样延迟，其作用是为了得到一个因果系统。（2）确定这个滤波器的单位脉冲响应为了得到一个长度为N的因果的线性相位FIR滤波器，我们令（3）用窗函数截取得到所设计FIR数字滤波器：3.窗函数的选择常用的窗函数有矩形（Rectangular）窗，汉宁（Hanning）窗，海明（Hamming）窗、布莱克曼（Blackman）窗、凯瑟（Kaiser）窗等。MATLAB提供了一些函数用于产生窗函数。如表7-1所示。表7-1MATLAB中产生窗函数的命令MATLAB函数窗函数MATLAB函数窗函数Boxcar矩形窗函数Blackman布莱克曼窗Hanning汉宁窗函数Kaiser凯瑟窗函数Hamming海明窗在设计过程中我们需要根据给定的滤波器技术指标，选择滤波器长N和窗函数。表7-2列出了常用窗函数的一些特性，可供设计参考。表7-2常用窗函数的特性窗函数窗函数频率特性加窗后滤波器指标旁瓣峰值dB主瓣宽度过渡带宽最小阻带衰减dB矩形窗-134π/N1.8π/N-21汉宁窗-318π/N6.2π/N-44海明窗-418π/N6.6π/N-53布莱克曼窗-5712π/N11π/N-74凯瑟窗是一种广泛在实际中广泛应用的窗函数，它由下式给定：其中是修正的零阶贝塞尔函数，参数控制最小阻带衰减，这种窗函数对于相同的N可以提供不同的过渡带宽。由于贝塞尔函数比较复杂，这种窗函数的设计方程很难推导，然而幸运的是，有一些经验设计方程可以直接使用。已知给定的指标，滤波器长度N和凯瑟窗参数可以按如下凯瑟窗方程给出过渡带带宽：频率取样法：基本原理：频率取样法从频域出发，把理想的滤波器Hd(ejw)等间隔采样得到Hd(k）,将Hd(k）作为实际设计滤波器的H（k）。H(k）=Hd(k）=H(ejw)|w=2πk/Nk=0,1,…N-1得到H(k）以后可以由H(k）唯一确定滤波器的单位脉冲响应h(n),H(ejw)也可以由H(k）求得：h(n)=IDFT[H（k）],H(ejw)=如果我们设计的是线性相位FTR滤波器，则H（k)的幅度和相位一定满足线性相位滤波器的约束条件。我们将H(k)表示成如下形式：H(k)=|H(k)|ejθ(k)=Hr(k)ejθ(k)当h(n)为实数，则H(k)=H*(N-k),由此得到：Hr(k)=Hr(N-k)三、实验内容1、设计一个数字低通FIR滤波器，其技术指标如下：分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。矩形窗wp=0.2*pi;wst=0.3*pi;rp=0.25;as=50;tr_width=wst-wp;N=ceil(1.8*pi/tr_width)+1;n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_boxcar=boxcar(N)';%产生窗函数，并对hd(n)加窗得到h(n)h=hd.*w_boxcar;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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>>汉宁窗：wp=0.2*pi;wst=0.3*pi;rp=0.25;as=50;tr_width=wst-wp;N=ceil(6.2*pi/tr_width)+1;n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_hanning=hanning(N)';%产生窗函数，并对hd(n)加窗得到h(n)h=hd.*w_hanning;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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>>海明窗：wp=0.2*pi;wst=0.3*pi;rp=0.25;as=50;tr_width=wst-wp;N=ceil(6.6*pi/tr_width)+1;n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_hamming=hamming(N)';%产生窗函数，并对hd(n)加窗得到h(n)h=hd.*w_hamming;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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>>布莱克曼窗：wp=0.2*pi;wst=0.3*pi;rp=0.25;as=50;tr_width=wst-wp;N=ceil(11*pi/tr_width)+1;n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_blackman=blackman(N)';%产生窗函数，并对hd(n)加窗得到h(n)h=hd.*w_blackman;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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>>凯瑟窗：wp=0.2*pi;wst=0.3*pi;rp=0.25;as=50;tr_width=wst-wp;N=ceil((as-7.95)/(2.285*tr_width))+1;n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc((wc/pi)*(n-alpha));w_kaiser=kaiser(N)';%产生窗函数，并对hd(n)加窗得到h(n)h=hd.*w_kaiser;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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>>2、设计一个数字带通FIR滤波器，其技术指标如下：下阻带边缘：下通带边缘：上通带边缘：上阻带边缘：用海明窗设计：wp1=0.35*pi;wst1=0.2*pi;wp2=0.65*pi;wst2=0.8*pi;tr_width=wp1-wst1;N=ceil(6.6*pi/tr_width)+1;n=0:(N-1);alpha=(N-1)/2;wc1=(wp1+wst1)/2;wc2=(wp2+wst2)/2;hd=(wc2/pi)*sinc((wc2/pi)*(n-alpha))-(wc1/pi)*sinc((wc1/pi)*(n-alpha));w_hamming=hamming(N)';h=hd.*w_hamming;subplot(221);stem(n,hd,'filled');axistight;xlabel('n');ylabel('hd(n)');[Hr,w1]=zerophase(h);subplot(222);plot(w1/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');gridon3.采用频率取样法设计FIR数字低通滤波器，满足以下指标：Wp=0.2π,Rp=0.25dBWst=0.3π,As=50dB取N=20，过渡带没有样本取N=40，过渡带有一个样本，T=0.39取N=60，过渡带有两个样本，T1=0.5925，T2=0.1099代码及图形：（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);w=[0:500]*pi/500;H=freqz(h,l,w);[Hr,wr]=zerophase(h);subplot(221);plot(wdl,Hdr,wl(1:11)/pi,Hrs(1:11),'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(wr/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]);xlabel('\omega(\pi)');ylabel('Hr(ω)');subplot(224);plot(w/pi,20*log10((abs(H)/max(abs(H)))));axis([0,1,-50,5]);grid;xlabel('\omega(\pi)');ylabel('dB');（2）取N=40，过渡带有一个样本，T=0.39N=40;alpha=(N-1)/2;l=0:N-1;wl=(2*pi/N)*l;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];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;H=freqz(h,l,w);[Hr,wr]=zerophase(h);subplot(221);plot(wdl,Hdr,wl(1:21)/pi,Hrs(1:21),'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(wr/pi,Hr,wl(1:21