包含java實現fft變換的詞條

本文目錄一覽：

• 1、求快速傅里葉變換的演算法實現，C/C++/JAVA都行，要求適用於所有的整數，謝謝！
• 2、z7 7020中如何實現fft
• 3、java繪製wav波形
• 4、一維複數序列的快速傅里葉變換（FFT）
• 5、FFT的公式是什麼和演算法是怎樣實現

求快速傅里葉變換的演算法實現，C/C++/JAVA都行，要求適用於所有的整數，謝謝！

#includegraphics.h

#includestdio.h

#includemath.h

#includeconio.h

#includedos.h

float ar[2048],ai[2048];

float a[4100];

void testdata()

{

int i;

for(i=0;i512;i++)

{

ar[i]=50*cos(i*3.14159/32)+60*cos(i*3.14159/16)+53*cos(i*3.14159/64)+24*cos(i*3.14159/8)+10*cos(i*3.14159/128);

ai[i]=0;

}

}

void fft(int nn,int t)

{

int n1,n2,i,j,k,l,m,s,l1;

float t1,t2,x,y;

float w1,w2,u1,u2,z;

float pi=3.14159;

switch(nn)

{

case 2048: s=11; break;

case 1024: s=10; break;

case 512: s=9; break;

case 256: s=8; break;

}

n1=nn/2; n2=nn-1;

j=1;

for(i=1;i=nn;i++)

{

a[2*i]=ar[i-1];

a[2*i+1]=ai[i-1];

}

for(l=1;ln2;l++)

{

if(lj)

{

t1=a[2*j];

t2=a[2*j+1];

a[2*j]=a[2*l];

a[2*j+1]=a[2*l+1];

a[2*l]=t1;

a[2*l+1]=t2;

}

k=n1;

while (kj)

{

j=j-k;

k=k/2;

}

j=j+k;

}

for(i=1;i=s;i++)

{

u1=1;

u2=0;

m=(1i);

k=m1;

w1=cos(pi/k);

w2=t*sin(pi/k);

for(j=1;j=k;j++)

{

for(l=j;lnn;l=l+m)

{

l1=l+k;

t1=a[2*l1]*u1-a[2*l1+1]*u2;

t2=a[2*l1]*u2+a[2*l1+1]*u1;

a[2*l1]=a[2*l]-t1;

a[2*l1+1]=a[2*l+1]-t2;

a[2*l]=a[2*l]+t1;

a[2*l+1]=a[2*l+1]+t2;

}

z=u1*w1-u2*w2;

u2=u1*w2+u2*w1;

u1=z;

}

}

for(i=1;i=nn/2;i++)

{

ar[i]=a[2*i+2]/nn;

ai[i]=-a[2*i+3]/nn;

a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]);

}

}

int main ()

{

int i;

int gdriver=DETECT,gmode;

initgraph(gdriver,gmode,””);

setfillstyle(SOLID_FILL,WHITE);

bar(0,0,639,479);

//模擬測試數據

testdata();

//波形顯示

setcolor(BLACK);

line(10,119,550,119); // x軸

line(10,10,10,239); // y軸

setcolor(BLUE);

for(i=0;i511;i++)

line(i+10,119-ar[i],i+10+1,119-ar[i+1]);

//FFT

fft(512,-1);

//頻譜顯示

setcolor(BLACK);

line(10,459,550,459); // x軸

line(10,259,10,459); // y軸

setcolor(RED);

for(i=0;i256;i++)

line(2*i+10,459-a[i],2*i+10,459);

getch();

closegraph();

return 0;

}

z7 7020中如何實現fft

自己重新編程。

