北理工数字信号处理实验报告

1、北京理工人学ffrynry* bixfitule of 丁/噌y本科实验报告实验名称：数字信号处理实验课程名称：数字信号处理实验实验时间：周二810节任课教师：实验地点：4 423实验教师：实验类型： 原理验证 综合设计 自主创新学生姓名：任科飞学/班级：1120121159/05111251组号：学 院：信息与电子同组搭档：专 业：信息工程成绩：信息与电子学院KMM。 INrfiMliUThQM M0 flUCTMHICt欢迎下载8实验1利用DFT分析信号频谱一、实验目的1 .加深对DFT原理的理解。2 .应用DFT分析信号频谱。3 .深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及

2、解决办法。二、实验原理1、DFT和DTFT的关系有限长序列x(n)的离散时间傅里叶变换X(ej )在频率区间(0个等分点x(0), x(1),- x(k), - x(N 1)上的N个取样值可以由下式表示:(2 1)DTFT在N个等iN 1j2 knX(ej )| 2 k x(n)e N X(k) 0 k N 1k 0由上式可知，序列x(n)的n点dftX (k),实际上就是x(n)序列的间隔频率点X(0), X(1), X(k)X(N 1)上样本X(k)。2、利用DFT求DTFT方法1：由X(k)恢复出X(ej )的方法如图2.1所示:mi卬FT*工Q?)DTFT图2.1.由N点DFT恢复频谱

3、 DTFT的流程由图2.1所示流程图可知：X(ej )x(n)eX(k)WN kn(22)由式2-2可以得到x(ej )NX(k)k 1(23)其中(x)为内插函数sin(N 2)?eN sin(2 4)方法2:然而在实际 MATLAB十算中，上诉插值公式不见得是最好的方法。由于DFT是DTFT的取样值，其相邻的两个频率样本点的间距为2/N ,所以如果我们增加数据的长度N,使得得到的DFT谱线就更加精细， 其包络就越接近 DTFT的结果，这样可以利用 DFT来近 似计算DTFT如果没有更多的数据，可以通过补零来增加数据长度。3、利用DFT分析连续时间信号的频谱采用计算机分析连续时间信号的频谱，

4、第一步就是把连续时间信号离散化，这里需要进行连个操作：一是采样，二是截断。对于连续非周期信号 X (t),按采样间隔T进行采样，截取长度为 M,那么M 1X (j ) x (t)e j tdt T x (nT)e j nT (2 5) n 0对X (j )进彳r n点的频率采样，得到M 1jknX (j )|2 T x (nT)e N TXm (k)(2 6)kNT n 0因此，可以将利用 DFT分析连续非周期信号频谱的步骤归纳如下：(1)确定时域采样间隔 t,得到离散序列 x(n);(2)确定截取长度 m得到m点离散序列xM (n) x(n)w(n),这里的w(n)为窗函 数。(3)确定频域

5、采样点数 N,要求N M。(4)利用FFT计算离散序列的N点DFT,得到XM (k)。(5)根据式(2-6)由XM (k)计算X (j )采样点的近似值。采用上诉方法计算的频谱，需要注意如下三点问题：(1)频谱混叠。如果不满足采样定理的条件，频谱会很出现混叠误差。对于频谱无限宽 的信号，应考虑覆盖大部分主要频率的范围。(2)栅栏效应和频谱分辨率。使用 DFT计算频谱，得到的结果只是 N个频谱样本值，样 本值之间的频谱是未知的，就像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成正比，提高频谱分辨率，就要增加记录时间。(3)频谱泄露。对于信号截断会把窗函数的频谱会引入到信号频谱中，

6、造成频谱泄露。解决这问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。因此，要合理选取采样间隔和截取长度，必要时还需考虑适当的窗。对于连续周期信号， 我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上诉方法近似计算。4、可能用到MATLAB®数与代码实验中的DFT运算可以采用 MATLA讣提供的FFT来实现。DTFT可以利用MATLAB巨阵运算的方法进行计算。三、实验内容1. x(n) 2, 1,1,1,完成如下要求：(1)计算其DTF1并画出,区间的波形。(2)计算4点DFT,并把结果显示在(1)所画的图形中。(3)对x(n)补零，计算6

