FFT信號流圖：

程序實現是這樣：

程序流程如下圖：

首先進行位逆轉，其實很簡單，就是把二進制的位逆轉過來：

Matlab的位逆轉程序：

function a=bitreverse(Nbit, num)

%Nbit = 4;

%num = 8;

a = 0;

b = bitshift(1,Nbit-1);

for i = 1:Nbit;

if((bitand(num,1)) == 1)

a = bitor(a,b);

end

num = bitshift(num,-1);

b = bitshift(b,-1);

end;

說明：Nbit是逆轉位是幾位，num是逆轉的數即變量。

三個循環，第一個循環是進行N階的FFT運算

第二個循環其實就是，每一階FFT的時候，有多少組DFT對象，拿8點來說，第一階的時候，有4組DFT對象，到了第二階，就有2組，到了第三，就是最后一階，只有一組。

第三個循環，其實是在每一組DFT里邊，執行多少次蝶形運算！8點DIT FFT來說，第一階每組有一個蝶形，第二階每組有2個，第三階每組有4個蝶形。所以很容易得到三者的關系

i , j, k

三者，反別表示三層循環，然后得出循環次數的關系

stages = log2(PointNum)

i

從 0到stages – 1 ！

j

從 0 到其實就是

k

從0 到

旋轉因子W的選擇：

因為根據8點DIT-FFT圖，從第一階到最后一階，可以總結出一個規律：

都是 N是每組蝶形數據個個數，比如第一階每組有2個元素，N就是2，第二階每組4個元素，N就是4等。然后x往往都是從0開始到N/2 – 1;

根據旋轉因子的性質，其實可以有每階段每組都是：

蝶形運算設計：

根據信號流圖，得出以下算式：

完成了蝶形運算!

全部的matlab程序有：

PointNum = 512;

PointBitNum = 9;

fs = 1024*2;

t = 0:1:PointNum - 1;

%for u = 1:1:PointNum;

sampletab = cos(2*pi*543*t/fs) +

cos(2*pi*100*t/fs) + 0.2 + cos(2*pi*857*t/fs) +

cos(2*pi*222*t/fs);

%end

zeros(1,PointNum);

sampletab1 = sampletab;

index = 0;

for i = 1:PointNum

k =

i - 1

index = bitreverse(PointBitNum,i - 1)

sampletab(i) = sampletab1(index + 1);

end

%sampletab1

%sampletab

REX = sampletab;

IMX = zeros(1,PointNum);

i = 0; �T Loop for Log2N stages

j = 0; �T loop for leach sub-DFT

k = 0; �T Loop for each butterfly

stages = log2(PointNum);

for i = 0 : stages - 1

lenNum = 0;

for j = 0 : 2^(stages - (i + 1)) -

1

for k = 0 : 2^i - 1

R1 = REX(lenNum + 2^i + 1) * cos(2*pi*k*2^(stages - (i +

1))/PointNum);

R2 = IMX(lenNum + 2^i + 1) * sin(2*pi*k*2^(stages - (i +

1))/PointNum);

T1 = REX(lenNum + 2^i + 1) * sin(2*pi*k*2^(stages - (i +

1))/PointNum);

T2 = IMX(lenNum + 2^i + 1) * cos(2*pi*k*2^(stages - (i +

1))/PointNum);

REX(lenNum + 2^i + 1) = REX(lenNum + 1) - R1 - R2;

IMX(lenNum + 2^i + 1) = IMX(lenNum + 1) + T1 - T2;

REX(lenNum + 1) = REX(lenNum + 1) + R1 + R2;

IMX(lenNum + 1) = IMX(lenNum + 1) - T1 + T2 ;

lenNum = lenNum + 1;

endNum

= lenNum + 2^i;

end

lenNum = endNum;

end

end

subplot(3,1,1);

fft(sampletab1, PointNum);

x1 = abs(fft(sampletab1,

PointNum));

plot([0 : PointNum/2 - 1],

x1(1:PointNum/2));

grid on

subplot(3,1,2);

% [REX IMX]

am = sqrt(abs(REX.*REX) +

abs(IMX.*IMX));

plot(0:1:PointNum/2 - 1,

am(1:PointNum/2));

grid on

subplot(3,1,3);

plot(t, sampletab);

grid on

我還做了與MATLAB原來帶有的FFT做比較：

畫出的圖如下：

第一個是MATLAB自帶的FFT函數頻譜圖

第二個是我自己設計的FFT頻譜圖

第三個是信號的時域波形

思想已經有了，我以前也改過人家的FFT的C程序但是不是很理解，打算有機會用C語言實現定點FFT，因為在嵌入式上多數用定點FFT，相應的C++版本應該也會寫。

下面是網上的一些設計FFT的資料：