_ft是內建函數，不是matlab寫的，看不到源代碼的。下面是我寫的一個fft,可以用functionxn=myfft(x)N=length(x);M=log2(N);xtmp=zeros(1,N);value=zeros(1,M);for i=0:N-1repr=i;fort=1:1:Mrepr=bitshift(i,1-t);value(t)=bitand(repr,1);endpos=0;for k=1:1:Mpos=pos+value(k)*2^(M-k);endxtmp(pos+1)=x(i+1);endfor i=1:deepth=2^(i-1);width=2^(M-i);for t=1:2^i:Nfor k=1:deepthtmp=xtmp(t+k-1);wn=width*(k-1);xtmp(t+k-1)=tmp+exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1);xtmp(t+k+deepth-1)=tmp-exp(-j*2*pi*wn/N)*xtmp(t+k+deepth-1);endendendxn=xtmp;

__FT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下： 先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。

java繪製wav波形

本文講述如何實時錄音，以及將錄音波形頻譜實時顯示的方法。Windows提供了一個多媒體控制介面(MCI)，用它可以錄音，很方便。但是這種方法不能實時給出錄音的原始數據，因此要顯示波形和頻譜都是不可能實現的。要達到實時的效果，就要使用Windows提供的另一套函數，即低級音頻函數。其中和錄音有關的是以waveIn開頭的一組函數。

低級音頻函數的使用比較繁瑣，大致要有以下幾個步驟。

用waveInOpen打開設備，並設置回調。因為要保證實時性，所以不能用查詢的方式，而必須設置回調。

為設備分配足夠的內存做緩衝區，動態分配或靜態數組都可以。為了保證實時性，程序用了雙緩衝技術，在處理一個緩衝區數據的同時另一個緩衝區用於錄音。為了便於說明寫成Buffer1、Buffer2。

將Buffer1關聯到設備上去，waveInPrepareBuffer、waveInAddBuffer。

開始錄音，waveInStart

當驅動程序填滿這個緩衝區(Buffer1)時就會產生回調(消息為WIM_DATA)，這時立刻將Buffer2關聯到設備上繼續錄音，然後處理Buffer1，當驅動程序填滿Buffer2時又會產生回調，這是再將Buffer1關聯到設備上，而去處理Buffer2，如此反覆就使得錄音能夠實時的進行下去。

停止錄音，waveInStop

關閉設備，waveInClose

上面就是實時錄音的一個基本流程，產生回調時，其lParam就是一個指向WAVEHDR結構的指針，通過這個指針就可以得到原始數據(注意了：我這個程序只是單聲道8位的PCM格式，因此正好一個位元組就是一個採樣點的值，對於其他格式的我就不知道了)。有了原始數據就可以很容易的畫出波形了，而頻譜也就是多一次FFT變換而已，沒什麼太難的。

這裡是我用LCC寫的例子，非常簡單，而且相信大家都可以很容易的移植到BCB上來。另外，我偷了一下懶，FFT用了我以前封裝到DLL中的函數(現在已經忘了怎麼寫了^_^)，可以在此下載 (不要覺得奇怪，就是這個時鐘，將其中的DEMO.DLL拿來用就是了)。當然你也可以自己寫一個FFT放到程序里去。

-*-*-PATCH-*-*-2002-06-11

上次那個程序還不好，原因是雙緩衝還不能滿足實時的要求，因此這次用了8緩衝循環技術。新的程序也是用LCC寫的。

開始時依次將0－6號緩衝區關聯到設備上，而7號緩衝區不關聯，然後開始錄音，這是驅動程序開始向0號緩衝區輸出數據，當0號緩衝區數據滿時產生回調，這時驅動程序會自動向1號緩衝區輸出數據，但這時應將7號緩衝區關聯到設備上去，然後處理0號緩衝區的數據。然後當1號緩衝區數據滿時又會產生回調，這時驅動程序會自動向2號緩衝區輸出數據，而將0號緩衝區關聯到設備上去。這樣如此反覆，始終保持1個緩衝區用於驅動輸出數據，1個緩衝區用於處理數據，6個緩衝區處於等待狀態。

而如果只有兩個緩衝區，就只能一個用於驅動輸出，一個用於處理數據。這樣在產生回調的時候就會出現沒有處於等待狀態的緩衝區的情況，這樣驅動程序就不能自動向緩衝區輸出數據。而要等waveInPrepareBuffer、waveInAddBuffer調用之後，這樣就會出現一個小小的間隙，而導致錄音不連續。

