数字信号处理实验报告

1、目 录实验1 离散时间信号的频域分析-2实验2 FFT算法与应用-7实验3 IIR数字滤波器的设计-12实验4 FIR数字滤波器的设计-17实验1 离散时间信号的频域分析一实验目的信号的频域分析是信号处理中一种有效的工具。在离散信号的时域分析中，通常将信号表示成单位采样序列（n）的线性组合，而在频域中，将信号表示成复变量e或 e的线性组合。通过这样的表示，可以将时域的离散序列映射到频域以便于进一步的处理。在本实验中，将学习利用MATLAB计算离散时间信号的DTFT和DFT,并加深对其相互关系的理解。2、 实验原理(1)DTFT和DFT的定义及其相互关系。序列x(n)DTFT定义为= 它是关于自

2、变量的复函数，且是以2为周期的连续函数。可以表示为,其中，和分别是实部和虚部；还可以表示为 =，其中， 和分别是的幅度函数和相位函数；它们都是的实函数，也是以2为周期的周期函数。 序列的N点DFT定义为，是周期为N的序列。与的关系：是对）在一个周期中的谱的等间隔N点采样，即 ，而可以通过对内插获得，即 (2)使用到的MATLAB命令有基于DTFT离散时间信号分析函数以及求解序列的DFT函数。 1）基于 DTFT离散时间信号分析函数有：freqz,real,imag,abs,angle。函数freqz可以用来计算一个以e的有理分式形式给出的序列的DTFT值。freqz的形式多样，常见的有H=fr

3、eqz(num,den,w)，其中num表示序列有理分式DTFT的分子多项式系数，den表示分母多项式系数（均按z的降幂排列），矢量w表示在02中给定的一系列频率点集合。freqz函数的其他形式参见帮助文件。在求出DTFT值后，可以使用函数real,imag,abs和angle分别求出并绘出起实部、虚部和相位谱。例如）=利用函数freqz计算出，然后利用函数abs和angle分别求出幅频特性与相位特性最后利用plot命令绘出曲线。2） 求解序列DFT的函数有：fft,ifft。函数fft（x）可以计算R点序列的R点DFT值；而fft（x，N）则计算R点序列的N点DFT，若R>N，则直接截

