本文將根據離散傅里葉變換的快速算法編寫代碼

1，首先列出DFT的變換式和快速計算公式

2，假定N=2^L，那么采用基拆分具有如下樹狀規律

3，下面在表格里面列出快速算法的一般規律：

?	左子樹	右子樹	序列分組	組內分對	W(N,k)
第1層	b[L-1]=0	b[L-1]=1	2^(L-1)	2^(1-1)	N=2^1,k=[0~2^(1-1)-1]
第2層	b[L-2]=0	b[L-2]=1	2^(L-2)	2^(2-1)	N=2^2,k=[0~2^(2-1)-1]
...	...	...	...	...	...
第i層	b[L-i]=0	b[L-i]=1	2^(L-i)	2^(i-1)	N=2^i,k=[0~2^(i-1)-1]
...	...	...	...	...	?
第L層	b[L-L]=0	b[L-L]=1	2^(L-L)	2^(L-1)	N=2^L,k=[0~2^(L-1)-1]
?	?	?	?	?	?

?

?

?

?

?

?

?

?

?

?

4，運算的時候，從序號的高位進行，向低位滑行。需要注意的是取k值，即是藍色區域的數值需要按位逆序運算。完整的代碼如下：

function X_k=fft2_dit(x_n)
? ? N1 =length(x_n);
? ? L ?=ceil(log2(N1));
? ? N ?=2^L;
? ? X_t=zeros(1,N);
? ? %如果序列不是2整數次冪，那么進行補零處理
? ? for m=1:N
? ? ? ? if(m<=N1)
? ? ? ? ? ? X_t(m)=x_n(m);
? ? ? ? else
? ? ? ? ? ? X_t(m)=0;
? ? ? ? end
? ? end
? ? for m=1:L %層級循環
? ? ? ? N ? ? =2^m;%W(N,k)中N值
? ? ? ? k_max =N/2 -1;%W(N,k)中最大k值
? ? ? ? g_max =2^(L-m);%一層有多少組獨立的蝶形圖
? ? ? ?for g=1:g_max %層內分組 ? ?
? ? ? ? ? ? for p=0:k_max %組內分對 ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? k=0; %k值逆序運算
? ? ? ? ? ? ? ? for n=1:m-1
? ? ? ? ? ? ? ? ? ? k=bitshift(k,1) + bitget(p,n); ? ? ?
? ? ? ? ? ? ? ? end ? ?
? ? ? ? ? ? ? ? num_1 ? ? =(2*p+0)*g_max+ g;
? ? ? ? ? ? ? ? num_2 ? ? =(2*p+1)*g_max+ g;?
? ? ? ? ? ? ? ? X1 ? ? ? ?=X_t(num_1) + X_t(num_2)*exp(-1j*2*pi*k/N);
? ? ? ? ? ? ? ? X2 ? ? ? ?=X_t(num_1) - X_t(num_2)*exp(-1j*2*pi*k/N);
? ? ? ? ? ? ? ? X_t(num_1)=X1;
? ? ? ? ? ? ? ? X_t(num_2)=X2; ? ? ? ? ? ? ? ??
? ? ? ? ? ? end
? ? ? ?end
? ? end
? ?
? ? X_k=zeros(1,N); ? ?
? ? for m=0:2^L-1
? ? ? ? k=0; %輸出逆序運算 ? ?
? ? ? ? for n=1:L
? ? ? ? ? ? k=bitshift(k,1) + bitget(m,n); ? ? ?
? ? ? ? end
? ? ? ? X_k(m+1)=X_t(k+1);
? ? end
end

5，下面測試按時間基2拆分的正確性

x_n =[0 1 2 3 4 5 6 7];
Xk ?=fft2_dit(x_n);%測試輸出
Xp ?=fft(x_n);%使用matlab自帶的函數
Xk =Xk.'
Xp =Xp.'

測試結果如下
Xk =

? 28.0000 ? ? ? ? ?
? -4.0000 + 9.6569i
? -4.0000 + 4.0000i
? -4.0000 + 1.6569i
? -4.0000 ? ? ? ? ?
? -4.0000 - 1.6569i
? -4.0000 - 4.0000i
? -4.0000 - 9.6569i

