数字信号处理zhang 实验三.doc

实验名称：数字信号处理实验 实验项目：用FFT作谱分析 指导老师： 李秋菊 班级 ： 10电科 姓名： 张卫娟 学号：201000804084 成绩： 一、实验目的：1、 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。2、 熟悉应用FFT对典型信号进行频谱分析的方法。熟悉FFT算法原理和FFT子程序的应用。3、 学习用FFT对连续信号和时域离散信号进行谱分析的方法。了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。2、 实验原理：（一）、在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散傅里叶变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：反变换为：有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。 在信号处理中，DFT的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要通过DFT来实现。然而，当很大的时候，求一个点的DFT要完成次复数乘法和次复数加法，其计算量相当大。1965年J.W.Cooley和J.W.Tukey巧妙地利用因子的周期性和对称性，构造了一个DFT快速算法，即快速傅立叶变换(FFT)。 FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT是以2为基数的，其长度 。它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。（二）、在运用DFT进行频谱分析的过程中可能的产生混叠误差 序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。（三）、matlab函数应用： MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有Fft、Ifft、Fft2 、Ifft2, Fftn、ifftn和Fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。1、 fft和Ifft函数（1）调用方式：Y=fft(X)参数说明如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换；如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；如果X是（N维数组，则是对第一个非单元素的维进行离散傅立叶变换；（2） Yfft(X,N)参数说明N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。（3）Yfft(X,dim) 或Yfft(X,N,dim)参数说明在参数dim指定的维上进行离散傅立叶变换；当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。函数Ifft的参数应用与函数Fft完全相同。应用说明【实例1】fft的应用X=2 1 2 8;Y=fft(X，4)运行结果Y13.0000 0+7.0000i -5.0000 0-7.0000【实例2】fft(X,N,dim)的应用A=2 5 7 8; 1 4 0 5; 3 8 5 1; 9 1 2 7;Z=fft(A,1)2、 Fftshift和Ifftshift函数调用方式Z=fftshift(Y) 此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0处得值）移到频谱的中间位置。【实例3】 fftshift的应用X=rand(5,4);y=fft(X);z=fftshift(y);%只将傅立叶变换结果y中的直流成分移到频谱的中间位置.运行结果： y= 3.2250 2.5277 1.4820 1.63140.3294+0.2368i 0.0768+0.3092i 0.6453+0.4519i -0.7240-0.4116i-0.2867-0.6435i 0.5657+0.4661i -0.5515+0.2297i -0.0573-0.0881i-0.2867+0.6435i 0.5657-0.4661i -0.5515-0.2297i -0.0573+0.0881i0.3294-0.2368i 0.0768-0.3092i 0.6453-0.4519i -0.7240+0.4116iZ=-0.5515-0.2297i -0.0573+0.0881i -0.2867+0.6435i 0.5657-0.4661i0.6453-0.4519i -0.7240+0.4116i 0.3294-0.2368i 0.0768-0.3092i1.4820 1.6314 3.2250 2.52770.6453+0.4519i -0.7240-0.4116i 0.3294+0.2368i 0.0768+0.3092i-0.5515+0.2297i -0.0573-0.0881i -0.2867-0.6435i 0.5657+0.4661i【实例4】fft在信号分析中的应用使用频率分析方法对模拟信号x(t) x=sin(2*pi*100*t)进行频谱分析。采样频率fs=1000，采样点数N=512。并画出信号的时域波形及FFT变换后的幅频响应与相频响应。程序：fs=512;%采样频率N=512;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;%采样时间序列f0=100;%信号频率x=sin(2*pi*f0*t);subplot(3,1,1);plot(t,x);xlabel(t);ylabel(sin(2*pi*100*t);title(时域信号);Y=fft(x,N);%对信号进行FFT变换magY=abs(Y);%求得FFT变换后的幅度angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f,magY);%画出幅频响应xlabel(f);ylabel(幅度);title(N=512);grid;subplot(3,1,3);plot(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=512);grid;结果分析：假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。（例如f=100Hz点处，A=1*N/2=1*512/2=256,而第一个点就是直流分量，它的模值就是直流分量的N倍)。magY=abs(Y);%求得FFT变换后的幅度如果要换算为实际的幅度需要改为：magY=magY/(N/2);magY(1)=magY(1)/2;每个点的相位,就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（假设的第N+1个点）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。第n个点所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。前面的程序由：plot(f(1:N/2),magY(1:N/2);%画出幅频响应三、实验内容：1、 对被白噪声污染的信号进行频谱分析，从中鉴别出有用的信号。要求：将信号的幅度换算成实际的幅度，信号的频率换算成实际的频率程序： 10clear;fs=512;%采样频率N=512;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;%采样时间序列f0=100;%信号频率x=2+3*cos(100*pi*t-pi/6)+1.5*cos(150*pi*t+pi/2)+randn(1,N);subplot(3,1,1);plot(t,x);xlabel(t);ylabel(sin(2*pi*100*t);title(时域信号);Y=fft(x,N);%对信号进行FFT变换magY=abs(Y)/N*2;%求得FFT变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=512);grid;subplot(3,1,3);plot(f(1:N/2),angY(1:N/2);%画出相频响应 xlabel(f);ylabel(相位)title(N=512);grid;图形： 2、 对连续的单一频率周期信号, ,信号频率f=1Hz按采样频率fs=8Hz采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。幅频N=20程序：clear;fs=8;%采样频率N=20;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;%采样时间序列f0=100;%信号频率x=sin(2*pi*t);subplot(3,1,1);stem(t,x);xlabel(t);ylabel(sin(2*pi*100*t);title(时域信号);Y=fft(x,N);%对信号进行FFT变换magY=abs(Y)/(N/2);%求得FFT变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=20);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=20);grid;N=16程序：clear;fs=16;%采样频率N=16;%数据点数n=0:N-1;x=1,2,3,4,4,3,2,1,zeros(1,8);subplot(3,1,1);stem(x);xlabel(n);ylabel(sin(2*pi*100*t);title(时域信号);Y=fft(x);%对信号进行FFT变换magY=abs(Y)/N*2;%求得FFT变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=16);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=16);grid;图形：3、 对三角波序列，进行频谱分析设n=16。程序：clear;fs=16;%采样频率N=16;%数据点数n=0:N-1;x=1,2,3,4,4,3,2,1,zeros(1,8);subplot(3,1,1);stem(x);xlabel(n);ylabel(sin(2*pi*100*t);title(时域信号);Y=fft(x);%对信号进行FFT变换magY=abs(Y)/N*2;%求得FFT变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=16);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=16);grid;图形：4、 对矩形波进行频谱分析。矩形波为脉冲宽度为2s，持续时间为-55s，采用矩形脉冲函数rectpuls(t,2);程序：clear;fs=64;%采样频率N=64;%数据点数n=0:N-1;t=-5:1/fs:5;x=rectpuls