《数字信号处理》实验报告

1、物理与信息工程学院数字信号处理实验报告年级：2011级班级：信通4班 姓名：朱明贵学号：111100443老师：李娟福州大学2013 年 11 月实验一 快速傅里叶变换(FFT)及其应用一、实验目的1. 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。2. 熟悉应用FFT对典型信号进行频谱分析的方法。3. 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。4. 熟悉应用FFT实现两个序列的线性卷积和相关的方法。二、实验类型演示型三、实验仪器装有MATLAB语言的计算机四、实验原理 在各种信号序列中，有限长序列信号处理占有

2、很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：反变换为：     有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。    FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT是以2为基数的，其长度 。它的效率高，程

3、序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。  (一)在运用DFT进行频谱分析的过程中可能的产生三种误差1 混叠    序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信

4、号进行滤波。2泄漏    实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。3栅栏效应    DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样

5、，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。    减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。    (二)用FFT计算线性卷积    用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于

6、线性卷积的充要条件是FFT的长度NN1N2-1，对于长度不足N的两个序列，分别将他们补零延长到N。    当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。有两种方法：1 重叠相加法将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。2重叠保留法这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。 五、实验内容和要求    1、一个连续信号含两个频率分量，经采样得已知，分别为1/1

7、6和1/64，观察其频谱；当时，不变，其结果有何不同，为什么？代码：N=16;n=0:N-1;Df=1/16;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(245);stem(n,abs(Xk);xlabel('n');ylabel('X(k)');title('N=16,Df=1/16,频谱图');subplot(241);stem(n,abs(Xn);xlabel('n');ylabel('X(n)');title('N

8、=16,Df=1/16,时序图');Df=1/64;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(246);stem(n,abs(Xk);xlabel('n');ylabel('X(k)');title('N=16,Df=1/64,频谱图');subplot(242);stem(n,abs(Xn);xlabel('n');ylabel('X(n)');title('N=16,Df=1/64,时序图');N=1

9、28;n=0:N-1;Df=1/16;Xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(247);stem(n,abs(Xk);xlabel('n');ylabel('X(k)');title('N=128,Df=1/16,频谱图');subplot(243);stem(n,abs(Xn);xlabel('n');ylabel('X(n)');title('N=128,Df=1/16,时序图');Df=1/64;Xn=s

10、in(2*pi*0.125*n)+cos(2*pi*(0.125+Df)*n);Xk=fft(Xn,N);subplot(248);stem(n,abs(Xk);xlabel('n');ylabel('X(k)');title('N=128,Df=1/64,频谱图');subplot(244);stem(n,abs(Xn);xlabel('n');ylabel('X(n)');title('N=128,Df=1/64,时序图');频谱图：    2、用FFT分别实现的和

11、的线性卷积和10点、20点圆周卷积，记录其波形，并说明他们之间的关系。代码：function y = circonv( x1,x2,N )%UNTITLED3 Summary of this function goes here% Detailed explanation goes here if length(x1)>N error('N should higher than or equal to the length of x1') end if length(x2)>N error('') end x1=x1,zeros(1,N-length

12、(x1); x2=x2,zeros(1,N-length(x2); m=0:1:N-1; x2=x2(mod(-m,N)+1); H=zeros(N,N); for n=1:1:N H(n,:)=cirshift(x2,n-1,N); endy=x1*H'function y = cirshift( x,m,N )%UNTITLED4 Summary of this function goes here% Detailed explanation goes here if length(x)>N error('') end x=x zeros(1,N-length

13、(x); n=0:1:N-1; n=mod(n-m,N); y=x(n+1);n=0:11;x=0.8.n;h=1,1,1,1,1,1;N=17;x=x,zeros(1,5);h=h,zeros(1,11);X=fft(x);H=fft(h);Y=X.*H;y=ifft(Y);subplot(321);stem(x);xlabel('n');ylabel('x');title('序列x(n)');subplot(323);stem(h);xlabel('n');ylabel('h');title('序列h(

14、n)');subplot(325);stem(y);xlabel('n');ylabel('y');title('序列y(n),线性卷积');n=0:11;N=20;x=0.8.n;h=1,1,1,1,1,1;y=circonv(x,h,N);x=x,zeros(1,5);h=h,zeros(1,11);subplot(322);stem(x);xlabel('n');ylabel('x');title('序列x(n)');subplot(324);stem(h);xlabel('n

15、');ylabel('h');title('序列h(n)');subplot(326);stem(y);xlabel('n');ylabel('y');title('序列y(n)，圆周卷积');频谱图：说明他们之间的关系？周期卷积是线性卷积的周期延拓。 =+å六、思考题  1.实验中的信号序列和，在单位圆上的Z变换频谱和会相同吗？如果不同，说出哪一个低频分量更多一些，为什么？答：不同；其中，Xd(n)的低频分量较多，由图形可以看出，在低频处，Xd(n)的取值较多，呈递减趋势。2

16、. 对一个有限长序列进行DFT等价于将该序列周期延拓后进行DFS展开，因为DFS也只是取其中一个周期来计算，所以FFT在一定条件下也可以用以分析周期信号序列。如果实正弦信号用16点FFT来做DFS运算，得到的频谱是信号本身的真实谱吗？为什么？答：可以把有限长非周期序列假设为一无限长周期序列的一个主直周期，即对有限长非周期序列进行周期延拓，延拓后的序列完全可以采用DFS进行处理，即采用复指数。实验二 IIR数字滤波器的设计一、实验目的 1. 掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器

17、的MATLAB编程。 2. 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。 3. 熟悉Butterworth滤波器、Chebyshev滤波器和椭圆滤波器的频率特性。 二、实验类型设计型三、实验仪器装有MATLAB语言的计算机四、实验原理1 脉冲响应不变法用数字滤波器的单位脉冲响应序列模仿模拟滤波器的冲激响应,让正好等于的采样值，即，其中为采样间隔，如果以及分别表示的拉式变换及的Z变换，则2双线性变换法    S平面与z平面之间满足以下映射关系：s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面

18、的单位圆内。双线性变换不存在混叠问题。    双线性变换是一种非线性变换 ，这种非线性引起的幅频特性畸变可通过预畸而得到校正。   以低通数字滤波器为例，将设计步骤归纳如下：1. 确定数字滤波器的性能指标：通带临界频率、阻带临界频率、通带波动、阻带内的最小衰减、采样周期、采样频率 ；2.  确定相应的数字角频率，； 3.  计算经过预畸的相应模拟低通原型的频率，；4. 根据c和r计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数； 5. 用上面的双线性变换公式代入，求出所设计的传递函数； 6.  分

19、析滤波器特性，检查其指标是否满足要求。 五、实验内容和要求1用脉冲响应不变法设计一个巴特沃斯数字低通滤波器,要求在0-0.2内衰耗不大于3dB,在0.6 内衰耗不小于60dB，采样频率Fs=500 Hz. 代码：Omegap=0.2*pi*500;Omegas=0.6*pi*500;rp=3;as=60;fs=500;wp=Omegap/fs;ws=Omegas/fs;n,Omegac=buttord(Omegap,Omegas,rp,as,'s');b,a=butter(n,Omegac,'s');bz,az=impinvar(b,a,fs);w0=wp,ws

20、;hx=freqz(bz,az,w0);h,w=freqz(bz,az);plot(w*fs/(2*pi),abs(h);gridxlabel('频率hz')ylabel('频率响应幅度')频谱图：2 分别用脉冲响应不变法和双线性变换法分别设计一个巴特沃斯滤波器数字低通滤波器,使其特性逼近一个巴特沃斯模拟滤波器的性能指标如下：通带截止频率为2×2000（rad/s）,阻带截止频率为2×3000（rad/s），通带衰耗不大于3dB,阻带衰耗不小于15dB，采样频率Fs=100000 Hz，观察记录所设计数字滤波器的幅频特性曲线，记录带宽和衰减量

21、，检查是否满足要求。比较这两种方法的优缺点。代码（1）：omegap=2*pi*2000;omegas=2*pi*3000;rp=3;as=15;fs=10000;wp=omegap/fs;ws=omegas/fs;n,omegac=buttord(omegap,omegas,rp,as,'s');b,a=butter(n,omegac,'s');bz,az=impinvar(b,a,fs);w0=wp,ws;hx=freqz(bz,az,w0);h,w=freqz(bz,az);plot(w*fs/(2*pi),abs(h);gridxlabel('频

22、率hz')ylabel('频率响应幅度')屏幕图（1）：代码（2）：omegap=2*pi*2000;omegas=2*pi*3000;rp=3;as=15;fs=10000;wp=omegap/fs;ws=omegas/fs;omegap1=2*fs*tan(wp/2);omegas1=2*fs*tan(ws/2);n,omegac=buttord(omegap1,omegas1,rp,as,'s');b,a=butter(n,omegac,'s');bz,az=bilinear(b,a,fs);w0=wp,ws;hx=freqz(bz

23、,az,w0);h,w=freqz(bz,az);plot(w*fs/(2*pi),abs(h);gridxlabel('频率hz')ylabel('频率响应幅度')频谱图（2）：比较优缺点：脉冲响应不变法的优点：1，模拟频率到数字频率的转换时线性的。2，数字滤波器单位脉冲响应的数字表示近似原型的模拟滤波器单位脉冲响应，因此时域特性逼近好 缺点： 会产生频谱混叠现象，只适合带限滤波器双线性变换法优点： 克服多值映射得关系，可以消除频率的混叠 缺点： 是非线性的，在高频处有较大的失真。3 编写滤波器仿真程序，完成对实际采集的心电图信号序列（具体数据见下面）的总响应

24、序列，可直接调用MATLAB filter函数实现仿真。附：人体心电图采样信号在测量过程中往往受到工业高频干扰，所以，必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面的序列就是一个实际心电图信号采样序列样本，其中存在高频干扰，实验时，将其作为输入信号，滤除其中的干扰成分。   代码：clear all;wp=0.2*pi;ws=0.3*pi;rp=1;rs=15;n,wn=buttord(wp/pi,ws/pi,rp,rs);b,a=butter(n,wn);N=0.5/(0.02);figure(1);freqz(b,a,N);grid;xlabel

25、('频率hz');ylabel('频率响应幅度');xn=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;figure(2);subplot(2,1,1);stem(xn,'.');title('心电图信号采样序列');yn=filter(b,a,xn);subplot(2

26、,1,2)stem(yn,'.');title('滤波后的心电图信号');频谱图：六、思考题1. 双线性变换法中和之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。2.  能否利用公式 完成脉冲响应不变法的数字滤波器设计？为什么？答：IIR数字

27、滤波器的设计实际上是求解滤波器的系数和，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。但是它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。实验三 FIR数字滤波器的设计 一、实验目的1. 掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的MATLAB编程。 2. 熟悉线性相位FIR滤波器的幅频特性和相频特性。 3. 了解各种不同窗函数对滤波器性能的影响。 二、实验类型设计型三、实验仪器装有MATLAB语言的计算机四、实验原理&#

28、160;   线性相位实系数FIR滤波器按其N值奇偶和h(n)的奇偶对称性分为四种：    1h(n)为偶对称，N为奇数H(ej)的幅值关于=0，2成偶对称。     2h(n)为偶对称，N为偶数H(ej)的幅值关于=成奇对称，不适合作高通。    3h(n)为奇对称，N为奇数H(ej)的幅值关于=0，2成奇对称，不适合作高通和低通。    4h(n)为奇对称，N为偶数H(ej) =0、20，不适合作低通。FIR滤波器的常用设计方法：(一) 窗口法

29、60;   窗函数法设计线性相位FIR滤波器步骤 1. 确定数字滤波器的性能要求：临界频率、，滤波器单位脉冲响应长度N。 2. 根据性能要求，合理选择单位脉冲响应的奇偶对称性，从而确定理想频率响应的幅频特性和相频特性。 3. 求理想单位脉冲响应，在实际计算中，可对按M(M远大于N)点等距离采样，并对其求IDFT得，用代替。 4. 选择适当的窗函数，根据求所需设计的FIR滤波器单位脉冲响应。 5. 求，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。  窗函数的傅氏变换的主瓣决定了的过渡带宽。的旁瓣大小和多少决定了在通

30、带和阻带范围内波动的幅度。常用的几种窗函数有：(1) 矩形窗 w(n)=RN(n)； (2) Hanning窗 ； (3) Hamming窗 ； (4) Blackmen窗 ； (5) Kaiser窗 。 式中Io(x)为零阶贝塞尔函数。（二）频率采样法频率采样法是从频域出发，将给定的理想频率响应加以等间隔采样然后以此作为实际FIR数字滤波器的频率特性的采样值，即令由通过IDFT可得有限长序列将上式代入到Z变换中去可得其中()是内插函数五、实验内容及步骤1. N=50，编程并画出矩形窗、汉宁窗、海明窗和布莱克曼窗的时域波形和归一化的幅度谱，并比较各自的主要特点。代码：wvtool(boxcar

31、(50),hanning(50),hamming(50),blackman(50);grid;频谱图：并比较各自的主要特点?矩形窗优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象；汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降；海明窗的第一旁瓣衰减为一42dB海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB/（10oct），这比汉宁窗衰减速度慢；布莱克曼窗的幅度函数主要由五部分组成，他们的位移都不同，其幅度也是不同的WRg（w）使旁瓣再进一步抵消。

32、旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。设计程序时用backman函数调用。 2. N=15，带通滤波器的两个通带边界分别是，。用汉宁（Hanning）窗设计此线性相位带通滤波器，观察它的实际3dB和20dB带宽。N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N变化的影响。N=15汉宁代码：clc;clear alln=15;w1=0.3;w2=0.5;wn=w1,w2;b2=fir1(n,wn,hanning(n+1);freqz(b2,1);title('汉宁窗，N=15');N=15汉宁频谱：N=45汉宁代码：clc;clear alln=45;w1=0.3;w2=0.5;wn=w1,w2;b2=fir1(n,wn,hanning(n+1);