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),
3、计算 n 点 FFT，当 x 的长度大于 n 时，截断 x，否则补零。IFFT 可由 ifft 函数来计算。在信号处理中，DFT 的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要通过 DFT 来实现。然而，当 很大的时候，求一个 点的 DFT 要完成 次复数乘NNN法和 次复数加法，其计算量相当大。1965 年)1(J.W.Cooley 和 J.W.Tukey 巧妙地利用 因子的周期性和对NW称性，构造了一个 DFT 快速算法，即快速傅立叶变换(FFT) 。概念通过前面的知识，已经知道有限列长为 的序列 的N)(nx变换为DFT.nkNnWxkX10)()( 12,0Nk其逆变换为10
4、)()(Nknknx 1,0n由于 MATLAB 软件本身的特点，序列或向量元素下标从 1 开始记录，而不是从 0 开始。因此，上述两式在MATLAB 中相应的表达式为nkNnWxkX10)()( 12,0Nk1)()(knkNx ,n而下面所讨论使用的快速傅立叶变换 并不是与 不)(FTDFT同的另外一种变换，而是为减少 计算次数的一种快速有D效的算法。这种快速算法，主要是利用了 下面两个特性nkNW使长序列的 分解为更小点数的 所实现的。DFTFT利用 的对称性使 运算中有些项合并nkNW)()( knNknNnkNW利用 的周期性和对称性使长序列的 分解为更小点数nkN DFT的 DFT
5、)()( nNkNnknkNW快速傅立叶变换 算法正是基于这一基本思想而发展FT起来的。快速傅立叶变换算法形式很多，但是基本上可以分为两大类，即按时间抽取（DecimationInTime，简称.）法和按频率抽取（DecimationInFrequency）法。DIT在这里，以时间抽取 的 算法（库利图算法）为例，)(DITF简单说明一下 算法的算法原理。F为了讨论方便，设 ，其中 为整数。如果不满足这个2N条件，可以认为得加上若干零点来达到。由 的定义知DFTnkNnWxkX10)()( 12,0Nk其中 是列长为 的输入序列，把它按 的奇)1,0( n偶分成两个序列)()12(2rxrx
6、12,0Nr又由于 ，2NjNj WeW则 10 2110 )()()()()( n kNnkNkn XWxxkX为 奇 数为 偶 数上式表明了一个 N 点的 DFT 可以被分解为两个 N/2 点的DFT。同时，这两个 N/2 点的 DFT 按照上式又可以合成为一个 N 点的 DFT。为了要用点数为 N/2 点的 X (k)、X (k)来表达 N 点的 X(k)值12还必须要用 W 系数的周期性，即)(2krNrkW这样可得202021)(211 )()(NrNrrkkrxxkX.即 X （1)(21kXN同理可得 X （2)()2另外再加上 W 的对称性kNkNkNkNW2)2(就可以将 X
7、(k)的表达式分为前后两个部分：前半部分 )()(21kXkXN12,0Nr后半部分 )2()2()()2(1 kXWkXkkN1kN2,10Nr由以上分析可见，只要求出区间 内各个整数 k 值12,0N所对应的 X (k)、X (k)的值，即可求出 区间内的全部12 ,X（k）值，这一点恰恰是 FFT 能大量节省计算的关键所在。【实例-1】一被噪声污染的信号，很难看出它所包含的频率分量，如一个由 50Hz 和 120Hz 正弦信号构成的信号，受到均值随机噪声的干扰，数据采样率为 1000Hz，通过 FFT 来分析其信号频率成分，用 MATLAB 实现如下：t=0:0.001:0.6;x=si
8、n(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(221);plot(mtlb);title(原始语音信号);.y=fft(mtlb);subplot(222);plot(abs(y);title(FFT变换
9、 );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, .Fftn、ifftn 和 Fftshift、Ifftshift 等。当所处理的数据的长度为 2 的幂次时，采用基-2 算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为 2 的幂次或者用添零的方式来添
10、补数据使之成为 2 的幂次。Fft 和 Ifft 函数调用方式Y=fft(X)参数说明如果 X 是向量，则采用傅立叶变换来求解 X 的离散傅立叶变换；如果 X 是矩阵，则计算该矩阵每一列的离散傅立叶变换；如果 X 是（N 维数组，则是对第一个非单元素的维进行)D离散傅立叶变换；Yfft(X,N)参数说明N 是进行离散傅立叶变换的 X 的数据长度，可以通过对 X进行补零或截取来实现。（3）Yfft(X,dim) 或 Yfft(X,N,dim)参数说明在参数 dim 指定的维上进行离散傅立叶变换；当 X 为矩阵时，dim 用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2 表明变换按
11、行进行。.函数 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 在信号分析中的应用使用频率分析方法从受噪声污染的信号 x(t)中鉴别出有用的信号。程序：t=0:0.001:1; %采样周期为 0.001s,即采样频率为 1000Hz;%产生受噪声污染的正县正弦波信号；x=sin(2*pi*1
12、00*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 时域信号和频域信号的比较由上面的图 6-1 可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。但是通过对 x(t)作傅立叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100H
13、z 和 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/fs; %时间 序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速 Fourier
14、变换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频率之前随频率变化的振幅xlabel(频率 /Hz);ylabel(振幅 );title(N=128);grid on;%对信号采样数据为 1024 点的处理fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*s
15、in(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(振幅 );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 o
16、n;.例 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=128 ；（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=
17、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(振幅 );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; %真实频率s
18、ubplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出 Nyquist频率之前的振幅xlabel(频率 /Hz);ylabel(振幅 );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:
19、N/2)*2/N); %绘出 Nyquist频率之前的振幅xlabel(频率 /Hz);ylabel(振幅 );title(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频率之前的振幅xlabe
20、l(频率 /Hz);ylabel(振幅 );title(Ndata=136 Nfft=512);grid on;.结论：（1）当数据个数和 FFT 采用的数据个数均为 32 时，频率分辨率较低，但没有由于添零而导致的其他频率成分。（2）由于在时间域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。（3）FFT 程序将数据截断，这时分辨率较高。（4）也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT 振幅谱也基本不受影响。对信号进行频谱分析时，数据样本应有足够的长度，一般 FFT 程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的
21、质量，可减小因补零或截断而产生的影响。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
22、类似；在此，不再一一赘述。【实例 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 的参数应用与函数
23、Fftshift 完全相同。【实例 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.229
24、7i -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.4661i.0.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.
25、2297i -0.0573-0.0881i -0.2867-0.6435i 0.5657+0.4661i3. 5 离散傅立叶变换（DFT）逆变换 IFFT 为nkNNkWXnx10)()(12,0Nr通过下列修改，就可以用 FFT 算法来实现逆 FFT 运算：增加一个归一化因子 ；N1将用 其复数共轭 代替nkNWnk但是由因为，.*10* )(1)()( kXFTNWnXNnxkNk 所以，求 X（k）的逆 FFT 可以分为以下 3 个步骤：（1）取 X（k）的共扼，得到 ；)(*kX（2）求 的 FFT，得 N ；)(* )(*x（3）取 的共扼，并除以 N，即得到 x(n)。)(*nx采
26、用这种方法可以直接用 FFT 程序来计算逆 FFT。有关IFFT 的具体应用，与 FFT 一一对应，在此不再赘述。【实例-1】已知一模拟信号 ,现以采样率()()taxeu进行采样。用 DFT 计算当序列长度20sfHzL=100,L=20，N=200 点的幅度频谱样值并通过作图与理论上准确的频谱样值进行比较。解：原信号的傅里叶变换： 0 1()()jt jtaaXjxtededj其幅度值为21|()|aXj.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=
27、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);k2=(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);figu
28、re(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);figure(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 计算的结果混叠与泄漏的影响比较小，基本上接近原来信号的频谱。因为按给定的,相当于取信号的最高频率
29、，故在0, 频率20sfHz 10hfHzhf范围的能量为： 而信号2021|()|.495hEXjd的总能量为： 所以2|()|0.2xj，基本上满足频谱不混叠的要求。9%hxE.当序列长度为 20，进行 200 点 DFT 计算，由于截取 x(n)长度太短。 ，所以频谱因泄漏1()| 0.3790LTtxe出现较大的波动，以致与原信号频谱有较大差别。用 DFT 对离散信号进行谱分析 序列 x(n)在单位圆上的 Z 变换就是其傅里叶变换 ,即 。()jXe()(|jjzeeX对序列 x(n)进行 N 点 DFT 得到 X(k)，X(k)是 在区间()j上的 N 点等间隔采样。因此序列的傅里叶
30、变换可利用0,2DFT 来计算。实际工作中往往要把信号的观察时间限制在一定的时间间隔内，就需要取用有限个数据，即将信号截断，相当于将信号乘以窗函数。在时域的截断，相当于所研究的频谱与窗函数频谱的周期卷积，周期卷积的结果造成频谱从其正常的频谱扩展开来，称为泄漏。泄漏效应是 DFT 所固有的，可用加窗函数进行抑制。加窗函数后，信号成为： x()()wnxwn频谱为： ()()*jjjwXeWe频谱的等间隔取样值为： ()*()Xkk用窗函数加权使有限长度的输入信号周期延拓后，在边界上尽量减少不连续程度，从而达到抑制频谱泄漏的目的。.在实际应用中常用单边表示的窗函数，即窗函数这种将加权序列向右移动
31、N/2 点，只影响相w(n)0,1.N伴特性，并不影响振幅特性。【实例-2】信号 ,用加窗函数的 DFT 分析其频谱。0.5()nxefigure(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);title(矩形窗 X(k);figure(2);w=(hanning
32、(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);title(三角窗 x(n);Xk=fft(xn2,N);magXk=abs(Xk);.subplot(212);stem(k,magXk);title
33、(三角窗);.从图可见，加汉宁窗、三角窗的信号的频谱的泄漏要比矩形窗的泄漏小很多，通过选择合适的窗函数，达到抑制泄漏的目的。【实例-2】设 ()5(0.8)20nxn1)分解 x(n)成 和 （奇偶部分）;epop2)检验实序列的性质 ()ReX(k)DFTImepoxn%实序列的奇偶分解function xec,xoc=circevod(x)if any(imag(x)=0)error(x为非实数序列);End.N=length(x);n=0:N-1;xec=0.5*(x+x(mod(-n,N)+1);xoc=0.5*(x-x(mod(-n,N)+1);clc;clear;N=20;n=0:
34、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);subplot(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);