Xp =

? 28.0000 ? ? ? ? ?
? -4.0000 + 9.6569i
? -4.0000 + 4.0000i
? -4.0000 + 1.6569i
? -4.0000 ? ? ? ? ?
? -4.0000 - 1.6569i
? -4.0000 - 4.0000i
? -4.0000 - 9.6569i

---

6，首先列出DFT的變換式和快速計算公式

7,分析之后可以發現，每進行一層蝶形計算之后，序列對折一半。顯然對折的比特位從最高位開始，直到最低位結束，層層依次迭代。因為每進行一層運算，頻域值得序號進行一次奇偶交換，等于頻域值對折的比特位從最低位開始，直到最高位結束。需要注意的是k值不會有逆序的操作。最后得到的頻域值，進行一次逆序操作。程序如下：

function X_k=fft2_dif(x_n)
? ? N1 =length(x_n);
? ? L ?=ceil(log2(N1));
? ? N ?=2^L;
? ? X_t=zeros(1,N);
? ? %如果序列不是2整數次冪，那么進行補零處理
? ? for m=1:N
? ? ? ? if(m<=N1)
? ? ? ? ? ? X_t(m)=x_n(m);
? ? ? ? else
? ? ? ? ? ? X_t(m)=0;
? ? ? ? end
? ? end
? ? for m=1:L %層級循環
? ? ? ? N ? ? =2^(L-m+1);%W(N,k)中N值
? ? ? ? k_max =N/2-1;%W(N,k)中最大k值
? ? ? ? g_max =2^(m-1);%一層有多少組獨立的蝶形圖
? ? ? ?for g=0:g_max-1 %層內分組 ? ?
? ? ? ? ? ? for k=0:k_max %組內分對 ? ? ? ? ? ? ? ??
? ? ? ? ? ? ? ? num_1 ? ? =(2*g+0)*2^(L-m)+ k+1;
? ? ? ? ? ? ? ? num_2 ? ? =(2*g+1)*2^(L-m)+ k+1;?
? ? ? ? ? ? ? ? X1 ? ? ? ?=X_t(num_1) + X_t(num_2);
? ? ? ? ? ? ? ? X2 ? ? ? ?=(X_t(num_1) - X_t(num_2))*exp(-1j*2*pi*k/N);
? ? ? ? ? ? ? ? X_t(num_1)=X1;
? ? ? ? ? ? ? ? X_t(num_2)=X2; ? ? ? ? ? ? ? ??
? ? ? ? ? ? end
? ? ? ?end
? ? end
? ?
? ? X_k=zeros(1,N); ? ?
? ? for m=0:2^L-1
? ? ? ? k=0; %輸出逆序運算 ? ?
? ? ? ? for n=1:L
? ? ? ? ? ? k=bitshift(k,1) + bitget(m,n); ? ? ?
? ? ? ? end
? ? ? ? X_k(m+1)=X_t(k+1);
? ? end
end

8，下面測試按頻率基2拆分的正確性

x_n =[0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1];
Xk ?=fft2_dit(x_n);%測試輸出
Xp ?=fft(x_n);%使用matlab自帶的函數
Xk ?=Xk.'
Xp ?=Xp.'

測試結果如下

Xk =

? 64.0000 ? ? ? ? ?
?-26.2741 - 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -3.2398 - 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -1.4465 + 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -1.0396 + 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -1.0396 + 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -1.4465 + 0.0000i
? ? ? ? 0 ? ? ? ? ?
? -3.2398 + 0.0000i
? ? ? ? 0 ? ? ? ? ?
?-26.2741 - 0.0000i

Xp =

? ?64.0000
? -26.2741
? ? ? ? ?0
? ?-3.2398
? ? ? ? ?0
? ?-1.4465
? ? ? ? ?0
? ?-1.0396
? ? ? ? ?0
? ?-1.0396
? ? ? ? ?0
? ?-1.4465
? ? ? ? ?0
? ?-3.2398
? ? ? ? ?0
? -26.2741

---

自此，快速傅里葉算法研究到此結束，如果對你有所幫助，請點一個贊，謝謝。

總結

以上是生活随笔為你收集整理的快速傅里叶变换FFT的代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。