7、4点DFT,并显示结果。(4)是否可以由DFT计算DTFT如果可以，请编程实现。2 .考察序列x(n) cos(0.48 n) cos(0.52 n)(1) 0 n 10时，用DFT估atx(n)的频谱；将x(n)补零加长到长度为100点序列用DFT估af x(n)的频谱。要求画出相应波形。 0 n 100时，用DFT估方t x(n)的频谱，并画出波形。3 .已知信号 x(t) 0.15sin(2f1t) sin(2 f2t) 0.1sin(2 f3t)，其中 f11Hz ,f2 2Hz, f3 3Hz o从x(t)的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正

8、弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。4 .利用DFT近似分析连续时间信号的频谱(幅度谱)。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。四、实验代码及实验结果第一题实验代码x = 2,-1, 1, 1;n = 0 : 3;n1 = 0 : 63;z = zeros(1,60);x1 = x z;y = fft(x);y1 = fft(x1);w = -pi : 0.01 * pi : pi;X = x * exp(-j * n' * w);subplot(211);plot(w, abs(X);xlabel( '

9、;Omega/pi' );title( 'Magnitude' );hold on;stem(n, y, 'filled' );axis tight ;subplot(212);plot(w, angle(X)/pi);xlabel( 'Omega/pi' );title( 'Phase');axis tight ;figure;stem(n1, y1, 'filled' );实验结果MagnitudePhaseMagnitude分析可以由DFT计算DTFT)增加补零的个数，提高采样密度，可以由DFT近似计

10、算DTFT第二题实验代码n = 0 : 10;x = cos(0.48 * pi * n) + cos(0.52 * pi * n);n1 = 0 : 99;n2 = 0 : 100;z = zeros(1, 89);x1 = x z;x2 = cos(0.48 * pi * n2) + cos(0.52 * pi * n2);y = fft(x);y1 = fft(x1);y2 = fft(x2);subplot(211);stem(n, y, 'filled' );title( '0<=n<=10');subplot(212);stem(n1,

11、y1, 'filled' );title( '0<=n<99');figure;stem(n2, y2, 'filled' );title( '0<=n<=100');实验结果10-60Ckaul 口52589101(3)分析但是补零不能够增加分辨可以通过增大截取长度和增加补零的个数来提高频谱分辨率,力。第三题实验代码> > n=0:39;> > x1=0.15*sin(2*pi*n)+sin(4*pi*n)-0.1*sin(6*pi*n);> > X1=fft(x1,4

12、0);> > subplot(311)> > stem(n,abs(X1);>> x2=0.15*sin(2*pi*(1/6)*n)+sin(4*pi*(1/6)*n)-0.1*sin(6*pi*(1/6)*n);> > X2=fft(x2,40);> > subplot(312)> > stem(n,abs(X2);>> x3=0.15*sin(2*pi*0.1*n)+sin(4*pi*0.1*n)-0.1*sin(6*pi*0.1*n);> > X3=f3(x3,40);> > su

13、bplot(313)> > stem(n,abs(X3);实验结果欢迎下载分析：减小采样周期，可以使频谱分辨率提高。第四题 实验代码 n=0:100; x=exp(-0.1*n);X1=fft(x); k=0:100; stem(k,abs(X1); xlabel('k');实验结果实验代码n=0:2:100;x=exp(-0.1*n);X2=fft(x);k=0:50;stem(k,abs(X2); xlabel('k');实验结果实验代码>> n=0:50; x=exp(-0.1*n); X3=fft(x);k=0:50; stem(

14、k,abs(X3); xlabel('k');实验结果实验代码>> n=0:2:50;x=exp(-0.1*n);X4=fft(x); k=0:25; stem(k,abs(X4); xlabel('k');实验结果实验代码 >> n=0:150; x=exp(-0.1*n);X5=fft(x); k=0:150; stem(k,abs(X5); xlabel('k');实验结果实验代码>> n=0:2:150;x=exp(-0.1*n);X5=fft(x);k=0:75;stem(k,abs(X5);xlab

15、el('k');实验结果欢迎下载45分析：采样区间可以为 0到100,采样间隔为2。五、心得与体会通过本次实验，我们掌握并加深对DFT原理的理解并且学会应用DFT分析信号频谱，在此基础上利用 DFT分析信号频谱的原理，掌握了分析现实现象的步骤及办法。并且，通过这次的实验，我对DFT的原理确实有了更深的了解，对DFT和DTFT之间的关系也有了更多的认识：序列x(n)的N点DFT X(K),实际上就是x(n)序列的DTFT在N个等间隔频率点上的样本 X(K)。所以，我们可以通过增加数据的长度N或通过补零来使 DFT更加接近DTFT的结果。这样就可以利用DFT近似计算DTFT实验2利

16、用FFT计算线性卷积(选作)一、实验目的1 .掌握利用FFT计算线性卷积的原理及具体实现方法。2 .加深理解重叠相加法和重叠保留法。3 .考察利用FFT计算线性卷积各种方法的适用范围。二、实验原理1.线性卷积与圆周卷积设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为y1(n) x(n)* h(n) x(m)h(n m) (3-1) my1(n)的长度为L+M-1x(n)和h(n)的圆周卷积为 N 1y(n) x(n) e h(n) N x(m)h(n m)N Rn (n)(3-2)m 0圆周卷积与线性卷积相等而不产生交叠的必要条件为二二 - 一-(3-3)圆周卷积定理：

17、根据DFT性质，x(n)和h(n)的N点圆周卷积的 DFT等于它们的DFT的乘积:(3-4)DFTx(n)Oh(iL) = X(k)H(k)2 .快速卷积快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：必须选择N三L+M L；为了能使用基-2算法，要求N = 2 ：。采用补零的办法使得x(n)和h(n)的长度均为 N。 (2)计算 x(n)和 h(n)的 N 点 FFT。x(n)FFTX(k)h(n)FFTH(k)(3)组成乘积Y(k) X(k)H (k)(4)利用IFFT计算Y(k)的IDFT ,得到线性卷积 y(n)Y(k)

18、IFFT y(n)3 .分段卷积我们考察单位取样响应为 h(n)的线性系统，输入为 x(n),输出为y(n),则y(n) = x(n) * h(n)当输入序列x(n)极长时，如果要等x(n)全部集齐时再开始进行卷积，会使输出有较大延时；如果序列太长，需要大量存储单元。为此，我们把 x(n)分段，为别求出每段的卷积，合 在一起得到最后的总输出。这称为分段卷积。分段卷积可以细分为重叠保留法和重叠相加法。重叠保留法：设x(n)的长度为N, h(n)的长度为M。把序列x(n)分成多段N点序列再(口)，每段雨前一段重写 M-1个样本。并在第一个输入段前面补 M-1个零。计算每一段与 h(n)的 圆周卷积

19、，其结果中前 M-1个不等与线性卷积，应当舍去，只保留后面 N-M+1个正确的输 出样本，把它们合起来得到总的输出。利用FFT实现重叠保留法的步骤如下：(1)在x(n)前面填充M-1个零，扩大以后的序列为x(n) = £00 0,(2)将x(n)分为若干段N点子段，设L=N-M+1为每一段的有效长度，则第i段的数据为：%(门)=x(m) iL< m <1L + N -1, i > 0,0 < n < N - 1(3)计算每一段与h(n)的N点圆周卷积，利用FFT计算圆周卷积xi(n)FFT Xi(k)h(n) FFT H(k)Y(k) FFT Xi(k)

20、H(k)Y(k) IFFT y(n)(4)舍去每一段卷积结果的前M-1个样本，连接剩下的样本得到卷积结果y(n)o重叠相加法：设h(n)长度为M,将信号x(n)分解成长为L的子段。以羽(可表示没断信号 则：x(n)x/n)i 0x(n + iL),0 < n < L 1XlCn) = l 0其他 x(a) * h(n) = 2 看(力事 h(n)i=0每一段卷积居(口)的长度为L+M-1 ,所以在做求和时，相邻两段序列由M-1个样本重叠，即前一段的最后 M-1个样本和下一段前 M-1个样本序列重叠，这个重叠部分相加，再与不重 叠的部分共同组成 y(n)。利用FFT实现重叠保留法的步

21、骤如下：(1)将x(n)分为若干L点子段片(口)。(2)计算每一段与h(n)的卷积，根据快速卷积法利用FFT计算卷积。(3)将各段相加，得到输出 y(n)。y(n)yi(n iL)i 04、可能得到的 MATLAESi数实验中FFT运算可采用MATLAB提供白函数fft来实现。三、实验内容假设要计算序列 x(n)=u(n)-u(n-L),0 w nW L和h(n)=cos(0.2兀n),0 w nw M的线性卷积 完成以下实验内容。1 .设L=M,根据线性卷积的表达式和快速卷积的原理分别编程实现计算两个序列线性卷积的 方法，比较当序列长度分别为8,16,32,64,256,512,1024时两

22、种方法计算线性卷积所需时间。2 .当L=2048且M=256时比较直接计算线性卷积和快速卷积所需的时间，进一步考察当L=4096且M=256时两种算法所需的时间。3 .编程实现利用重叠相加法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与2题的结果进行比较。4 .编程实现利用重叠保留法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与2题的结果进行比较。四、实验代码及实验结果第一题实验代码L = input( 'length=' );n1 = 0 : 2 * L - 1;n2 = 0 : L - 1;x = ones(1,

23、L);h = cos(0.2 * pi * n2);z = zeros(1, L - 1);x1 = x z;hl = h z;tic;y1 = conv(x, h);toc;tic;X = fft(x1);H = fft(h1);Y = X .* H;y = ifft(Y);toc;实验结果(第一行为线性卷积，第二行为FFT)L=8M=8时间已过0.232076 秒。时间已过0.087391 秒。L=16M=16时间已过0.000332 秒。时间已过0.000098 秒。L=32M=32时间已过0.000150 秒。时间已过0.000051 秒。L=64M=64时间已过0.289400 秒。

24、时间已过0.000758 秒。L=256M=256时间已过0.000238 秒。时间已过0.000169 秒。L=512M=512时间已过0.000455 秒。时间已过0.000235 秒。L=1024M=1024时间已过0.001132 秒。时间已过0.000922 秒。分析：可见在在相同长度下，快速卷积比线性卷积快差不多一倍的时间。第二题实验代码实验结果：（第一行为线性卷积，第二行为FFT）L=2048M=256时间已过0.012061秒。时间已过0.000581秒。分析：快速卷积比线性卷积快。第三题Lx=input('L=');M=input('M=');

25、N=M+1;x=1;y=1;for i=1:L-1x=1 x;y=y cos(0.2*pi*i);endx=x,zeros(1,N-1);t=zeros(1,M-1);Y=zeros(1,Lx+M-1);a=floor(Lx/N);ticfor k=0:aA=x(k*N+1:k*N+N);y1=fft(A,Lx+M-1);y2=fft(h,Lx+M-1);y3=y1.*y2;q=ifft(y3,Lx+M-1);Y(k*N+1:k*N+M-1)=q(1:M-1)+t(1:M-1);Y(k*N+M:k*N+N尸q(M:N);t(1:M-1)=q(N+1:N+M-1);endtoc实验结果：L=20

26、48M=256时间已过0.007656秒。分析：较第2题结果快。第四题实验代码：Lx=input('L=');M=input('M=');N=M+1;x=1;y=1;for i=1:L-1x=1 x;y=y cos(0.2*pi*i);endM1=M-1;L=N+M1;h=h,zeros(1,N-M);x=zeros(1,M1),x,zeros(1,N-1);a=floor (Lx+M1-1)/(L)+1;Y=zeros(1,N);ticfor k=0:a-1xk=x(k*L+1:k*L+N);b=fft(xk,N);C=fft(h,N);Z=b.*C;Y(k+

27、1,:)=ifft(Z,N);endtocL=2048M=256实验结果：L=2048M=256时间已过0.002848秒。分析：较第2题结果快。五、心得与体会本次实验要求我们掌握利用FFT计算线性卷积的原理及具体实现方法，通过实验加深理解重叠相加法和重叠保留法并考察利用FFT计算线性卷积各种方法的适用范围。本次实验让我切实看到了 FFT算法的高效性及对于庞大的数据的处理实时性，同时对重叠保留法、重叠相加法也有了更深的认识，对以后的实验有了很好的理论基础，受益颇多。实验3 IIR数字滤波器设计一、实验目的1 .掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。2 .加深

28、理解数字滤波器和模拟滤波器之间的技术指标转化。3 .掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。二、实验内容1 .设采样频率为fs =4kHz ,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器， 其3dB截止频率为fc=1kHz。2 .设采样频率为fs=10kHz,设计数字低通滤波器，满足如下指标通带截止频率：fp=1kHz ,通带波动：Rp=1dB阻带截止频率：fst=1.5kHz ,阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫 I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应 不变法进行设计。结合实验结果，分别讨论采用上述设计的数字滤波器是否都

29、能满足给定指 标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。四、实验代码及实验结果实验代码：巴特沃斯：fp = 1000；fs = 1500;f = 10000;Wp = 2 * pi * fp;Ws = 2 * pi * fs;wp = 2 *f * tan(2 * pi * fp / (2 * f);ws = 2 * f * tan(2 * pi * fs / (2 * f);Rp = 1;As = 15;N1 = ceil(log10(10A(Rp/10)-1)/(10A(As/10)-1)/(2*log10(Wp/Ws);N2 = ceil(log10(10A(Rp

30、/10)-1)/(10A(As/10)-1)/(2*log10(wp/ws);OmegaCI = Wp / (10 A (Rp / 10) - 1) A (1 / (2 * N1);OmegaC2 = wp / (10 A (Rp / 10) - 1) A (1 / (2 * N2);z, p, k = buttap(N1);b1 = k * OmegaC1 A N1;a1 = poly(p * OmegaC1);H, w = freqs(b1, a1);subplot(221);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)' );y

31、label( subplot(223);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)' );ylabel( b2, a2 = butter(N2, OmegaC2, H, w = freqs(b2, a2);subplot(222);plot(w/pi, abs(H);axis tight ;grid on;xlabel( 'Omega(pi)' );ylabel( subplot(224);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)' );ylab

32、el(实验结果：'|H(jOmega)|' );title(? ? ?1?3?3?i o |2? ±?');'Phase of H(jOmega)|' );title('s')；'|H(jOmega)|' );title(Phase of H(jOmega)|? ? ?1?3?3?i o |2? ±?');? ? ?1?D?t ?'););title( ' ° i i ?1?D?t ?');巴特波斯脉冲脉冲响应不变巴特泱斯双线性变换巴特沃斯脉冲脉冲响应不变: 4

33、123口 x 104巴特法斯双线性变换_6?-0切比雪夫I型：fp = 1000；fs = 1500;f = 10000;Wp = 2 * pi * fp;Ws = 2 * pi * fs;wp = 2 *f * tan(2 * pi * fp / (2 * f);ws = 2 * f * tan(2 * pi * fs / (2 * f);Rp = 1;As = 15;e = (10 人(Rp / 10) - 1)人 0.5;N1 = ceil(acos(As * As) - 1) A 0.5 / e) / acos(Ws / Wp);N2 = ceil(acos(As * As) - 1)

34、 A 0.5 / e) / acos(ws / wp);b1, a1 = cheby1(N1, Rp, Ws,'s');H, w = freqs(b1, a1);subplot(221);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)' subplot(223);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)');ylabel( '|H(jOmega)|');ylabel( 'Phase of);title( '?D ±

35、; e ? - o 1?3?3?i o |2? 士？);H(jOmega)|' );title('?D ± e ? - o 1?3?3?i o |2? 士？);b2, a2 = cheby1(N2, Rp, ws,'s');H, w = freqs(b2, a2);subplot(222);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)' );ylabel( '|H(jOmega)|' );title( '?D ± e ? - 6 1?D»?'

36、);subplot(224);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)' );ylabel( 'Phase of H(jOmega)|');title( '?D ± e ? - 6 1?D?± ?');实验结果：切比雪夫1脉冲.脉冲响应不变切比雪夫1双线性变换§ 0.5切比雪夫1脉冲脉冲响应不变41231口N& 1口4切比雪夫1双线性变换2 0 2 OOJSW2 0-2 = 一cdh 一口用 Ed.4 ,；01234«in4-4023 d* 父 1

37、04切比雪夫n型：fp = 1000;fs = 1500;f = 10000;Wp = 2 * pi * fp;Ws = 2 * pi * fs;wp = 2 *f * tan(2 * pi * fp / (2 * f);ws = 2 * f * tan(2 * pi * fs / (2 * f);Rp = 1;As = 15;e = (10 人(Rp / 10) - 1)人 0.5;K = Wp / Ws;k1 = e /(As * As) - 1)人 0.5;k = wp / ws;N1 = ceil(ellipke(K) * ellipke(1 - k1) A 0.5) / (ellip

38、ke(kl) * ellipke(1 - K) A 0.5);'s');N2 = ceil(ellipke(k) * ellipke(1 - k1) A 0.5) / (ellipke(kl) * ellipke(1 - k) A 0.5);b1, a1 = ellip(N1, Rp, As, Ws, H, w = freqs(b1, a1);subplot(221);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)');ylabel('|H(jOmega)|);title('1 ?2?3?3?i o |2

39、? ±？');subplot(223);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)');ylabel('Phase of H(jOmega)|');title( ' i ?2?3?3?i o |2? ±?');b2, a2 = ellip(N2, Rp, As, ws,'s');H, w = freqs(b2, a2);subplot(222);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)');

40、ylabel('|H(jOmega)|);title('i?2?D?t ?');subplot(224);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)');ylabel('Phase of H(jOmega)|');title( ' i?2?D?t ?');实验结果:切比雪夫2脉冲脉冲响应不变切比雪夫2脉冲脉冲响应不变切比雪夫2双线性变换一sDrocliswcl椭圆型:fp = 1000;fs = 1500;f = 10000;Wp = 2 * pi * fp;Ws = 2

41、* pi * fs;wp = 2 *f * tan(2 * pi * fp / (2 * f);ws = 2 * f * tan(2 * pi * fs / (2 * f);Rp = 1;As = 15;e = (10 人(Rp / 10) - 1)人 0.5;K = Wp / Ws;k1 = e /(As * As) - 1)人 0.5;k = wp / ws;N1 = ceil(ellipke(K) * ellipke(1 - k1) A 0.5) / (ellipke(k1) * ellipke(1 - K) A 0.5);N2 = ceil(ellipke(k) * ellipke(1

42、 - k1) A 0.5) / (ellipke(k1) * ellipke(1 - k) A 0.5);b1, a1 = ellip(N1, Rp, As, Ws,'s');H, w = freqs(b1, a1);subplot(221);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)' subplot(223);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)');ylabel( '|H(jOmega)|' );title();ylabel

43、( 'Phase of H(jOmega)|? ?2?3?3?i o |2? ±？'););title( ' i ?2?3?3?i o |2? ±?');b2, a2 = ellip(N2, Rp, As, ws,'s');H, w = freqs(b2, a2);subplot(222);plot(w/pi, abs(H);grid on;xlabel( 'Omega(pi)' subplot(224);plot(w/pi, angle(H);grid on;xlabel( 'Omega(pi)

44、9;);ylabel( '|H(jOmega)|' );title();ylabel( 'Phase of H(jOmega)|? ?2?D?t ?'););title( ' i ?2?D?t ?');实验结果:椭圆脉冲脉冲响应不娈椭圆双线性变换口 切x 104椭圆脉冲脉冲响应不变： 4123口 X 11/椭周双线性受换分析：脉冲响应不变法的优点是频率坐标的变换是线性的，因此如果模拟滤波器的频响是限带于折叠频率以内的话，通过变换后的数字滤波器的频响可以不失真地反映原响应与频率之 间的关系。但是，其最大的缺点是有频谱的周期延拓效应，会产生混叠失真。

45、故，脉冲响应 不变法只能用于限带的频响特性，如衰减特性较好的低通，或带通，而且高频衰减越大，频 响的混叠效应就越小。 双线性变换法的优点是消除了脉冲响应不变法的固有的混叠失真，但缺点是频率变换之间的非常严重的非线性，会使变换后时域上的图像产生非常严重的失真。五、心得与体会通过实验，我学习了利用脉冲响应不变法设计IIR数字滤波器的基本原理为：从时域响应出发，使数字滤波器的单位脉冲响应 h（iO模仿模拟滤波器的单位冲激响应 hJtkhOi）等于hjt）的取样值。所以设计中，应先根据数字指标转换来的模拟指标，设计模拟滤波器，再利用impinvar函数实现脉冲响应不变法模拟滤波器到数字滤波器的变换。脉

46、冲响应不变法和双线性法设计IIR数字滤波器是非常实用的技能，也会时我对 DSP这门课程有更深刻的理解和认识。本次实验最大的收获是了解了多种模拟滤波器到数字滤波器的转换方法。通过实验， 确实巩固了课本知识。完成实验报告的过程中，查阅了相关资料，不仅开拓了知识面，更加深了对本部分知识的理解。实验4频率取样法设计FIR数字滤波器、实验目的掌握频率取样法设计FIR数字滤波器的原理及具体方法二、实验原理1、基本原理频率取样法从频域出发，把理想的滤波器 Hd(ej )等间隔采样得到 Hd(k),将Hd(k)作为实际设计滤波器的 H(k)。H(k) Hd(k) H(ej )12k 0,1,L ,N 1(8-

47、1)kN得到H(k)以后可以由H (k)来确定唯一确定滤波器的单位脉冲响应h(n), H (ej )可以由H(k)求得:h(n) IDFTH(k)N 1H(ej ) H(k)(k 02Tk)(8-2、3)其中(x)为内插函数sin(N /2)Nsin( /2)(8-4 )由H(k)求得的频率响应 H(ej )将逼近Hd(ej )如果我们设计的是线性相位FIR滤波器，则H (k)的幅度和相位满足线性相位滤波器的约束条件。我们将H (k)表示为如下形式:H (k)= H(k) ej (k)Hr(k)ej (k)(8-5)当h(n)为实数，则H (k) H (N k)由此得到Hr(k) Hr(N k

48、)(8-6)即Hr(k)以k N / 2为中心偶对称。在利用线性相位条件可知，对于 1型和2型线性相位滤波器(k)0,L ,(8-7)1,L N 1对于3型和4型线性相位滤波器(k)N 1 2m2=k(N k)N 121,L N 1(8-8 )其中，一 .L"一；"4._二ER设计步骤(1)由给定的理想滤波器给出 Hr(k)和(k)。(2)由 H(k)= H (k) ej (k) H.(k)ej (k)求得 H (k)根据H(k)求得h(n)或H(ej )三、实验内容1 .设计一个数字低通 FIR滤波器，其技术指标如下：co_p=0.2 兀，R_p=0.25dBco_st=

49、0.3 兀，A_s=50dB分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果分别讨论上述方法设计的数字滤波器是否符合指标。2 .设计一个数字带通 FIR滤波器，其技术指标如下：3 .采用频率取样设计法设计FIR数字低通滤波器，满足以下指标co_p=0.2 兀，R_p=0.25dBco_st=0.3 兀，A_s=50dB(1)取N=20,过渡带没有样本。(2)取N=40,过渡带有一个样本，T=0.39。(3)取 N=60,过渡带有两个样本，T1=0.5925, T2=0.1099 。(4)分别采用上述方法设计的数字滤波器是否都能满足给定的指标要求。4.采用频率取样技术

50、设计下面的高通滤波器co _st=0.6 兀，A_s=50dBw_p=0.8 兀,R_p=1dB对于高通滤波器，N必须为奇数(或1型滤波器)。选择N=33,过渡带有两个样本，过渡带 样本最优值为 T1=0.1095 , T2=0.598。四、实验代码及实验结果第一题矩形窗：实验代码：wp = 0.2 * pi; wst = 0.3 * pi;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) * sin

51、c(wc / pi) * (n - alpha);w_boxcar = boxcar(N)'h = hd .* w_boxcar;subplot(221);stem(n, hd,'filled' );axis tight ; xlabel( 'n' ); ylabel( 'hd(n)');Hr, w1 = zerophase(h););ylabel( 'H(omega)' );subplot(222); plot(w1/pi, Hr);axis tight ; xlabel('omega八pi'subplo

52、t(223); stem(n, h,'filled' );'h(n)');axis tight ; xlabel( 'n' ); ylabel( H, w = freqz(h, 1);subplot(224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel('omega八pi'); ylabel( 'dB');grid on;实验结果：分析：阻带的衰减小于50db,所以，不符合指标。 汉宁窗：实验代码：wp = 0.2 * pi; wst = 0.3 * pi

53、;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)'h = hd .* w_hanning;subplot(221);stem(n, hd,'filled' );axis tight ; xlabel( 'n' ); ylabel( 'hd(n)&#

54、39;);Hr, w1 = zerophase(h););ylabel( 'H(omega)' );subplot(222); plot(w1/pi, Hr);axis tight ; xlabel('omega八pi'subplot(223); stem(n, h,'filled' );'h(n)');axis tight ; xlabel( 'n' ); ylabel( H, w = freqz(h, 1);subplot(224); plot(w / pi, 20 * log10(abs(H) / max(abs(H);xlabel(