遞迴版的快速傅利葉轉換 FFT ，C語言實現

看了好多個版本的快速傅利葉轉換 (fft)，看懂了理論部分，但是卻無法把理論聯絡到**。也許是我的水平太差了吧。網上看到的版本中大量使用位操作符《和》。開始的時候還真的反應不過來。這裡，我不講fft的理論，只是實現一下遞迴版本的fft實現。有網友說fft是玄科學，我深表贊同。相同的**，在a電腦裡執行結果是正確的，用b電腦就錯誤。找了好久才發現只需要把兩個陣列申請很大的空間即可，但這是為什麼呢？反正沒人告訴我原因。下面是具體**：

#include
#include
#include
double pi =
3.1415926
;//定義複數 
struct complex
;//網上大部分都是採用c++實現fft，使用了操作符過載。
//但是其它語言沒有操作符過載，所以我們還是用函式的方法實現複數的加法、減法和乘法
struct complex add
(struct complex a,
struct complex b)
//複數的加法 
//複數的減法，不能用minus，好像是與關鍵字衝突 
struct complex myminus
(struct complex a,
struct complex b)
struct complex multiply
(struct complex a,
struct complex b)
//複數的乘法
void
fft(
int n, complex * arr,
int type)
int m = n /2;
//由於保證了n是足夠大的偶數，所以m一定是n的一半 
struct complex *arr1;
//儲存arr的偶數部分 
arr1 =
(struct complex *)(
malloc
(m *
sizeof
(struct complex)))
;//動態申請記憶體 
struct complex *arr2;
//儲存arr的奇數部分
arr2 =
(struct complex *)(
malloc
(m *
sizeof
(struct complex)))
;//動態申請記憶體 
for(
int i =
0; i < m; i++
)fft
(m, arr1, type)
;//對偶數部分進行快速傅利葉變換 
fft(m, arr2, type)
;//對奇數部分進行快速傅利葉變換
struct complex wn;
//把360度等分成n份，這時的複數。 
wn.x =
cos(
2.0* pi / n)
;//複數的x 
wn.y = type *
sin(
2.0* pi / n)
;//複數的y，注意逆運算的時候，type為-1 
struct complex w;
//x = 1，y = 0的複數 
w.x =1;
w.y =0;
//計算arr中的每乙個值 
for(
int i =
0; i < m; i++)}
intmain()
for(
int i =
0; i < m; i++
)//陣列b的初始化
int lim =1;
//lim為1, 2, 4, 8, 16, 32....除了1之外，其餘都是偶數
//lim為1時，fft直接就遞迴返回了，所以實際上lim都是偶數 
while
(lim < n + m -1)
//即兩函式的傅利葉變換的乘積等於它們卷積後的傅利葉變換，
//能使傅利葉分析中許多問題的處理得到簡化。
//f(g(x) * f(x)) = f(g(x)) * f(f(x))
//其中f表示的是傅利葉變換
fft(lim, a,1)
;//快速傅利葉轉換 
fft(lim, b,1)
;//快速傅利葉轉換
for(
int i =
0; i < lim; i++
)fft
(lim, c,-1
);//對c進行快速傅利葉轉換的逆轉換 
//取實數四捨五入，此時虛數部分應當為0或由於浮點誤差接近0
//理論上計算結果應該是虛部為0，但實際上虛部可能有一點誤差
//為了抵消這個誤差，用下面的方法進行處理
//最終顯示c的內容 
for(
int i =
0; i < n + m -
1; i++
)return0;
}

執行結果如下：

如果main函式中，struct complex * a、b、c申請的空間比較小 (100或更小)，則會出錯，不知道為什麼。

下次再寫迭代版本的fft吧。

遞迴版的快速傅利葉轉換 FFT ，C語言實現

迭代版的快速傅利葉轉換 FFT ，C語言實現

快速傅利葉

快速傅利葉（FFT）

遞迴版的快速傅利葉轉換 FFT ，C語言實現

迭代版的快速傅利葉轉換 FFT ，C語言實現

快速傅利葉

快速傅利葉（FFT）

相關推薦