FFT(離散傅氏變換的快速算法)

FFT(離散傅氏變換的快速算法)

本詞條是多義詞，共2個義項
更多義項 ▼ 收起列表 ▲

FFT（Fast Fourier Transformation）是離散傅氏變換（DFT）的快速算法。即為快速傅氏變換。它是根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。

基本介紹

  • 中文名：快速傅氏變換
  • 外文名：Fast Fourier Transformation
  • 表達式：FFT
  • 套用學科：數學
  • 適用領域範圍：數位訊號處理
算法,公式,源碼含義,使用方法,

算法

FFT是一種DFT的高效算法，稱為快速傅立葉變換（fast Fourier transform），它根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。FFT算法可分為按時間抽取算法和按頻率抽取算法，先簡要介紹FFT的基本原理。從DFT運算開始，說明FFT的基本原理。
DFT的運算為：
FFT
式中
FFT
由這種方法計算DFT對於
的每個K值，需要進行4N次實數相乘和（4N-2）次相加，對於N個k值，共需4N*N次實數相乘和（4N-2)*N次實數相加。改進DFT算法，減小它的運算量，利用DFT中
的周期性和對稱性，使整個DFT的計算變成一系列疊代運算，可大幅度提高運算過程和運算量，這就是FFT的基本思想。
FFT對傅氏變換的理論並沒有新的發現，但是對於在計算機系統或者說數字系統中套用離散傅立葉變換，可以說是進了一大步。
FFT算法圖(Butterfly算法)FFT算法圖(Butterfly算法)
設x(n)為N項的複數序列，由DFT變換，任一X（m）的計算都需要N次複數乘法和N-1次複數加法，而一次複數乘法等於四次實數乘法和兩次實數加法，一次複數加法等於兩次實數加法，即使把一次複數乘法和一次複數加法定義成一次“運算”（四次實數乘法和四次實數加法），那么求出N項複數序列的X（m）,即N點DFT變換大約就需要N^2次運算。當N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設N=2k,k為正整數），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以後，總的運算次數就變成N+2*（N/2)^2=N+（N^2）/2。繼續上面的例子，N=1024時，總的運算次數就變成了525312次，節省了大約50%的運算量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的DFT變換就只需要N/2log2N次的運算，N在1024點時，運算量僅有5120次，是先前的直接算法的近1/200，點數越多，運算量的節約就越大，這就是FFT的優越性。

公式

FFT
其中

源碼含義

