《数字信号处理实验报告_完整版.pdf》由会员分享,可在线阅读,更多相关《数字信号处理实验报告_完整版.pdf(57页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、word 文档 可自由复制编辑 实验 1 利用 DFT分析信号频谱 一、实验目的 1.加深对 DFT原理的理解。2.应用 DFT分析信号的频谱。3.深刻理解利用 DFT分析信号频谱的原理,分析实现过程中出现的现象及解决方法。二、实验设备与环境 计算机、MATLAB软件环境 三、实验基础理论 1.DFT与 DTFT的关系 有限长序列 的离散时间傅里叶变换 在频率区间 的 N 个等间隔分布的点 上的 N 个取样值可以由下式表示:212/0()|()()01NjknjNk NkX ex n eX kkN 由上式可知,序列 的 N点 DFT ,实际上就是 序列的 DTFT在 N个等间隔频率点 上样本
2、。2.利用 DFT求 DTFT 方法 1:由恢复出的方法如下:由图 2.1所示流程可知:101()()()Njj nknj nNnnkX ex n eX k WeN 由上式可以得到:IDFT DTFT ()word 文档 可自由复制编辑 12()()()NjkkX eX kN 其中为内插函数 12sin(/2)()sin(/2)NjNxeN 方法 2:实际在 MATLAB计算中,上述插值运算不见得是最好的办法。由于DFT是 DTFT的取样值,其相邻两个频率样本点的间距为 2 /N,所以如果我们增加数据的长度 N,使得到的 DFT谱线就更加精细,其包络就越接近 DTFT的结果,这样就可以利用 D
3、FT计算 DTFT。如果没有更多的数据,可以通过补零来增加数据长度。3.利用 DFT分析连续信号的频谱 采用计算机分析连续时间信号的频谱,第一步就是把连续信号离散化,这里需要进行两个操作:一是采样,二是截断。对于连续时间非周期信号,按采样间隔 T 进行采样,阶段长度 M,那么:10()()()Mj tj nTaaanXjx t edtTx nT e 对进行 N 点频域采样,得到 2120()|()()MjknNaaMknNTXjTx nT eTXk 因此,可以将利用DFT分析连续非周期信号频谱的步骤归纳如下:(1)确定时域采样间隔T,得到离散序列(2)确定截取长度M,得到 M 点离散序列,这里
4、为窗函数。(3)确定频域采样点数N,要求 N M。(4)利用 FFT计算离散序列的N 点 DFT,得到.(5)根据上式由计算采样点的近似值。采用上述方法计算信号的频谱需要注意如下三个问题:word 文档 可自由复制编辑(1)频谱混叠。如果不满足采样定理的条件,频谱会出现混叠误差。对于频谱无限宽的信号,应考虑覆盖大部分主要频率分量的范围。(2)栅栏效应和频谱分辨率。使用 DFT计算频谱,得到的结果只是 N 个频谱样本值,样本值之间的频谱是未知的,像通过一个栅栏观察频谱,称为“栅栏效应”。频谱分辨率与记录长度成反比,要提高频谱分辨率,就要增加记录时间。(3)频谱泄露。对信号截断会把窗函数的频谱引入
5、信号频谱,造成频谱泄露。解决这个问题的主要办法是采用旁瓣小的窗函数,频谱泄露和窗函数均会引起误差。因此,要合理选取采样间隔和截取长度,必要时还需考虑加适当的窗。对于连续时间周期信号,我们在采用计算机进行计算时,也总是要进行截断,序列总是有限长的,仍然可以采用上述方法近似计算。4.可能用到的 MATLAB函数与代码 实验中 DFT运算可采用 MATLAB中提供的函数 fft来实现。DTFT可采用 MATLAB矩阵运算的方法进行计算,如下式所示:NNjnjnjnNnnnnjnjeeenxnxnxenxeX211.,)(21 四、实验内容 1、已知1,1,1,2)(nx,完成如下要求:(1)计算其
6、DTFT,并画出,区间的波形。(2)计算 4 点 DFT,并把结果显示在(1)所画的图形中。(3)对 x(n)补零,计算 64 点 DFT,并显示结果。(4)根据实验结果,分析是否可以由 DFT计算 DTFT,如果可以,如何实现。解:(1)计算其 DTFT,并画出,区间的波形。n=0:3;x=2-1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n*w);word 文档 可自由复制编辑 subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis tight;subplot(212);plot(w,angl
7、e(X)/pi);xlabel(Omega/pi);title(Phase);axis tight;(2)计算 4 点 DFT,并把结果显示在(1)所画的图形中。Xk=fft(x);subplot(211);hold on;stem(n,abs(Xk),filled);plot(w,abs(X);axis tight;xlabel(Omega/pi);title(Magnitude);subplot(212);hold on;plot(w,angle(X)/pi);stem(n,angle(Xk),filled);axis tight;-3-2-10123123/Magnitude-3-2-1
8、0123-0.500.5/Phaseword 文档 可自由复制编辑 xlabel(Omega/pi);title(Phase);运行结果如下:(3)对 x(n)补零,计算 64 点 DFT,并显示结果。x=2-1 1 1 zeros(1,60);n=0:63;Xk=fft(x);subplot(211);stem(n,abs(Xk),filled);subplot(212);stem(n,angle(Xk),filled);-3-2-101230123/Magnitude-3-2-10123-1-0.500.51/Phaseword 文档 可自由复制编辑 (4)根据实验结果,分析是否可以由 D
9、FT计算 DTFT,如果可以,如何实现。可以由 DFT 计算 DTFT,序列补零后,长度越长,DFT 点越多,其 DFT 越逼近 DTFT连续波形。所以令序列补零至足够长时其 DFT 序列的波形与 DTFT的波形在一定的分辨率下已经相同。2、考察序列 x(n)=cos(0.48 n)+cos(0.52 n)(1)0=n=10时,用 DFT估计 x(n)的频谱;将 x(n)补零加长到长度为 100点序列用 DFT估计 x(n)的频谱,要求画出相应波形。(2)0=n=100时,用 DFT估计 x(n)的频谱。并画出波形。(3)根据实验结果,分析怎样提高频谱分辨率 解:(1)0=n n=0:10;x
10、=cos(0.48*pi.*n)+cos(0.52*pi.*n);Xk=fft(x);subplot(211);01020304050607001234010203040506070-2-1012word 文档 可自由复制编辑 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);012345678910-2024680102030405060708090100-10-50510word 文档 可自由复制编辑(2)0=n n=0:100;x=co
11、s(0.48*pi.*n)+cos(0.52*pi.*n);Xk=fft(x);stem(Xk,filled);axis tight;(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做频谱分析,确定适合的参数,使得到的频谱的频率分辨率符合需要。1020
12、30405060708090100-30-20-100102030word 文档 可自由复制编辑 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);word 文档 可自由复制编辑 n=0:200;x=0.15*sin(0.1*pi.*n)+sin(0.2*pi.*n)-0.1*sin
13、(0.3*pi.*n);Xk=fft(x);stem(abs(Xk),filled);word 文档 可自由复制编辑 结果分析:上图为 x(t)信号截取过后的连续时间信号的傅里叶变换幅频特性曲线,截取周期为 100s(即为采样时间 M),根据(j)kXMc 的特点,知 1hz处的幅值为 0.15(j2*1)1007.52X,(j2*2)5X,(j2*3)50X与图像相符。4、利用 DFT近似分析连续时间信号 x(t)=e-0.1 u(t)的频谱(幅度值)。分析采用不同的采样间隔和截取长度进行计算的结果,并最终确定合适的参数。n=0:10;x=exp(-0.1.*n);Xk=fft(x);ste
14、m(Xk,filled);n=0:20;x=exp(-0.05.*n);Xk=fft(x);stem(Xk,filled);word 文档 可自由复制编辑 n=0:40;x=exp(-0.025.*n);Xk=fft(x);stem(Xk,filled);word 文档 可自由复制编辑 五、心得体会 通过本次实验,加深了对 DFT原理的理解。学会了应用 DFT分析信号的频谱。深刻理解到利用 DFT分析信号频谱的原理,能够分析实现过程中出现的现象及解决方法。word 文档 可自由复制编辑 实验二 利用 FFT计算线性卷积 一、实验目的 1.掌握利用 FFT计算线性卷积的原理及具体实现方法。2.加
15、深理解重叠相加法和重叠保留法。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 算法,要求 。采
16、用补零的办法使得 x(n)和 h(n)的长度均为 N。(2)计算 x(n)和 h(n)的 N 点 FFT。word 文档 可自由复制编辑(3)组成乘积 (3)利用 IFFT计算 Y(k)的 IDFT,得到线性卷积 y(n)3.分段卷积 我们考察单位取样响应为 h(n)的线性系统,输入为 x(n),输出为 y(n),则 如果 x(n)极长时,如果要等 x(n)全部集齐时再开始进行卷积,会使输出有较大延时;如果序列太长,需要大量存储单元。为此,我们把 x(n)分段,为别求出每段的卷积,合在一起得到最后的总输出。这称为分段卷积。分段卷积可以细分为重叠保留法和重叠相加法。重叠保留法:设 x(n)的长度
17、为 ,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计算圆周卷积 word 文档 可自由复制编辑(4)设每一段
18、卷积结果的前 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)。三、实验内容 假设要计算序列()()()
19、,0 x nu nu nLnL和()cos(0.2),0h nnnM的线性卷积,完成以下实验内容:设 L=M,根据线性卷积的表达式和快速卷积的原理,分别编程实现计算两个序列线性卷积的方法,比较当序列长度分别为 8,16,32,64,256,512,1024时,两种计算方法计算线性卷积所需时间。for i=1:7 L=input(L:);n=0:L;word 文档 可自由复制编辑 x=heaviside(n)-heaviside(n-L);h=cos(0.2*pi.*n);tic y=conv(x,h);toc end L:8 时间已过 0.000050 秒。L:16 时间已过 0.000037
20、 秒。L:32 时间已过 0.000056 秒。L:64 时间已过 0.000065 秒。L:256 时间已过 0.000101 秒。L:512 时间已过 0.000188 秒。L:1024 时间已过 0.000254 秒。for i=1:7 L=input(L:);n=0:L;x=heaviside(n)-heaviside(n-L);h=cos(0.2*pi.*n);tic Xk=fft(x,L+1);word 文档 可自由复制编辑 Hk=fft(h,L+1);Yk=Xk.*Hk;y=ifft(Yk);toc end L:8 时间已过 0.000048 秒。L:16 时间已过 0.0000
21、44 秒。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时两种算法所需时间。for i=1:2 L=input(L:);M=input(M:);n1=0:L;n2=0:M;x=heaviside(n1)-heaviside(n1-L);h=cos(0.2*pi.*n2);word 文档 可自由复制编辑 tic Xk
22、=fft(x,L+1);Hk=fft(h,L+1);Yk=Xk.*Hk;y=ifft(Yk);toc end L:2048 M:256 时间已过 0.001139 秒。L:4096 M:256 时间已过 0.001503 秒。for i=1:2 L=input(L:);M=input(M:);n1=0:L;n2=0:M;x=heaviside(n1)-heaviside(n1-L);h=cos(0.2*pi.*n2);tic y=conv(x,h);toc end L:2048 M:256 时间已过 0.000263 秒。L:4096 word 文档 可自由复制编辑 M:256 时间已过 0.
23、000368 秒。3.L=input(L:);M=input(M:);n2=0:M;k=L/M;h=cos(0.2*pi.*n2);Y=0;tic for i=1:k n1=(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;end toc L:2048 M:256 时间已过 0.025870 秒。L=input(L:);M=input(M:);n2=0:M;k=L/M;word 文档 可自由复制编辑 h=cos(0.2*pi.*n2);Y=0;ti
24、c for i=1:k n1=(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;end toc L:4096 M:256 时间已过 0.022300 秒。L=4096;n1=0:L;x=heaviside(n1)-heaviside(n1-L);n2=0:256;h=cos(0.2*pi.*n2);tic y1=conv(x,h);toc tic Xk=fft(x);Hk=fft(h);Yk=Xk.*Hk;word 文档 可自由复制编辑 y2=if
25、ft(Yk);toc 编程实现利用重叠相加法计算两个序列的线性卷积,考察 L=2048且 M=256时计算线性卷积的时间,与第二题结果进行比较。clear tic N=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);for i=Lenx:K*L-1 x(i+1)=0;end Y=zeros(K,N);YY=zeros(1,(K-1)*L+N);for k=0:K-1 xk
26、=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,:);end toc word 文档 可自由复制编辑 时间已过 0.028816 秒。编程实现利用重叠保留法计算两个序列的线性卷积,考察 L=2048且 M=256时计算线性卷积的时间,与第二题结果进行比较。clc clear tic N=512;m=0:256;h=cos(0.2*pi*m);n=0:2048;x=heaviside(n)-heaviside(n-2048);Lenx=length(x);M=
27、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);for k=0:K-1 xk=fft(x1(k*L+1:k*L+N);Y(k+1,:)=(ifft(xk.*h);end Z=reshape(Y(:,M:N),1,);toc 时间已过 0.044424 秒。word 文档 可自由复制编辑 四、实验心得 通过本次实验,掌握了利用 FFT计算线性卷积的原理及具体实现方法,加深了理解重叠相加法和重叠保留法。word 文档 可自
28、由复制编辑 实验 3 IIR 数字滤波器设计 一、实验目的 掌握利用脉冲响应不变法设计 IIR数字滤波器的原理及具体方法。加深理解数字滤波器和模拟滤波器之间的技术指标转化。掌握脉冲响应不变法设计 IIR数字滤波器的优缺点及适用范围。二、实验内容 设采样频率为sf 10kHz,设计数字低通滤波器,满足如下指标 通带截止频率:pf 1kHz,通带波动:pR 1dB:阻带截止频率:sf 1.5kHz,阻带衰减:sR 15dB:要求分别采用巴特沃斯、切比雪夫型、切比雪夫型和椭圆模拟原型滤波器,并分别结合脉冲响应不变法和双线性变换法进行设计。结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足
29、给定指标要求,分析脉冲响应不变法和双线性变换法设计 IIR数字滤波器的优缺点及适用范围。1 巴特沃斯 脉冲响应不变法: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(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)word 文档 可自由复制编辑 Omegac1=
30、Wp/(10(Rp/10)-1)(1/(2*N1)Omegac2=wp/(10(Rp/10)-1)(1/(2*N1)z,p,k=buttap(N1)b1=k*Omegac1N1 a1=poly(p*Omegac1)H,w=freqs(b1,a1)subplot(221)plot(w/pi,abs(H)grid on subplot(223)plot(w/pi,angle(H)grid on b2,a2=butter(N2,Omegac2,S)H,w=freqs(b2,a2)subplot(222)plot(w/pi,abs(H)axis tight grid on subplot(224)pl
31、ot(w/pi,angle(H)grid on word 文档 可自由复制编辑 双线性法: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(2
32、21);word 文档 可自由复制编辑 plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);plot(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subplot(224);plot(w/pi,grd);grid on;2 切比雪夫 I 型 脉冲响应不变法:wp=0.2*pi;word 文档 可自由复制编辑 ws=0.3*pi;Rp=1;As=15;F=10000;T=1/F;e=(10(Rp/10)-1)(1
33、/2);A=10(As/20);wp1=wp/T;ws1=ws/T;N=ceil(acosh(A2-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);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);plot(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subpl
34、ot(224);plot(w/pi,grd);grid on;word 文档 可自由复制编辑 双线性变换法: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(A2-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);word 文档 可自由复制编辑 subplot(221);pl
35、ot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);plot(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subplot(224);plot(w/pi,grd);grid on;3 切比雪夫 II 型 脉冲响应不变法:wp=0.2*pi;word 文档 可自由复制编辑 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=w
36、p/T;ws1=ws/T;N=ceil(acosh(A2-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);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);plot(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subplot(224);plot(w/pi,grd
37、);grid on;word 文档 可自由复制编辑 双线性变换法: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(A2-1)(1/2)/e)/(acosh(ws1/wp1)wc=ws/pi b,a=cheby2(N,As,wc);word 文档 可自由复制编辑 w=0:500*pi/500;H,w=freqz(b,a);subplot(221);plot(w/pi,abs(H);grid o
38、n;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);plot(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subplot(224);plot(w/pi,grd);grid on;4 椭圆原型模拟滤波器 word 文档 可自由复制编辑 脉冲响应不变法: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/ws
39、1;k1=e/(A2-1)(1/2);K1=ellipke(k);K2=ellipke(1-k12)(1/2)K3=ellipke(k1)K4=ellipke(1-k2)(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);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;subplot(223);word 文档 可自由复制编辑 plo
40、t(w/pi,angle(H)/pi);grid on;grd=grpdelay(b,a,w);subplot(224);plot(w/pi,grd);grid on;双线性法: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);word 文档 可自由复制编辑 A=10(As/20);k=wp1/ws1;k1=e/(A2-1)(1/2);K1=ellipke(k);K2=ellipke(1-k12)(1/2)K3=ellipke(k1)K4=e
41、llipke(1-k2)(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);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,ang
42、le(H)/pi);grid on;xlabel(Omega(pi);ylabel(Phase of H(ejOmega)(pi);grd=grpdelay(b,a,w);subplot(224);word 文档 可自由复制编辑 plot(w/pi,grd);grid on;xlabel(Omega(pi);ylabel(Group delay);00.5100.51()|H(ej)|00.20.40.60.81-20-15-10-50()|H(ej)|(dB)00.51-1-0.500.51()Phase of|H(ej)|()00.51-1000100200()Group delay 三、
43、结果分析:(1)对于设计要求的符合情况 根据以 dB 为单位的幅频特性图可以看出,能满足设计要求的滤波器有:脉冲响应不变法的切比雪夫型,以及双线性变换法的全部四种滤波器。而脉冲响应不变法的巴特沃斯、切比雪夫型和椭圆数字滤波器不满足设计要求。(2)脉冲响应不变法的优缺点及适用范围 优点:幅频特性与模拟滤波器线性对应。脉冲响应不变法的一个重要特点是频率坐标的变换是线性的,=,与 是线性关系。因此如果模拟滤波的频响带限于折叠频率以内的话,通过变换后滤波器的频响可不失真地反映原响应与频率的关系。稳定性好。如果 Ha(s)是稳定的,即其极点在 S 左半平面,映射到H(Z)也是稳定的。缺点:设计结果不固定
44、。在设计时,往往不能直接根据要求的参数得到合适的word 文档 可自由复制编辑 滤波器,需要设计者多次修改设计参数。可能产生混叠失真。因频谱周期延拓效应,脉冲响应不变法只能用于带限的频响特性,如衰减特性很好的低通或带通。而高频衰减越大,频响的混淆效应越小,至于高通和带限滤波器,由于它们在高频部分不衰减,因此将完全混淆在低频响应中。所以用脉冲响应不变法实现高通和带限滤波器时,应增加一保护滤波器,滤掉高于折叠频率以上的频带,然后再用脉冲响应不变法转换为数字滤波器,这会增加设计的复杂性和滤波器的阶数。适用范围:根据幅频特性与模拟滤波器相对应的特点,在要求时域脉冲响应能模仿模拟滤波器的场合,一般使用脉
45、冲响应不变法。(3)双线性变换法的优缺点及适用范围 优点:消除了脉冲响应不变法所固有的混叠误差。靠频率的严重非线性关系得到 S 平面与 Z 平面的单值一一对应关系,整个 j 轴单值对应于单位圆一周,其中 和 为非线性关系。设计结果具有一般性。双线性变换法根据设计要求和计算流程计算得到的滤波器,一般能够满足原来的要求,得到的频率特性较好。双线性变换比脉冲响应法的设计计算更直接和简单。由于 s 与 z 之间的简单代数关系,所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。缺点:频率响应有畸变。与 的非线性关系,导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变,(使数字滤波器与模
46、拟滤波器在响应与频率的对应关系上发生畸变)。例如,一个模拟微分器,它的幅度与频率是线性关系,但通过双线性变换后,就不可能得到数字微分器。一个线性相位的模拟滤波器经双线性变换后,滤波器就不再有线性相位特性。适用范围:双线性变换目前仍是使用得最普遍、最有成效的一种设计工具。这是因为大多数滤波器都具有分段常数的频响特性,如低通、高通、带通和带阻等,它们在通带内要求逼近一个衰减为零的常数特性,在阻带部分要求逼近一个衰减为的常数特性,这种特性的滤波器通过双线性变换后,虽然频率发生了非线性变化,但其幅频特性仍保持分段常数的特性。四、实验心得 通过本次实验,学会了掌握利用脉冲响应不变法设计 IIR数字滤波器
47、的原理及具体方法。加深了理解数字滤波器和模拟滤波器之间的技术指标转化。掌握到脉冲响应不变法设计 IIR数字滤波器的优缺点及适用范围 word 文档 可自由复制编辑 实验 4 FIR 数字滤波器设计 一、实验目的 掌握窗函数法设计 FIR数字滤波器的原理及具体方法。二、实验设备与环境 计算机、MATLAB软件环境。三、实验内容 1.设计一个数字低通滤波器 FIR滤波器,其技术指标如下 p0.2,pR0.25dB st0.3,sA50dB 分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。(1)矩形窗 程序代
48、码如下: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%根据性能指标确定滤波器长度 N n=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);%画出时域和频域图形,以确定设计出来的滤波器是否符合给定的指标
49、 word 文档 可自由复制编辑 stem(n,hd,filled);Hr,wl=zerophase(h);subplot(222);plot(wl/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;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);grid on;%set(gca,xt
50、ick,(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(0 1-2 0);运行结果如下:N=19 word 文档 可自由复制编辑 05101520-0.100.10.20.300.20.40.60.800.51/H()05101500.10.2nh(n)00.51-100-500/dB 05101520-0.100.10.20.300.20.40.60.800.51/H()05101500.10.2nh(n)00.20.40