4、取R点DFT的前N点，若R<N,则x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算R点的谱序列的R点IDFT值；而ifft(X,N)同fft(x,N)的情况。三、实验设备计算机、MATLAB软件四、实验内容（1)编程计算并画出下面DTFT的实部，虚部、幅度和相位谱。程序如下：num=0.1313 -0.1553 0.1313 0.0518;den=1 1.2828 1.0388 0.3418;w=0:0.001:2*pi;H=freqz(num,den,w);figuresubplot(221)plot(w,real(H);title('实部')gri

5、d on;subplot(222)plot(w,imag(H)title('虚部')gridsubplot(223)plot(w,abs(H)title('幅度')gridsubplot(224)plot(w,abs(H)title('相位谱')grid on;波形如下：图1-1 DTFT的实部，虚部、幅度和相位谱(2)计算32点序列x(n)=cosn,0n31的32点和64点DFT，分别绘出幅度谱图形，并绘出该序列的DTFT图形。程序如下：n1=0:15;n2=0:31;x1=cos(13*pi*n1)/32)x2=cos(13*pi*n2)/

6、32)X1=fft(x1);X2=fft(x2);subplot(211)plot(n1,abs(X1)title('32点DFT')grid on;subplot(212)plot(n2,abs(X2)title('64点DFT')grid on; 波形如下：图1-2 序列DFT幅度谱图形和DTFT图形DTFT程序如下：A=1;n=0:31;x=cos(13*pi*n)/32);B=x;w=0:0.01:2*pi*2;H=freqz(B,A,w);magH=abs(H);subplot(2,1,1);plot(w,magH);grid;ylable('

7、Magnitude');subplot(2,1,2);plot(w,phaH);grid;xlable('w');ylable('phase');grid on;波形如下：图1-3 DTFT图形实验2 FFT算法与应用一、实验目的在理论学习的基础上，通过本次实验。加深对快速傅里叶变换的理解，熟悉FFT算法。熟悉应用FFT对典型信号进行频谱分析的方法。了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。二实验原理在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。无限长序列往往也可以用有限长序列来逼近。对于有限长序列可

8、以使用离散傅里叶变换（DFT），这一变换可以很好地反映序列的频域特性，并且容易利用快速算法在计算机的上的实现。当序列的长度是N时，定义离散傅里叶变换为 =式中，的反变换定义为 根据上式令，则有所以，X（k）是在z变换在单位圆上的等间隔采样，或者说是序列傅里叶变换的等间隔采样。时域采样在满足乃奎斯特定理时，就是不会发生频谱混叠；同样，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列的混叠。DFT是对序列傅里叶变换的等间距采样，因此可以利用于序列的频谱分析。在运用于DFT进行频谱分析的时候可能有三种误差，分析如下：（）混叠现象序列的频谱是采样信号频谱的周期延拓，周期是，因此当采样速率

9、不满足乃奎斯特定理时，即采样频率小于两倍的信号频率时，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实地反映原信号的频谱。（）泄露现象实际中的信号序列往往很长，甚至是无限长序列，为了方便，往往用截断的序列来近似它们。这样可以使用较短的DFT来对信号进行频谱分析。这种截断等价于给原信号序列乘以一个矩形窗函数，而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。（3）栅栏响应因为DFT是对单位圆上z变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度来看，用DFT来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上

10、看到真实的频谱。减小栅栏效应的一个方法是在源程序的末端补一些零值，从而变动DFT的点数。函数fft（x）可以计算R带你序列的R点的DFT值；而fft（x，N）则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R<N，则x先进行补零扩展为N点序列再求N点DFT。函数ifft（X）可以计算R点的谱序列的R点IDFT值；而ifft（X,N）同fft（x，N）的情况。函数conv（x1，x2）可以计算两序列的线性卷积。三 实验设备计算机，MATLAB软件。四 实验内容编制信号产生子程序及本实验的频谱分析主程序。实验中需要用到的基本信号包括：（1） 三角波序列： x(n)

11、=n+13,0 13-n ,4 0,其它程序如下：n2=0:7;x2=13,14,15,16,9,8,7,6;subplot(221);stem(n2,x2);title('x2序列');grid on;k1=0:7;y21=fft(x2,8);magy21=abs(y21);subplot(222);stem(k1,magy21);title('x2的8点FFT');grid on;k2=0:15;y22=fft(x2,16);magy22=abs(y22);subplot(224);stem(k2,magy22);title('x2的16点FFT&#

12、39;);grid on;波形如下：图2-1 三角波序列以及它的8点和16点FFT（2） 利用DFT计算下面两序列的线性卷积x(n)=3,-4,6,0,6,-4,h(n)=1,-4,3,4程序如下：clear all; close all; clc; N=9;g=3 -4 6 0 6 -4; h=1 -4 3 4; x=conv(g,h) gf=fft(g,N); hf=fft(h,N); xx=gf.*hf X=ifft(xx,N); grid on;波形如下： 图2-2 用DFT计算的序列x(n)与h(n)线性卷积图（3） 已知某序列x（n）在单位圆上的N=64等分样点的z变换为X（Z）=

13、X（k）=，k=0,1,2,63。试分别用N点IFFT程序和N点FFT程序计算,绘出及。程序如下：k=0:63;N=64;X=1./(1-0.13*exp(-j*2*pi*k/N);x=ifft(X,64)stem(abs(x);xlabel('n');ylabel('x(n)');title('25FF.单位圆上等分样点的Z变换图')grid on;波形如下： 图2-3 单位圆等分样点Z变换图实验3 IIR数字滤波器的设计一 实验目的 从理论上讲，任何的线性时不变（LTI）离散时间系统都可以看做一个数字滤波器，因此设计数字滤波器实际就是设计离散

14、时间系统。数字滤波器包括IIR（无限冲击响应）和FIR（有限冲击响应）型，在设计时通常采用不同的方法。 本实验通过使用MATLAB函数对数字滤波器进行设计和实现，要求掌握IIR数字巴特沃斯滤波器,数字切比雪夫滤波器的设计原理，设计方法和设计步骤；能根据给定的滤波器指标进行滤波器设计；同时也加深学生对数字滤波器的常用指标和设计过程的理解。二 实验原理在IIR滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器Ha(s),然后根据一定的准则将比模拟原型滤波器转换为数字滤波器Hd(z),即我们需要设计的数字滤波器。IIR滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定

15、主要是在模拟原型滤波器设计中进行的。IIR数字滤波器的设计方法如下：（1）冲击响应不变法。对满足设计要求的模拟原型滤波器Ha(s)进行部分分式展开为 ;基于=，可得(2)双线性变换法。对设计要求中给出的边界频率进行预畸处理，然后用得到的频率进行模拟滤波器设计，得到模拟原型滤波器Ha(s);用双线性变换法求出数字滤波器。一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。冲击响应不变法一个重要特点是频率坐标的变换是线性的，因此如果模拟滤波器的频率响应带限于折叠频率的话，则通过变换后滤波器的频率响应可不失真地反映原响应与频率的关系。与冲激响应不变法比较，双线性变换的主要优点

16、是靠频率的非线性关系得到s平面与z平面的单值一一对应关系，整个值对应于单位圆一周。所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。MATLAB提供了一组标准的数字滤波器设计函数，大大简化了滤波器的设计过程。（1） butter。功能：Butterworth模拟/数字滤波器设计如图；格式：b,a=butter(n,wn,ftype,s),b,a=butter(n,wn,ftype)。（2） Cheby1、cheby2。功能：chebyshev、chebyshev型模拟/数字滤波器设计如图所示。格式：b,a=cheby1(n,Rp,wn,ftype,),b,a=cheby2(n,Rs

17、,wn,ftype)IIR滤波器设计函数还包括：buttord,chebwin,cheb1ord,cheb2ord,ellip,ellipord等。还可以用N,Wn=buttord(Wp,Ws,Rp,Rs)的MATLAB命名估算一个Butterworth滤波器的阶数。三实验设备计算机，MATLAB软件。四实验内容利用MATLAB编程方法或利用MATLAB中的fdatool工具设计不同功能的IIR数字滤波器。（1）基于chebyshev型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求参数为通带截止频率=0.213；通带最大衰减=1dB;阻带截止频率=0.413；阻带最小衰减=35dB。程

18、序如下：clear;wp=0.213*pi;ws=0.413*pi;Ap=1;As=35;T=1; Rip=10(-Ap/20);Atn=10(-As/20);OmegaP=wp*T;OmegaS=ws*T;n,Wn=cheb1ord(OmegaP,OmegaS,1,35,'s'); b,a=cheby1(n,1,Wn,'low','s'); freqs(b,a); bz,az=impinvar(b,a,T);H,W=freqz(bz,az,512,T); ma=20*log10(abs(H),pha=20*log10(unwrap(angle(

19、H),hi=impz(bz,az);ni=step(bz,az);subplot(2,2,1),plot(W,ma);title('幅频特性(dB)');xlabel('w(/pi)');ylabel('dB');subplot(2,2,2),plot(W,pha);title('相频特性');xlabel('w(/pi)');ylabel('pha(/pi)');subplot(2,2,3),plot(hi);title('单位冲击响应');xlabel('n');

20、ylabel('h(n)');subplot(2,2,4),plot(ni); title('单位阶跃响应'); xlabel('n');ylabel('h(n)');bz,az; 波形如下：图3-1 基于chebyshevI模型的数字滤波器(2)基于Butterworth型模拟滤波器原型使用双线性变换方法设计数字滤波器的，要求参数为通带截止频率=0.24；通带最大衰减=1.3dB;阻带截止频率=0.13；阻带最小衰减=40dB。程序如下：clear;wp=0.4*pi; ws=0.13*pi; Ap=1.3; As=40; fs

21、=1; ts=1/fs; wp2=2*fs*tan(wp/2*ts); ws2=2*fs*tan(ws/2*ts); n,wn=buttord(wp2,ws2,Ap,As,'s') z,p,k=buttap(n) bap,aap=zp2tf(z,p,k) b,a=lp2lp(bap,aap,wn) bz,az=bilinear(b,a,fs); h,w=freqz(bz,az); subplot(2,1,1);plot(w/pi,abs(h);grid onxlabel('频率');ylabel('幅度');subplot(2,1,2);plot

22、(w/pi,20*log10(abs(h);grid onxlabel('频率');ylabel('幅度'); 波形如下:图3-2 基于Butterworth模型的数字滤波器 实验4 FIR数字滤波器的设计一实验目的掌握用窗函数法设计FIR数字滤波器的原理与其设计步骤；熟悉线性相位数字滤波器的特性。学习编写数字滤波器的设计程序的方法，并能进行正确编程；根据给定的滤波器指标，给出设计步骤。二实验原理如果系统的冲击响应为已知，则系统的输入-输出关系为: 对于低通滤波器，只要设计出低通滤波器的冲击响应函数，就可以由式得到系统的输出了。假设所希望的数字滤波器的频率响应为

23、H(e)，它是频域的周期函数，周期为2，那么它与H(e)相对应的傅里叶系数为: 但是将h（n）作为滤波器脉冲响应有两个问题：一是它是无限长的，与FIR数字滤脉冲响应有限长这一前提不一致；二是它是非因果的，。对此，要采取两项措施：一是将截短；二是将其往右移。由此得到的实际频域响应，与理想频域响应相近，但不完全一致。理论证明上述现象是对进行简单截短处理的必然结果，一般称为吉布斯现象，为尽可能的减少吉布斯现象，应对进行加窗截取，即以作为FIR滤波器的系数。其中N为窗口的长度。窗口函数的形状和窗口长度N决定了窗函数法设计出的FIR滤波器的性能。设计时，要根据阻带的最小衰减和过渡带宽来选择恰当的窗函数类

24、型和窗长度N。常用窗函数有矩形窗、海明窗、和布莱克曼窗等。窗函数设计FIR滤波器步骤如下：（1）给定理想频率响应的幅频特性和相频特性；（2）求理想单位脉冲响应，在实际计算中，可对采样，并对其求IDFT的，用代替；（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗长度N；（4）根据,求所需设计的FIR滤波器单位脉冲响应；（5）求，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。三实验设备计算机，MATLAB软件。四实验内容(1)用海明窗设计一个13阶的FIR的低通滤波器，wc=0.13程序如下：clear all;%海明窗N=13;wc=0.13*pi;h=fir1(N-1,wc/pi,'b