在C環境下的源碼
源碼（1）：
註：親測，這個版本無法運行，作者刪改了重要內容請參考源碼（2）
//快速傅立葉變換// 入口參數：// l: l=0, 傅立葉變換;l=1, 逆傅立葉變換// il: il=0,不計算傅立葉變換或逆變換模和幅角；il=1,計算模和幅角// n: 輸入的點數，為偶數，一般為32，64，128，...,1024等// k: 滿足n=2^k(k>0),實質上k是n個採樣數據可以分解為偶次冪和奇次冪的次數// pr[]: l=0時，存放N點採樣數據的實部// l=1時, 存放傅立葉變換的N個實部// pi[]: l=0時，存放N點採樣數據的虛部// l=1時, 存放傅立葉變換的N個虛部//// 出口參數：// fr[]: l=0, 返回傅立葉變換的實部// l=1, 返回逆傅立葉變換的實部// fi[]: l=0, 返回傅立葉變換的虛部// l=1, 返回逆傅立葉變換的虛部// pr[]: il=1,i=0 時，返回傅立葉變換的模// il=1,i=1 時，返回逆傅立葉變換的模// pi[]: il=1,i=0 時，返回傅立葉變換的輻角// il=1,i=1 時，返回逆傅立葉變換的輻角void fft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il) {    int it,m,is,i,j,nv,l0;    double p,q,s,vr,vi,poddr,poddi;    for(it=0; it<=n-1; m=it++) {        is=0;        for(i=0; i<=k-1; i++) {            j=m/2;            is=2*is+(m-2*j);            m=j;        }        fr[it]=pr[is];        fi[it]=pi[is];    }//----------------------------    pr[0]=1.0;    pi[0]=0.0;    p=6.283185306/n;    pr[1]=cos(p);    pi[1]=-sin(p);    if (l)        pi[1]=-pi[1];    for(i=2; i<=n-1; i++) {        p=pr[i-1]*pr[1];        q=pi[i-1]*pi[1];        s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);        pr=p-q;        pi=s-p-q;    }    for(it=0; it<=n-2; it+=2) {        vr=fr[it];        vi=fi[it];        fr[it]=vr+fr[it+1];        fi[it]=vi+fi[it+1];        fr[it+1]=vr-fr[it+1];        fi[it+1]=vi-fi[it+1];    }    m=n/2;    nv=2;    for(l0=k-2; l0>=0; l0--) {        m/=2;        nv<<=1;        for(it=0; it<=(m-1)*nv; it+=nv)            for(j=0; j<=(nv/2)-1; j++) {                p=pr[m*j]*fr[it+j+nv/2];                q=pi[m*j]*fi[it+j+nv/2];                s=pr[m*j]+pi[m*j];                s*=(fr[it+j+nv/2]+fi[it+j+nv/2]);                poddr=p-q;                poddi=s-p-q;                fr[it+j+nv/2]=fr[it+j]-poddr;                fi[it+j+nv/2]=fi[it+j]-poddi;                fr[it+j]+=poddr;                fi[it+j]+=poddi;            }    }    if(l)        for(i=0; i<=n-1; fr/=n,fi[i++]/=n);    if(il)        for(i=0; i<=n-1; i++) {            pr=sqrt(fr*fr+fi*fi);            if(fabs(fr)<0.000001*fabs(fi))                pi=fi*fr>0?90.0-90.0;            else                pi=atan(fi/fr)*360.0/6.283185306;        }    return;}
源碼（2）
ps:可以運行的
// The following line must be defined before including math.h to correctly define M_PI#define _USE_MATH_DEFINES#include <math.h>#include <stdio.h>#include <stdlib.h>#define PI M_PI /* pi to machine precision, defined in math.h */#define TWOPI (2.0*PI)/*FFT/IFFT routine. (see pages 507-508 of Numerical Recipes in C)Inputs:data[] : array of complex* data points of size 2*NFFT+1.data[0] is unused,* the n'th complex number x(n), for 0 <= n <= length(x)-1, is stored as:data[2*n+1] = real(x(n))data[2*n+2] = imag(x(n))if length(Nx) < NFFT, the remainder of the array must be padded with zerosnn : FFT order NFFT. This MUST be a power of 2 and >= length(x).isign: if set to 1,computes the forward FFTif set to -1,computes Inverse FFT - in this case the output values haveto be manually normalized by multiplying with 1/NFFT.Outputs:data[] : The FFT or IFFT results are stored in data, overwriting the input.*/void four1(double data[], int nn, int isign) {    int n, mmax, m, j, istep, i;    double wtemp, wr, wpr, wpi, wi, theta;    double tempr, tempi;    n = nn << 1;    j = 1;    for (i = 1; i < n; i += 2) {        if (j > i) {            tempr = data[j];            data[j] = data[i];            data[i] = tempr;            tempr = data[j+1];            data[j+1] = data[i+1];            data[i+1] = tempr;        }        m = n >> 1;        while (m >= 2 && j > m) {            j -= m;            m >>= 1;        }        j += m;    }    mmax = 2;    while (n > mmax) {        istep = 2*mmax;        theta = TWOPI/(isign*mmax);        wtemp = sin(0.5*theta);        wpr = -2.0*wtemp*wtemp;        wpi = sin(theta);        wr = 1.0;        wi = 0.0;        for (m = 1; m < mmax; m += 2) {            for (i = m; i <= n; i += istep) {                j =i + mmax;                tempr = wr*data[j] - wi*data[j+1];                tempi = wr*data[j+1] + wi*data[j];                data[j] = data[i] - tempr;                data[j+1] = data[i+1] - tempi;                data[i] += tempr;                data[i+1] += tempi;            }            wr = (wtemp = wr)*wpr - wi*wpi + wr;            wi = wi*wpr + wtemp*wpi + wi;        }        mmax = istep;    }}
在C++環境下的源碼（經過檢驗以下代碼得到的結果有問題）
bool FFT(complex<double> * TD, complex<double> * FD, int r) {//一維快速Fourier變換。//complex<double> * TD ——指向時域數組的指針; complex<double> * FD ——指向頻域數組的指針; r ——2的冪數，即疊代次數    LONG count; // Fourier變換點數    int i,j,k; // 循環變數    int bfsize,p; // 中間變數    double angle; // 角度    complex<double> *W,*X1,*X2,*X;    count = 1 << r; // 計算Fourier變換點數為1左移r位    W = new complex<double>[count / 2];    X1 = new complex<double>[count];    X2 = new complex<double>[count]; // 分配運算所需存儲器// 計算加權係數（旋轉因子w的i次冪表）    for(i = 0; i < count / 2; i++) {        angle = -i * PI * 2 / count;        W[ i ] = complex<double> (cos(angle), sin(angle));    }// 將時域點寫入X1    memcpy(X1, TD, sizeof(complex<double>) * count);// 採用蝶形算法進行快速Fourier變換    for(k = 0; k < r; k++) {        for(j = 0; j < 1 << k; j++) {            bfsize = 1 << (r-k);            for(i = 0; i < bfsize / 2; i++) {                p = j * bfsize;                X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2] * W[i * (1<<k)];                X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2] * W[i * (1<<k)];            }        }        X = X1;        X1 = X2;        X2 = X;    }// 重新排序for(j = 0; j < count; j++) {    p = 0;    for(i = 0; i < r; i++) {        if (j&(1<<i)) {            p+=1<<(r-i-1);        }    }    FD[j]=X1[p];}// 釋放記憶體delete W;delete X1;delete X2;return true;}
在Matlab環境下的源碼
function X=myfft(x)
N=length(x);
if N==1
X=x;
return;
end
%myfft函式 用遞歸實現
t=log2(N);
t1=floor(t);
t2=ceil(t);
if t1~=t2; %若x的長度N不為2的整數次冪，則補0至最接近的2的整數次冪
x=[x zeros(1,2^t2-N)];
N=2^t2;
end
w0=exp(-1i*2*pi/N);
X=zeros(1,N);
n=1:N/2;
xo(n)=x(2*n-1);
xe(n)=x(2*n);
XO=myfft(xo); %遞歸調用
XE=myfft(xe);
for n=1:N/2
X(n)=XO(n)+XE(n)*(w0^(n-1));
X(n+N/2)=XO(n)-XE(n)*(w0^(n-1));
end

使用方法

一.調用方法

X=FFT(x)；
X=FFT(x，N)；
x=IFFT(X);
x=IFFT(X,N)
用MATLAB進行譜分析時注意：
（1）函式FFT返回值的數據結構具有對稱性。
例：
N=8;
n=0:N-1;
xn=[4 3 2 6 7 8 9 0];
Xk=fft(xn)
→
Xk =
39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i
Xk與xn的維數相同，共有8個元素。Xk的第一個數對應於直流分量，即頻率值為0。
（2）做FFT分析時，幅值大小與FFT選擇的點數有關，但不影響分析結果。在IFFT時已經做了處理。要得到真實的振幅值的大小，只要將得到的變換後結果乘以2除以N即可。
二.FFT套用舉例

例1：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變換
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*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('振幅');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;
運行結果：
fs=100Hz，Nyquist頻率為fs/2=50Hz。整個頻譜圖是以Nyquist頻率為對稱軸的。並且可以明顯識別出信號中含有兩種頻率成分：15Hz和40Hz。由此可以知道FFT變換數據的對稱性。因此用FFT對信號做譜分析，只需考察0~Nyquist頻率範圍內的幅頻特性。若沒有給出採樣頻率和採樣間隔，則分析通常對歸一化頻率0~1進行。另外，振幅的大小與所用採樣點數有關，採用128點和1024點的相同頻率的振幅是有不同的表現值，但在同一幅圖中，40Hz與15Hz振動幅值之比均為4：1，與真實振幅0.5：2是一致的。為了與真實振幅對應，需要將變換後結果乘以2除以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=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; %真實頻率
subplot(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: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頻率之前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
結論：
（1）當數據個數和FFT採用的數據個數均為32時，頻率解析度較低，但沒有由於添零而導致的其他頻率成分。
（2）由於在時間域內信號加零，致使振幅譜中出現很多其他成分，這是加零造成的。其振幅由於加了多個零而明顯減小。
（3）FFT程式將數據截斷，這時解析度較高。
（4）也是在數據的末尾補零，但由於含有信號的數據個數足夠多，FFT振幅譜也基本不受影響。
對信號進行頻譜分析時，數據樣本應有足夠的長度，一般FFT程式中所用數據點數與原含有信號數據點數相同，這樣的頻譜圖具有較高的質量，可減小因補零或截斷而產生的影響。
例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)
（1）數據點過少，幾乎無法看出有關信號頻譜的詳細信息；
（2）中間的圖是將x(n)補90個零，幅度頻譜的數據相當密，稱為高密度頻譜圖。但從圖中很難看出信號的頻譜成分。
（3）信號的有效數據很長，可以清楚地看出信號的頻率成分，一個是0.24Hz，一個是0.26Hz，稱為高解析度頻譜。
可見，採樣數據過少，運用FFT變換不能分辨出其中的頻率成分。添加零後可增加頻譜中的數據個數，譜的密度增高了，但仍不能分辨其中的頻率成分，即譜的解析度沒有提高。只有數據點數足夠多時才能分辨其中的頻率成分。

相關詞條

熱門詞條

聯絡我們