一維複數序列的快速傅里葉變換（FFT）

設x（N）為N點有限長離散序列，代入式（8－3）、式（8－4），並令 其傅里葉變換（DFT）為

地球物理數據處理基礎

反變換（IDFT）為

地球物理數據處理基礎

兩者的差異只在於W的指數符號不同，以及差一個常數1/N，因此下面我們只討論DFT正變換式（8－5）的運算量，其反變換式（8－6）的運算是完全相同的。

一般來說，W是複數，因此，X（j）也是複數，對於式（8－5）的傅里葉變換（DFT），計算一個X（j）值需要N次複數乘法和N－1次複數加法。而X（j）一共有N個值（j＝0，1，…，N－1），所以完成整個DFT運算總共需要N2次複數乘法和N（N－1）次複數加法。

直接計算DFT，乘法次數和加法次數都是與N2成正比的，當N很大時，運算量會很大，例如，當N＝8時，DFT需64次複數乘法；而當N＝1024時，DFT所需乘法為1048576次，即一百多萬次的複數乘法運算，對運算速度要求高。所以需要改進DFT的計算方法，以減少運算次數。

分析Wjk，表面上有N2個數值，由於其周期性，實際上僅有N個不同的值W0，W1，…，WN－1。對於N＝2m時，由於其對稱性，只有N/2個不同的值W0，W1，…，

地球物理數據處理基礎

因此可以把長序列的DFT分解為短序列DFT，而前面已經分析DFT與N2成正比，所以N越小越有利。同時，利用ab＋ac＝a（b＋c）結合律法則，可以將同一個Wr對應的係數x（k）相加後再乘以Wr，就能大大減少運算次數。這就是快速傅里葉變換（FFT）的演算法思路。

下面，我們來分析N＝2m情況下的FFT演算法。

1.N＝4的FFT演算法

對於m＝2，N＝4，式（8－5）傅里葉變換為

地球物理數據處理基礎

將式（8－7）寫成矩陣形式

地球物理數據處理基礎

為了便於分析，將上式中的j，k寫成二進位形式，即

地球物理數據處理基礎

代入式（8－7），得

地球物理數據處理基礎

分析Wjk的周期性來減少乘法次數

地球物理數據處理基礎

則 代回式（8－9），整理得

地球物理數據處理基礎

上式可分層計算，先計算內層，再計算外層時就利用內層計算的結果，可避免重複計算。寫成分層形式

地球物理數據處理基礎

則X（j1 j0）＝X2（j1 j0）。

上式表明對於N＝4的FFT，利用Wr的周期關係可分為m＝2步計算。實際上，利用Wr的對稱性，仍可以對式（8－11）進行簡化計算。考慮到

地球物理數據處理基礎

式（8－11）可以簡化為

地球物理數據處理基礎

令j＝j0；k＝k0，並把上式表示為十進位，得

地球物理數據處理基礎

可以看到，完成上式N＝4的FFT計算（表8－1）需要N·（m－1）/2＝2次複數乘法和N·m＝8次複數加法，比N＝4的DFT演算法的N2＝16次複數乘法和N·（N－1）＝12次複數加法要少得多。

表8－1 N＝4的FFT演算法計算過程

註：W0＝1；W1＝－i。

［例1］求N＝4樣本序列1，3，3，1的頻譜（表8－2）。

表8－2 N＝4樣本序列

2.N＝8的FFT演算法

類似N＝4的情況，用二進位形式表示，有

地球物理數據處理基礎

寫成分層計算的形式：

地球物理數據處理基礎

則X（j2 j1 j0）＝X3（j2 j1 j0）。

對式（8－14）的X1（k1 k0 j0）進行展開，有

地球物理數據處理基礎

還原成十進位，並令k＝2k1＋k0，即k＝0，1，2，3，有

地球物理數據處理基礎

用類似的方法對式（8－14）的X2（k0 j1 j0），X3（j2 j1 j0）進行展開，整理得