N點基-2 FFT算法的實現方法

從圖4我們可以總結出對于點數為N=2^L的DFT快速計算方法的流程：

1.對于輸入數據序列進行倒位序變換。

該變換的目的是使輸出能夠得到X(0)~X(N-1)的順序序列，同樣以8點DFT為例，該變換將順序輸入序列x(0)~x(7)變為如圖4的x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其實現方法是：假設順序輸入序列一次村在A(0)~A(N-1)的數組元素中，首先我們將數組下標進行二進制化(例：對于點數為8的序列只需要LOG2(8)

=

3位二進制序列表示，序號6就表示為110)。二進制化以后就是將二進制序列進行倒位，倒位的過程就是將原序列從右到左書寫一次構成新的序列，例如序號為6的二進制表示為110，倒位后變為了011，即使十進制的3。第三步就是將倒位前和倒位后的序號對應的數據互換。依然以序號6為例，其互換過程如下：

temp = A(6);?A(6) = A(3);?A(3) = A(6);

實際上考慮到執行效率，如果對于每一次輸入的數據都需要這個處理過程是非常浪費時間的。我們可以采用指向指針的指針來實現該過程，或者是采用指針數組來實現該過程。

2.蝶形運算的循環結構。

從圖4中我們可以看到對于點數為N =

2^L的fft運算，可以分解為L階蝶形圖級聯，每一階蝶形圖內又分為M個蝶形組，每個蝶形組內包含K個蝶形。根據這一點我們就可以構造三階循環來實現蝶形運算。編程過程需要注意旋轉因子與蝶形階數和蝶形分組內的蝶形個數存在關聯。

3.浮點到定點轉換需要注意的關鍵問題?上邊的分析都是基于浮點運算來得到的結論，事實上大多數嵌入式系統對浮點運算支持甚微，因此在嵌入式系統中進行離散傅里葉變換一般都應該采用定點方式。對于簡單的DFT運算從浮點到定點顯得非常容易。根據式(1)，假設輸入x(n)是經過AD采樣的數字序列，AD位數為12位，則輸入信號范圍為0~4096。為了進行定點運算我們將旋轉因子實部虛部同時擴大2^12倍，取整數部分代表旋轉因子。之后，我們可以按照(1)式計算，得到的結果與原結果成比例關系，新的結果比原結果的2^12倍。但是，對于使用蝶形運算的fft我們不能采用這種簡單的放大旋轉因子轉為整數計算的方式。因為fft是一個非對稱迭代過程，假設我們對旋轉因子進行了放大，根據蝶形流圖我們可以發現其最終的結果是，不同的輸入被放大了不同的倍數，對于第一個輸入x(0)永遠也不會放大。舉一個更加形象的例子，還是以圖4為例。從圖中可以看出右側的X(0)可以直接用下式表示：?

從上式我們可以看到不同輸入項所乘的旋轉因子個數(注意這里是個數，就算是wn^0，也被考慮進去了，因為在沒有放大時wn^0等于1，放大后所有旋轉因子指數模均不為1，因此需要考慮)。這就導致輸入不平衡，運算結果不正確。經查閱相關資料，比較妥善的做法是，首先對所有旋轉因子都放大2^Q倍，Q必須要大于等于L，以保證不同旋轉因子的差異化。旋轉因子放大，為了保證其模為1，在每一次蝶形運算的乘積運算中我們需要將結果右移Q位來抵消這個放大，從而得到正確的結果。之所以采用放大倍數必須是2的整數次冪的原因也在于此，我們之后可以通過簡單的右移位運算將之前的放大抵消，而右移位又代替了除法運算，大大節省了時間。?4.計算過程中的溢出問題?最后需要注意的一個問題就是計算過程中的溢出問題。在實際應用中，AD雖然有12位的位寬，但是采樣得到的信號可能較小，例如可能在0~8之間波動，也就是說實際可能只有3位的情況。這種情況下為了在計算過程中不丟失信息，一般都需要先將輸入數據左移P位進行放大處理，數據放大可能會導致溢出，從而使計算錯誤，而溢出的極限情況是這樣：假設我們數據位寬為D位(不包括符號位)，AD采樣位數B位，數字放大倍數P位，旋轉因此放大倍數Q位，FFT級聯運算帶來的最大累加倍數L位。我們得到:?

假設AD位寬12，數據位寬32，符號位1位，因此有效位寬31位，采樣點數N，那么我們可以得到log2(N)+P+Q<=19,假設點數128,又Q>=L可以得到放大倍數P<=5。

總結

以上是生活随笔為你收集整理的matlab 倒位序fft程序,[转载]MATLAB的一个FFT程序的全部內容，希望文章能夠幫你解決所遇到的問題。

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