数字信号处理实验报告6_离散傅里叶变换及快速算法

1、实验六 离散傅立叶变换及其快速算法1、 实验目的：掌握快速傅立叶变换的应用方法；掌握离散余弦变换的应用方法；掌握Z变换的应用方法；了解Chip z变换的基本概念；掌握Hilbeit变换的初步应用；了解倒谱变换的基本概念。二、实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。三、实验内容：3.1快速傅立叶变换（FFT）DFT是信号分析与处理中的一种重要变换。但是直接计算DFT的运算量与变换的长度N的平方成正比，当N较大时，计算量太大。在快速傅里叶变换（简称FFT）出现之前，直接用DFT算法进行谱分析和信号的实时处理是不实际的。FFT使得DFT的运算效率大大提高，为数字信号处理技术应用于

2、各种信号的实时处理创造了条件，推动了数字信号处理技术的发展。FFT算法的MATLAB实现MATLAB提供fft函数来计算x(n)的DFT,fft函数是有机器语言，而不是以MATLAB指令格式写成的，因此它的执行速度很快。格式：y=fft(x)，计算信号x的快速傅里叶变换y。当y为矩阵（多通道信号）时，计算x中每一列信号的离散傅里叶变换。当x的长度为2的幂时，有基2算法，否则采用较慢的分裂基算法。Y=fft(x,n)，计算n点FFT，当x的长度大于n时，截断x否则补零。Y=fft(x,n),计算n点FFT，当x的长度大于n时，截断x，否则补零。IFFT可由ifft函数来计算。 在信号处理中，DF

3、T的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要通过DFT来实现。然而，当很大的时候，求一个点的DFT要完成次复数乘法和次复数加法，其计算量相当大。1965年J.W.Cooley和J.W.Tukey巧妙地利用因子的周期性和对称性，构造了一个DFT快速算法，即快速傅立叶变换(FFT)。概念通过前面的知识，已经知道有限列长为的序列的变换为 其逆变换为 由于MATLAB软件本身的特点，序列或向量元素下标从1开始记录，而不是从0开始。因此，上述两式在MATLAB中相应的表达式为 而下面所讨论使用的快速傅立叶变换并不是与不同的另外一种变换，而是为减少计算次数的一种快速有效的算法。这种快速算法，

4、主要是利用了下面两个特性使长序列的分解为更小点数的所实现的。利用的对称性使运算中有些项合并 利用的周期性和对称性使长序列的分解为更小点数的 快速傅立叶变换算法正是基于这一基本思想而发展起来的。快速傅立叶变换算法形式很多，但是基本上可以分为两大类，即按时间抽取（DecimationInTime，简称）法和按频率抽取（DecimationInFrequency）法。在这里，以时间抽取的算法（库利图算法）为例，简单说明一下算法的算法原理。 为了讨论方便，设，其中为整数。如果不满足这个条件，可以认为得加上若干零点来达到。由的定义知 其中是列长为的输入序列，把它按的奇偶分成两个序列 又由于，则 上式表明

5、了一个N点的DFT可以被分解为两个N/2点的DFT。同时，这两个N/2点的DFT按照上式又可以合成为一个N点的DFT。 为了要用点数为N/2点的X(k)、X(k)来表达N点的X(k)值还必须要用W系数的周期性，即 这样可得 即 X（同理可得 X（另外再加上W的对称性 就可以将X(k)的表达式分为前后两个部分：前半部分 后半部分 由以上分析可见，只要求出区间内各个整数k值所对应的X(k)、X(k)的值，即可求出区间内的全部X（k）值，这一点恰恰是FFT能大量节省计算的关键所在。【实例-1】一被噪声污染的信号，很难看出它所包含的频率分量，如一个由50Hz和120Hz正弦信号构成的信号，受到均值随机

6、噪声的干扰，数据采样率为1000Hz，通过FFT来分析其信号频率成分，用MATLAB实现如下：t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+randn(1,length(t);Y=fft(y,512);subplot(211);plot(x);title('受噪声污染的信号');n=0:511;f=1000*n/512;subplot(212);plot(f,abs(Y);title('FFT');【实例-2】用FFT分析语音信号的频谱。结果如下，用MATLAB实现如下。load mtlb;subplot(

7、221);plot(mtlb);title('原始语音信号');y=fft(mtlb);subplot(222);plot(abs(y);title('FFT变换');y(abs(y)<50)=0;x=ifft(y);subplot(223);plot(abs(y);title('去掉幅值小于1的FFT变换');subplot(224);plot(real(x);title('重构语音信号');3.2函数应用MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有Fft、Ifft、Fft2 、Ifft2

8、, Fftn、ifftn和Fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。Fft和Ifft函数调用方式Y=fft(X)参数说明·如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换；·如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；·如果X是（N维数组，则是对第一个非单元素的维进行离散傅立叶变换；Yfft(X,N)参数说明N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。（3）Yff

9、t(X,dim) 或Yfft(X,N,dim)参数说明·在参数dim指定的维上进行离散傅立叶变换；·当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。函数Ifft的参数应用与函数Fft完全相同。应用说明【实例6-1】fft的应用X=2 1 2 8;Y=fft(X)运行结果Y13.0000 0+7.0000i -5.0000 0-7.0000【实例6-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)【实例6-3】fft在信号分析中的应用使用频率分

10、析方法从受噪声污染的信号x(t)中鉴别出有用的信号。程序：t=0:0.001:1; %采样周期为0.001s,即采样频率为1000Hz; %产生受噪声污染的正县正弦波信号；x=sin(2*pi*100*t)+sin(2*pi*200*t)+rand(size(t);subplot(2,1,1)plot(x(1:50); %画出时域内的信号；Y=fft(x,512); %对X进行512点的傅立叶变换；f=1000*(0:256)/512; %设置频率轴（横轴）坐标，1000为采样频率；subplot(2,1,2)plot(f,Y(1:257); %画出频域内的信号运行结果：图6-1 时域信号和频

11、域信号的比较 由上面的图6-1可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。但是通过对x(t)作傅立叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100Hz和200Hz的两个频率分量。【实例6-4】x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。%fft.m%x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。clf;fs=100;N=128; %采样频率和数据点数n=0:N-1;t=n/f

12、s; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求得Fourier变换后的振幅f=n*fs/N; %频率序列subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;subplot(2,2,2),plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅

13、xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;%对信号采样数据为1024点的处理fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求取Fourier变换的振幅f=n*fs/N;subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel(&

14、#39;振幅');title('N=1024');grid on;subplot(2,2,4)plot(f(1:N/2),mag(1:N/2); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：（1）数据个数N=32，FFT所用的采样点数NFFT=32；（2）N=32，NFFT=128；（3）N=136，NFFT=12

15、8；（4）N=136，NFFT=512。clf;fs=100; %采样频率Ndata=32; %数据长度N=32; %FFT的数据长度n=0:Ndata-1;t=n/fs; %数据对应的时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %时间域信号y=fft(x,N); %信号的Fourier变换mag=abs(y); %求取振幅f=(0:N-1)*fs/N; %真实频率subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel(&

16、#39;振幅');title('Ndata=32 Nfft=32');grid on;Ndata=32; %数据个数N=128; %FFT采用的数据长度n=0:Ndata-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:N-1)*fs/N; %真实频率subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅'

17、);title('Ndata=32 Nfft=128');grid on;Ndata=136; %数据个数N=128; %FFT采用的数据个数n=0:Ndata-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:N-1)*fs/N; %真实频率subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');title(&

18、#39;Ndata=136 Nfft=128');grid on;Ndata=136;%数据个数N=512; %FFT所用的数据个数n=0:Ndata-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);mag=abs(y);f=(0:N-1)*fs/N; %真实频率subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');title('Ndata

19、=136 Nfft=512');grid on;结论：（1）当数据个数和FFT采用的数据个数均为32时，频率分辨率较低，但没有由于添零而导致的其他频率成分。（2）由于在时间域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。（3）FFT程序将数据截断，这时分辨率较高。（4）也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT振幅谱也基本不受影响。     对信号进行频谱分析时，数据样本应有足够的长度，一般FFT程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的质量，可减小因补零或截断

20、而产生的影响。3.3 Fft2和Ifft函数调用方式（1）Y=fft2(X)参数说明·如果X是向量，则此傅立叶变换即变成一维傅立叶变换fft；·如果X是矩阵，则是计算该矩阵的二维快速傅立叶变换；·数据二维傅立叶变换fft2（X）相当于fft(fft(X)，即先对X的列做一维傅立叶变换，然后再对变换结果的行做一维傅立叶变换。（2） Yfft2(X,M,N) 通过对X进行补零或截断，使得成为（M*N）的矩阵。函数Ifft2的参数应用与函数Fft2完全相同。Fftn.、Ifftn是对数据进行多维快速傅立叶变换，其应用与Fft2、Ifft2类似；在此，不再一一赘述。【实例

21、6-4】fft2、ifft2的应用A=2 5 7 8 9; 1 3 7 5 0; 2 6 1 4 9; 8 1 5 2 6;Y=fft2(A);B=ifft2(Y);运行结果：3.4 Fftshift和Ifftshift函数调用方式Z=fftshift(Y) 此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0处得值）移到频谱的中间位置。参数说明·如果X是向量，则变换Y的左右两边；·如果X是矩阵，则交换Y的一、三象限和二、四象限；·如果Y是多维数组，则在数组的每一维交换其“半空间”。 函数Iffshift的参数应用与函数Fftshift完全相同。【实

22、例6-5】 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

23、.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.0881

24、i -0.2867-0.6435i 0.5657+0.4661i3. 5离散傅立叶变换（DFT） 逆变换IFFT为 通过下列修改，就可以用FFT算法来实现逆FFT运算：·增加一个归一化因子；·将用其复数共轭代替但是由因为， 所以，求X（k）的逆FFT可以分为以下3个步骤：（1）取X（k）的共扼，得到；（2）求 的FFT，得N；（3）取 的共扼，并除以N，即得到x(n)。采用这种方法可以直接用FFT程序来计算逆FFT。有关IFFT的具体应用，与FFT一一对应，在此不再赘述。【实例-1】已知一模拟信号,现以采样率进行采样。用DFT计算当序列长度L=100,L=20，N=200点

25、的幅度频谱样值并通过作图与理论上准确的频谱样值进行比较。解：原信号的傅里叶变换：其幅度值为MatLab程序如下：clc;clear;fs=20;L=100;N=200;n=0:L-1;t1=n/fs;xn1=exp(-t1);xn=xn1,zeros(1,N-1);Xk1=fft(xn,N);magXk1=abs(Xk1);k1=(0:length(magXk1)'-1)*N/length(magXk1);L=20;N=200;n=0:L-1;t2=n/fs;xn2=exp(-t2);xn=xn2,zeros(1,N-L);Xk2=fft(xn,N);magXk2=abs(Xk2);k

26、2=(0:length(magXk2)'-1)*N/length(magXk2);figure(1);subplot(211);plot(t1,xn1);title('xa(t) t=5s');subplot(212);plot(t2,xn2);title('xa(t) t=1s');figure(2);subplot(211);plot(k1,magXk1);title('X(k) L=100 N=200');subplot(212);plot(k2,magXk2);title('X(k) L=20 N=200');fi

27、gure(3);Omeger=0:0.1:20*pi;Xa=1./(1+Omeger.2);subplot(111);plot(Omeger/pi,Xa);title('|Xa(j/Omega)|');计算结果如图所示，结论如下：当序列长度为100，进行200点DFT计算的结果混叠与泄漏的影响比较小，基本上接近原来信号的频谱。因为按给定的,相当于取信号的最高频率，故在0,频率范围的能量为：而信号的总能量为：所以，基本上满足频谱不混叠的要求。当序列长度为20，进行200点DFT计算，由于截取x(n)长度太短。，所以频谱因泄漏出现较大的波动，以致与原信号频谱有较大差别。用DFT对离

28、散信号进行谱分析 序列x(n)在单位圆上的Z变换就是其傅里叶变换,即。对序列x(n)进行N点DFT得到X(k)，X(k)是在区间上的N点等间隔采样。因此序列的傅里叶变换可利用DFT来计算。实际工作中往往要把信号的观察时间限制在一定的时间间隔内，就需要取用有限个数据，即将信号截断，相当于将信号乘以窗函数。在时域的截断，相当于所研究的频谱与窗函数频谱的周期卷积，周期卷积的结果造成频谱从其正常的频谱扩展开来，称为泄漏。泄漏效应是DFT所固有的，可用加窗函数进行抑制。加窗函数后，信号成为：频谱为：频谱的等间隔取样值为：用窗函数加权使有限长度的输入信号周期延拓后，在边界上尽量减少不连续程度，从而达到抑制

29、频谱泄漏的目的。在实际应用中常用单边表示的窗函数，即窗函数这种将加权序列向右移动N/2点，只影响相伴特性，并不影响振幅特性。【实例-2】信号,用加窗函数的DFT分析其频谱。figure(1);L=20;N=200;n=0:L-1;xn=exp(-0.05*n);xn1=xn,zeros(1,N-L);Xk=fft(xn1,N);magXk=abs(Xk);k=(0:length(magXk)'-1)*N/length(magXk);subplot(211);stem(k,xn1);title('矩形窗x(n)');subplot(212);stem(k,magXk);t

30、itle('矩形窗X(k)');figure(2);w=(hanning(L)'xn1=xn.*w;xn2=xn1,zeros(1,N-L);subplot(211);stem(k,xn2);title('汉宁窗x(n)');Xk=dft(xn2,N);magXk=abs(Xk);subplot(212);stem(k,magXk);title('汉宁窗X(k)');figure(3);w=(triang(L)'xn1=xn.*w;xn2=xn1,zeros(1,N-L);subplot(211);stem(k,xn2);titl

31、e('三角窗x(n)');Xk=fft(xn2,N);magXk=abs(Xk);subplot(212);stem(k,magXk);title('三角窗'); 从图可见，加汉宁窗、三角窗的信号的频谱的泄漏要比矩形窗的泄漏小很多，通过选择合适的窗函数，达到抑制泄漏的目的。【实例-2】设1) 分解x(n)成和（奇偶部分）;2) 检验实序列的性质%实序列的奇偶分解function xec,xoc=circevod(x)if any(imag(x)=0) error('x为非实数序列');EndN=length(x);n=0:N-1;xec=0.5*

32、(x+x(mod(-n,N)+1);xoc=0.5*(x-x(mod(-n,N)+1);clc;clear;N=20;n=0:N-1;x=5*(0.8).n;xec,xoc=circevod(x);X=fft(x,N);Xec=fft(xec,N);Xoc=fft(xoc,N);figure(1);subplot(221);stem(n,x);title('x(n)');subplot(222);stem(n,abs(X);title('|X(k)|');subplot(223);stem(n,real(X);title('ReX(k)');su

33、bplot(224);stem(n,imag(X);title('ImX(k)');figure(2);subplot(221);stem(n,xec);title('xec(n)');subplot(222);stem(n,xoc);title('xoc(n)');subplot(223);stem(n,real(Xec);title('DFTxec(n)');subplot(224);stem(n,imag(Xoc);title('DFTxoc(n)');从此例可见，实序列的偶分量关于N/2点对称，奇分量关于N

34、/2点反对称；偶分量的DFT等于实序列的DFT的实部，奇分量的DFT等于实序列的DFT的虚部。3.6 圆周移位与圆周卷积圆周移位：function y = cirshftt(x,m,N) % 长度为 N 的x序列: (时域)作m采样点圆周移位 % - % y = cirshftt(x,m,N) % y = 包含圆周移位的输出序列 % x = 长度 % m = 移位采样数 % N = 圆周缓冲器长度 % 方法: y(n) = x(n-m) mod N) % Check for length of x if length(x) > N error('N 必须 >= x的长度&#

35、39;) end x = x zeros(1,N-length(x); n = 0:1:N-1; n = mod(n-m,N); y = x(n+1); 圆周卷积：function y = circonvt(x1,x2,N) % 在x1 和 x2: (时域)之间的N点圆周卷积 % - % y = circonvt(x1,x2,N) % y = 包含圆周卷积的输出序列 % x1 = 长度 N1 % x2 = 长度 N2 % N = 圆周缓冲器的大小 % 方法 y(n) = sum (x1(m)*x2(n-m) mod N) % Check for length of x1 if length(x

36、1) > N error('N 必须 >= x1的长度') end % Check for length of x2 if length(x2) > N error('N 必须 >= x2的长度') end x1=x1 zeros(1,N-length(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,:) = cirshftt(x2,n-1,N); end y = x1*H'

37、【例】序列x(n)=9,8,7,6,5,4,3,2,1,求分别移1位、3位、5位、7位和9位的圆周移位。n=0:8;x=9,8,7,6,5,4,3,2,1;y1=cirshftt(x,1,9);y3=cirshftt(x,3,9);y5=cirshftt(x,5,9);y7=cirshftt(x,7,9);y9=cirshftt(x,9,9);subplot(321);stem(n,x);ylabel('x(n)');subplot(322);stem(n,y1);ylabel('y1(n)');subplot(323);stem(n,y3);ylabel(&#

38、39;y3(n)');subplot(324);stem(n,y5);ylabel('y5(n)');subplot(325);stem(n,y5);ylabel('y7(n)');subplot(326);stem(n,y9);ylabel('y9(n)');x1=1,2,2,3;x2=1,2,3,4,3,2;N=length(x1)+length(x2)-1;n=0:N-1;n1=0:N-3;n2=0:N-2;y1=circonvt(x1,x2,N-2);y2=circonvt(x1,x2,N-1);y3=circonvt(x1,x2,N);x1=x1 zeros(1,N-length(x1);x2=x2 zeros(1,N-length(x2);X1=dft(x1,N);X2=dft(x2,N);X=X1.*X2;x=idft(X,N);x=real(x);subplot(231);stem(n,x1);title('x1(n)');subplot(232);