地球物理數據處理基礎

用式（8－16）、式（8－17）逐次計算到X3（j）＝X（j）（j＝0，1，…，7），即完成N＝23＝8的FFT計算，其詳細過程見表8－3。

表8－3 N＝8的FFT演算法計算過程

註：對於正變換 對於反變換 所

［例2］求N＝8樣本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的頻譜。

表8－4 N＝8樣本序列

3.任意N＝2m的FFT演算法

列出N＝4，N＝8的FFT計算公式，進行對比

地球物理數據處理基礎

觀察式（8－18）、式（8－19），不難看出，遵循如下規律：

（1）等式左邊的下標由1遞增到m，可用q＝1，2，…，m代替，則等式右邊為q－1；

（2）k的上限為奇數且隨q的增大而減小，至q＝m時為0，所以其取值範圍為k＝0，1，2，…，（2m－q－1）；

（3）j的上限為奇數且隨q的增大而增大，且q＝1時為0，其取值範圍為j＝0，1，2，…，（2q－1－1）；

（4）k的係數，在等式左邊為2q，等式右邊為2q－1（包括W的冪指數）；

（5）等式左邊序號中的常數是2的乘方形式，且冪指數比下標q小1，即2q－1；等式右邊m對式子序號中的常數都是定值2m－1。

歸納上述規則，寫出對於任意正整數m，N＝2m的FFT演算法如下：

由X0（p）＝x（p）（p＝0，1，…，N－1）開始：

（1）對q＝1，2，…，m，執行（2）～（3）步；

（2）對k＝0，1，2，…，（2m－q－1）及j＝0，1，2，…，（2q－1－1），執行

地球物理數據處理基礎

（3）j，k循環結束；

（4）q循環結束；由Xm（p）（p＝0，1，…，N－1）輸出原始序列x（p）的頻譜X（p）。

在計算機上很容易實現上述FFT演算法程序，僅需要三個複數數組，編程步驟如下：

（1）設置複數數組X1（N－1），X2（N－1）和 （數組下界都從0開始）；

（2）把樣本序列x賦給X1，即X1（k）＝x（k）（k＝0，1，…，N－1）；

（3）計算W，即正變換 反變換

（4）q＝1，2，…，m，若q為偶數，執行（6），否則執行第（5）步；

（5）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循環，作

X2（2qk＋j）＝X1（2q－1k＋j）＋X1（2q－1k＋j＋2m－1）

X2（2qk＋j＋2q－1）＝［X1（2q－1k＋j）－X1（2q－1k＋j＋2m－1）］W（2q－1k）

至k，j循環結束；

（6）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循環，作

X1（2qk＋j）＝X2（2q－1k＋j）＋X2（2q－1k＋j＋2m－1）

X1（2qk＋j＋2q－1）＝［X2（2q－1k＋j）－X2（2q－1k＋j＋2m－1）］W（2q－1k）

至k，j循環結束；

（7）q循環結束，若m為偶數，輸出X1（j），否則輸出X2（j）（j＝0，1，…，N－1），即為所求。

FFT的公式是什麼和演算法是怎樣實現

二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：

先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：

for (int i=0; iM; i++)

FFT_1D(ROW[i]，N);

for (int j=0; jN; j++)

FFT_1D(COL[j]，M);

其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。

所以，關鍵是一維FFT演算法的實現。下面討論一維FFT的演算法原理。

【1D-FFT的演算法實現】

設序列h(n)長度為N，將其按下標的奇偶性分成兩組，即he和ho序列，它們的長度都是N/2。這樣，可以將h(n)的FFT計算公式改寫如下 ：

(A)

由於

所以，(A)式可以改寫成下面的形式：

按照FFT的定義，上面的式子實際上是：

其中，k的取值範圍是 0~N-1。

我們注意到He(k)和Ho(k)是N/2點的DFT，其周期是N/2。因此，H(k)DFT的前N/2點和後N/2點都可以用He(k)和Ho(k)來表示