离散傅立叶变换及其快速算法.doc

第5章 离散傅立叶变换及其快速算法本章介绍了数字信号变换技术的详细内容，主要有：离散傅立叶变换、快速傅立叶变换（FFT）、离散余弦变换、Z变换5.1试验目的 理解信号变换的基本概念 理解离散傅立叶变换的基本概念 掌握快速傅立叶变换的应用方法 掌握离散余弦变换的应用方法 掌握Z变换的应用方法 了解Chip z变换的基本概念 掌握Hilbeit变换的初步应用 了解倒谱变换的基本概念5.2 信号变换概述 信号是数字信号处理领域中最基本、最重要的概念。而数字信号变换技术，又是对信号进行处理操作的最基本的有效途径之一。因此，数字信号变换技术，便成为数字信号处理领域中专业人员所必须要张我的一项最基本的技能。 简单地说，数字信号变换技术就是为了处理操作上的方便和可能，通过数学变换，将一个域内的信号变换映射倒另一个域内的信号的方法。常用的数字信号变换主要有：傅立叶变换、离散余弦变换（DCT）、Z变换、Chirp z变换、Hilbert变换等。这些变换，都有着各自的理论和其应用背景。MATLAB中的工具箱对这几种典型的变换，都提供了相对应的、具体的应用函数。这可以使得工程人员大大节省无谓的工作量，从而将主要精力放到新技术的创新和研发上面。下面将对这几种变换的含义和应用进行具体的介绍。5.3 离散傅立叶变换傅立叶变换是信号分析和处理的重要工具。有限长序列作为离散信号的一种，在数字信号处理种占有着极其重要的位置。对于有限长序列，离散傅立叶变换不仅在理论上有着重要的意义，而且有快速计算的方法快速傅立叶变换。所以在各种数字信号处理的运算方法中，越来越起到核心的作用。下面，就对离散傅立叶变换及其MATLAB函数应用，结合实际工程实例做说明。5.3.1 傅立叶变换的几种形式1、非周期连续时间信号的傅立叶变换非周期连续时间信号的傅立叶变换可以表示为逆变换为在这里，是模拟角频率。可以看到，时域的连续函数造成频域的非周期谱，时域的非周期性造成频域的连续谱。结论：非周期连续时间函数对应于一非周期连续频域变换函数。2、周期连续时间信号的傅立叶变换周期为的周期性连续时间信号傅立叶变换是离散频域函数，可表示为逆变换为这就是经常称之为傅立叶级数的变换形式。在这里，也是模拟角频率。可以看到，时域的连续函数造成频率域的非周期谱，频域函数的离散造成时域函数的周期性。结论：周期连续时间函数对应于一非周期离散频域变换函数。3、非周期离散时间信号的傅立叶变换可以表示为逆变换为在这里，是数字频率，它和模拟角频率的关系为。可以看到，时域的取样对应于频域的周期延拓，而时域函数的非周期性造成频域的离散谱。结论：非周期离散时间函数对应于一周期连续频域变换函数。4、周期离散时间信号的傅立叶变换周期离散时间信号的傅立叶变换离散傅立叶变换，可以表示为逆变换为可以看到，时域的取样对应于频域的周期延拓，而时域函数的周期性造成频域的离散谱。结论：周期离散时间函数对应于一周期离散频域变换函数。5.3.2 离散傅立叶变换离散傅立叶级数变换是周期序列，仍不便于计算机计算。但离散傅立叶级数虽是周期序列，却只有个独立的数值，所以它的许多特性可以通过有限长序列延拓来得到。对于一个长度为的有限长序列，也即只在个点上有非零值，其余皆为零，即把序列以为周期进行周期延拓得到周期序列，则有所以，有限长序列的离散傅立叶变换（DFT）为逆变换为若将DFT变换的定义写成矩阵形式，则得到X=Ax，其中DFT变换矩阵A为Dftmtx 函数：用来计算DFT变换矩阵A的函数调用方式（1） Adftmta（n）：返回nn的DFT变换矩阵A。若x为给定长度的行向量，则yx*A，返回x的DFT变换y。（2） Aiconj（dftmtx（n）/n；返回nn的IDFT变换矩阵Ai。应用说明【实例5-1】 A=dftmtx(4) Ai=conj(dftmtx(4)/4运行结果A = 1.0000 1.0000 1.0000 1.0000 1.0000 0 - 1.0000i -1.0000 0 + 1.0000i 1.0000 -1.0000 1.0000 -1.0000 1.0000 0 + 1.0000i -1.0000 0 - 1.0000iAi = 0.2500 0.2500 0.2500 0.2500 0.2500 0 + 0.2500i -0.2500 0 - 0.2500i 0.2500 -0.2500 0.2500 -0.2500 0.2500 0 - 0.2500i -0.2500 0 + 0.2500i【实例5-2】如果是一个N16的有限序列，用MATLAB求其DFT的结果，并画出其结果图，如图51所示。图 51 有限长序列的DFT结果图程序N=16;n=0:1:N-1; %时域采样xn=sin(n*pi/8)+sin(n*pi/4);k=0:1:N-1; %频域采样WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;subplot(2,1,1)stem(n,xn);subplot(2,1,2)stem(k,abs(Xk);运算结果Xk = Columns 1 through 5 0.0000 -0.0000 - 8.0000i -0.0000 - 8.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 6 through 10 -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 11 through 15 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 + 8.0000i Column 16 0.0000 + 8.0000iDFT的性质两个序列和都是N点有限长序列，设1、 线性式中a，b为任意常数。2、 圆周移位一个有限长序列的圆周移位定义式中，表示的周期延拓序列的移位有限长序列圆周移位后的DFT为【实例5-3】求有限长序列的圆周移位。并画出其结果图，如图52所示。图 5-2 有限长序列的圆周移位结果图程序N=20;m=10;n=0:1:N-1;x=8*(0.4).n;n1=mod(n+m),N);xm=x(n1+1);subplot(2,1,1)stem(n,x);title(Original Sequence);xlabel(n);ylabel(x(n);subplot(2,1,2)stem(n,xm);title(Circular Shift Sequence);xlabel(n);ylabel(x(n+10)mod20);输出结果：x = Columns 1 through 8 8.0000 3.2000 1.2800 0.5120 0.2048 0.0819 0.0328 0.0131 Columns 9 through 16 0.0052 0.0021 0.0008 0.0003 0.0001 0.0001 0.0000 0.0000 Columns 17 through 20 0.0000 0.0000 0.0000 0.00003、圆周卷积假设则有用表示圆周卷积，则上式可化简为MATLAB内部用于计算圆周卷积的函数Circonvt程序如下：function y=circonvt(x1,x2,N)if(length(x1)N) error(N should bigger than or equal to the length of x1!)endif(length(x2)N) error(N should bigger than or equal to the length of x2!)x1=x1,zeros(1,N-length(x1)x2=x2,zeros(1,N-length(x2)m=0:1:N-1for n=1:1:Nendy=x1*H4.共轭对称性 令的共轭复数序列为，则 用和分别表示序列的实部和虚部，即 用和分别表示实部和虚部序列的，即 而且可以证明得到 通常称为的共轭偶部，为的共轭奇部。所以说，对于时域、频域的对应关系来说，序列实部对应于的共轭偶部，序列的虚部对应于的共轭奇部。5.序列乘积 6.DFT形式下的帕塞瓦尔定理 5.4 快速傅立叶变换（FFT） 在信号处理中，DFT的计算具有举足轻重的地位，信号的相关、滤波、谱估计等都要通过DFT来实现。然而，当很大的时候，求一个点的DFT要完成次复数乘法和次复数加法，其计算量相当大。1965年J.W.Cooley和J.W.Tukey巧妙地利用因子的周期性和对称性，构造了一个DFT快速算法，即快速傅立叶变换(FFT)。5.4.1概念通过前面的知识，已经知道有限列长为的序列的变换为 其逆变换为 由于MATLAB软件本身的特点，序列或向量元素下标从1开始记录，而不是从0开始。因此，上述两式在MATLAB中相应的表达式为 而下面所讨论使用的快速傅立叶变换并不是与不同的另外一种变换，而是为减少计算次数的一种快速有效的算法。这种快速算法，主要是利用了下面两个特性使长序列的分解为更小点数的所实现的。（1） 利用的对称性使运算中有些项合并 （2） 利用的周期性和对称性使长序列的分解为更小点数的 快速傅立叶变换算法正是基于这一基本思想而发展起来的。快速傅立叶变换算法形式很多，但是基本上可以分为两大类，即按时间抽取（DecimationInTime，简称）法和按频率抽取（DecimationInFrequency）法。在这里，以时间抽取的算法（库利图算法）为例，简单说明一下算法的算法原理。 为了讨论方便，设，其中为整数。如果不满足这个条件，可以认为得加上若干零点来达到。由的定义知 其中是列长为的输入序列，把它按的奇偶分成两个序列 又由于，则 上式表明了一个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能大量节省计算的关键所在。5.4.2 函数应用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） （3）Yfft(X,dim) 或Yfft(X,N,dim)参数说明在参数dim指定的维上进行离散傅立叶变换；当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。函数Ifft的参数应用与函数Fft完全相同。应用说明【实例5-4】fft的应用X=2 1 2 8;Y=fft(X)运行结果Y13.0000 0+7.0000i -5.0000 0-7.0000【实例5-5】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)【实例5-6】fft在信号分析中的应用使用频率分析方法从受噪声污染的信号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); %画出频域内的信号运行结果图5-3 时域信号和频域信号的比较 由上面的图5-3可以看出，从受噪声污染信号的时域形式中，很难看出正弦波的成分。但是通过对x(t)作傅立叶变换，把时域信号变换到频域进行分析，可以明显看出信号中100Hz和200Hz的两个频率分量。2. 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类似；在此，不再一一赘述。【实例5-7】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. Fftshift和Ifftshift函数调用方式Z=fftshift(Y) 此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0处得值）移到频谱的中间位置。参数说明如果X是向量，则变换Y的左右两边；如果X是矩阵，则交换Y的一、三象限和二、四象限；如果Y是多维数组，则在数组的每一维交换其“半空间”。 函数Iffshift的参数应用与函数Fftshift完全相同。【实例5-8】 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.286