基2時抽8點FFT的matlab實現流程及FFT的內部機理


前言

本來想用verilog描述FFT算法,雖然是8點的FFT算法,但寫出來的資源用量及時延也不比調用FFT IP的好,

還是老實調IP吧,了解內部機理即可,無需重復發明輪子。

 參考

https://wenku.baidu.com/view/6f5862997c1cfad6185fa725.html

https://blog.csdn.net/shengzhadon/article/details/46737517

流程

FFT能做什么在此就不贅述了,只了解數據的運算流程。

1.FFT的基本公式:

第一眼看這個公式,肯定是腦袋瞬間宕機。

2.旋轉因子:記住旋轉因子具有可約性,對稱性,周期性。

表示方法有兩種,通過歐拉公式轉換,本質上是一致的。

Wn=exp(-j*2*k*pi/N) ,N表示FFT點數,k表示第幾個旋轉因子。

Wn = cos(2*pi*k/N) - i*sin(2*pi*k/N)

第二次腦袋瞬間宕機。

3.蝶形圖:

好在數學家為了普通人類能理解公式,繪制了帥氣的蝴蝶漫畫圖,8點FFT的如下:

不直觀,添加幾條輔助線再看:可以看到分為三級蝶形運算。

比如第一級的蝶形運算結果:x0’=x[0]+x[4]*w0,x1’=x[0]-x[4]*w0。其他節點以此類推。

注意-1的位置和旋轉因子的位置。注意數據和旋轉因子都是復數,這就是說蝶形圖中的乘法和加減都是復數運算。

而所謂代碼實現,不管什么代碼,本質上就是對各級的公式進行實現,從而得到結果。

覺得講得不清楚,那么看下圖更直觀:當然圖中少標識了-1。

4.數據輸入倒序:

 從上圖左側可以看到,序列按照了一定的規則進行了倒序,如果按照順序進行數據輸入,肯定是不正確的。十進制可能看不出來,但使其轉換為二進制表示就可以知道:

 

5.Matlab驗證算法:

到這一步,就可以把蝶形結構用matlab語言描述出來了。蝶形因子進行了2^16次放大,數據經過了兩級放大,結果需除去放大因子。

x序列為fs=500hz采樣下的125hz且有直流分量的8點采樣信號。

clc;
clear all;
close all;
%放大了2^16次的系數
w0 = 65536;
w1 = 46341 - 46341*i;
w2 = -65536*i;
w3 = -65536 - 65536*i;
% w0 = 1;
% w1 = 0.7071 - 0.7071*1i;
% w2 = -1i;
% w3 = -0.7071 - 0.7071*1i;
x = [4916,11469,4916,-1638,4916,11469,4916,-1638];
%%第一級蝶形運算,沒有放大
x1(1)=x(1)+x(5);  
x1(2)=x(1)-x(5);  
x1(3)=x(3)+x(7);  
x1(4)=x(3)-x(7);  
x1(5)=x(2)+x(6);  
x1(6)=x(2)-x(6); 
x1(7)=x(4)+x(8);  
x1(8)=x(4)-x(8);  
%%第二級蝶形運算放大了65536,因為系數放大了2^16,其他部分應當相應的放大
x2(1)=x1(1)*65536+x1(3)*65536;  
x2(2)=x1(2)*65536+x1(4)*(w2);  
x2(3)=x1(1)*65536-x1(3)*65536;  
x2(4)=x1(2)*65536-x1(4)*(w2); 
x2(5)=x1(5)*65536+x1(7)*65536;  
x2(6)=x1(6)*65536+x1(8)*(w2); 
x2(7)=x1(5)*65536-x1(7)*65536; 
x2(8)=x1(6)*65536-x1(8)*(w2);  
%%第三級蝶形運算放大了65536,因為系數放大了2^16,其他部分相應的進行放大
y(1)=x2(1)*65536+x2(5)*65536;
y(2)=x2(2)*65536+x2(6)*(w1);  
y(3)=x2(3)*65536+x2(7)*(w2); 
y(4)=x2(4)*65536+x2(8)*(w3);  
y(5)=x2(1)*65536-x2(5)*65536;  
y(6)=x2(2)*65536-x2(6)*(w1); 
y(7)=x2(3)*65536-x2(7)*(w2); 
y(8)=x2(4)*65536-x2(8)*(w3);  
% plot(abs(y/(65536^3))-abs(fft(x)))
figure;
plot(x);
title('x value');
figure;
plot(abs(y/(65536^2)));
title('旋轉因子放大取整計算結果');
figure;
plot(abs(fft(x)));
title('matlab自帶fft函數計算結果');

6.查看一下勞動成果:可以看到matlab自帶的FFT和手動蝶形算出的FFT結果是一致的。

7.轉成verilog描述,無非就是對各級的蝶形公式進行相關的實現。

注意:(1)乘法和加減法為復數運算。

(2)各級位寬需要注意,避免溢出。

看到蝶形圖及相關公式,可以看到還是有點算法復雜度的。

雖然可以手敲實現,但FPGA廠商已經提供了足夠好用的FFT IP core,資源量和計算延遲都很nice,

所以還是老實用IP吧。哈哈

8.FFT的一些性質:

(1)采樣速率和點數的關系:

頻譜分辨率△f=fs/N。點數和采樣速率共同決定了FFT的頻譜分辨率。

某一個點的頻率關系:f=k*fs/N。注意FFT計算結果的第一點為直流分量。

(2)柵欄效應及補零處理:

頻譜分辨率決定了透過柵欄窗子看真實頻譜的真實度。補零可使得離散譜外觀更加平滑,同時增長序列長度。

(3)FFT變換后的頻域模幅度對應關系:

FFT計算出的第0個點為直流分量,其模值為直流分量的N倍。其余位置求得的模值需要除以N/2,才為真實的模值。

第0點:模值/N

其他頻率點:模值/(N/2)

(4)頻域幅度及相位計算:

某點(im,re)的幅度信息為:sqrt(im^2+re^2),即實部平方加虛部平方開根號。

相位為:atan(im/re),反三角算即可,即為本頻譜時域的初相。

 

以上。

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM