实验三 快速傅立叶变换及其应用

实验三快速傅立叶变换及其应用一、实验目的了解计算DFT算法存在的问题及改进途径。掌握几种DFT算法（时间抽取算法DIT算法，频率抽取算法DIF算法，线性调频Z变换即CZT法）。学习并掌握FFT的应用。二、实验原理有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列.但其计算量太大(与N的平方成正比）,很难实时地处理问题,因此引出了快速傅里叶变换(FFT)。FFT并不是一种新的变换形式,它只是DFT的一种快速算法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法.DFT的快速算法—FFT是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。时间抽选FFT算法的理论推导和流图详见《数字信号处理》教材。该算法遵循两条准则：（1）对时间奇偶分；（2）对频率前后分。这种算法的流图特点是：（1）基本运算单元都是蝶形任何一个长度为N=2M的序列，总可通过M次分解最后成为2点的DFT计算。如图所示：WNk称为旋转因子计算方程如下：Xm+1(p)=Xm(p)+WNkXm(q)Xm+1(q)=Xm(p)-WNkXm(q)（2）同址（原位）计算这是由蝶形运算带来的好处，每一级蝶形运算的结果Xm+1(p)无须另外存储，只要再存入Xm(p)中即可，Xm+1(q)亦然。这样将大大节省存储单元。（3）变址计算输入为“混序”（码位倒置）排列，输出按自然序排列，因而对输入要进行“变址”计算（即码位倒置计算）。“变址”实际上是一种“整序”的行为，目的是保证“同址”。FFT的应用凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。FFT主要应用在1、快速卷积2、相关3、频谱分析快速傅立叶变换的MATLAB实现提供fft函数计算DFT格式X=fft(x)X=fft(x,N)如果x的长度小于N，则在其后填零使其成为N点序列，若省略变量N，则DFT的长度即为x的长度。如果N为2的幂，则得到高速的基-2FFT算法；若N不是2的乘方，则为较慢的混合算法。如果x是矩阵，则X是对矩阵的每一列向量作FFT。由题目可得x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)fs=100N=128/1024例：已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成，数据采样频率为100Hz,试绘制N=128点DFT的幅频图。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);f=(0:length(y)-1)'*fs/length(y);mag=abs(y);stem(f,mag);title(‘N=128点’)利用FFT进行功率谱的噪声分析已知带有测量噪声信号其中f1=50Hz,f2=120Hz,为均值为零、方差为1的随机信号，采样频率为1000Hz，数据点数N=512。试绘制信号的频谱图和功率谱图。12()sin(2)sin(2)2()xtftftt()tt=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t));Y=fft(y,512);P=Y.*conj(Y)/512;%求功率f=1000*(0:255)/512;subplot(2,1,1);plot(y);subplot(2,1,2);plot(f,P(1:256));序列长度和FFT的长度对信号频谱的影响。已知信号其中f1=15Hz,f2=40Hz,采样频率为100Hz.在下列情况下绘制其幅频谱。Ndata=32,Nfft=32;Ndata=32,Nfft=128;12()0.5sin(2)2sin(2)xtftftfs=100;Ndata=32;Nfft=32;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,1)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=32’)Nfft=128;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)’*fs/length(y);subplot(2,1,2)plot(f(1:Nfft/2),mag(1:Nfft/2))title(‘Ndata=32,Nfft=128’)快速傅立叶逆变换（IFFT）函数调用格式y=ifft(x)y=ifft(x,N)当N小于x长度时，对x进行截断，当N大于x长度时，对x进行补零。对信号进行DFT，对其结果进行IDFT，并将IDFT的结果和原信号进行比较。f1=40Hzf2=15HzFs=100Hz12()sin(2)sin(2)xtftftfs=100;N=128;n=0:N-1;t=n/fs;x=sin(2*pi*40*t)+sin(2*pi*15*t);subplot(2,2,1)plot(t,x)title('originalsignal')y=fft(x,N);mag=abs(y);f=(0:length(y)-1)'*fs/length(y);subplot(2,2,2)plot(f,mag)title('FFTtooriginalsignal')xifft=ifft(y);magx=real(xifft);ti=[0:length(xifft)-1]/fs;subplot(2,2,3)plot(ti,magx);title('signalfromIFFT')yif=fft(xifft,N);mag=abs(yif);subplot(2,2,4)plot(f,mag)title('FFTtosignalfromIFFT')线性卷积的FFT算法在MATLAB实现卷积的函数为CONV，对于N值较小的向量，这是十分有效的。对于N值较大的向量卷积可用FFT加快计算速度。由DFT性质可知，若DFT[x1(n)]=X1(k),DFT[x2(n)]=X2(k)则若DFT和IDFT均采用FFT和IFFT算法，可提高卷积速度。121212()*()[()()][[()][()]]xnxnIDFTXkXkIDFTDFTxnDFTxn计算x1(n)和x2(n)的线性卷积的FFT算法可由下面步骤实现计算X1(k)=FFT[x1(n)];计算X2(k)=FFT[x2(n)];计算Y(k)=X1(k)X2(k);计算x1(n)*x2(n)=IFFT[Y(k)].用函数conv和FFT计算同一序列的卷积，比较其计算时间。L=5000;N=L*2-1;n=1:L;x1=0.5*n;x2=2*n;t0=clock;yc=conv(x1,x2);conv_time=etime(clock,t0)t0=clock;yf=ifft(fft(x1,N).*fft(x2,N));fft_time=etime(clock,t0)clock函数读取瞬时时钟etime(t1,t2)函数计算时刻t1,t2间所经历的时间。利用FFT完成序列的互相关函数的计算已知，x(n)={1,3,-1,1,2,3,3,1}，y(n)={2,1,-1,1,2,0,-1,3}。MATLAB实现x=[1,3,-1,1,2,3,3,1];y=[2,1,-1,1,2,0,-1,3];k=length(x);xk=fft(x,2*k);yk=fft(y,2*k);rm=real(ifft(conj(xk).*yk));rm=[rm(k+2:2*k)rm(1:k)];m=(-k+1):(k-1);stem(m,rm)四、实验报告要求简述实验目的、原理论述FFT的特点。与离散卷积结果相比较，讨论快速卷积方法的优越性。