數(shù)字信號(hào)處理 - 基4FFT算法編程

實(shí)驗(yàn)三基4－FFT算法編程實(shí)驗(yàn)?zāi)康暮鸵驠FT是快速計(jì)算DFT的一類(lèi)算法的總稱(chēng)。通過(guò)序列分解，用短序列的DFT代替長(zhǎng)序列的DFT，使得計(jì)算量大大下降?；?-FFT是混合基FFT的一個(gè)特例。通過(guò)編寫(xiě)基4-FFT算法程序，加深對(duì)FFT思路、算法結(jié)構(gòu)的理解。實(shí)驗(yàn)內(nèi)容與步驟編寫(xiě)16點(diǎn)基4-FFT算法的MATLAB程序（studentname.m文件）。產(chǎn)生16點(diǎn)輸入序列x，用自己的學(xué)好作為前10點(diǎn)的抽樣值，后面補(bǔ)6個(gè)零值抽樣。算出16點(diǎn)頻譜序列X，用stem(X)顯示頻譜圖形?；?－FFT算法簡(jiǎn)述

基4與基2的算法思路是一樣的，區(qū)別僅僅在于每一組的元素個(gè)數(shù)不同，但是在算法流圖方面差別還是比較大的。假設(shè)十六點(diǎn)序列為：x(k)

0

<

k

<15

基4－FFT算法簡(jiǎn)述將序列進(jìn)行如下分組：x(0)=

x(4n)(1)xx(2)(3)=

x(4n+2)=

x(4n+

3)其中n取值從零到4X(k)=∑x(4l)W+∑x(4l

+1)W16(4l+1)k

+∑x(4l

+2)W+∑x(4l

+

3)W16(4l+3)k

=

基4－FFT算法簡(jiǎn)述原序列的DFT可以變成：

3l=0

4lk16

3l=0

3l=0

(4l+2)k16

3l=0k

2

3

X(0)((k))4

+W16X(1)((k))4

+W16kX(2)((k))4

+W16kX(3)((k))4可表示為為四個(gè)子序列的4點(diǎn)DFT的和的形式基4－FFT算法簡(jiǎn)述上式可以寫(xiě)為：如果令0≤

k

≤

3k

2

3kkk

X(k)=

X(0)(k)+W16X(1)(k)+W16kX(2)(k)+W16kX(3)(k)

X(k

+

4)=

X(0)(k)?

jW16X(1)(k)?W162kX(2)(k)+

jW163kX(3)(k)

X(k

+8)=

X(0)(k)?W16X(1)(k)+W162kX(2)(k)?W163kX(3)(k)

X(k

+12)=

X(0)(k)+

jW16X(1)(k)?W162kX(2)(k)?

jW163kX(3)(k)這樣就用四個(gè)子序列的DFT完全表示了16點(diǎn)序列的16點(diǎn)DFT流圖結(jié)構(gòu)簡(jiǎn)述

時(shí)域基2-FFT可以把4點(diǎn)DFT轉(zhuǎn)化為兩步蝶形運(yùn)算，第一步為求兩個(gè)子序列的兩點(diǎn)DFT，第二步蝶形運(yùn)算根據(jù)子序列的DFT求出原序列的4點(diǎn)DFT。

16點(diǎn)的基4－FFT應(yīng)該是類(lèi)似的，也可以分成兩步蝶形運(yùn)算，只是每一個(gè)蝶形運(yùn)算的運(yùn)算對(duì)象變成了四個(gè)。

第一步應(yīng)該是求四個(gè)子序列的4點(diǎn)DFT，第二步根據(jù)第一步的結(jié)果得到16點(diǎn)序列的DFT流圖如下:(因?yàn)樘^(guò)于繁瑣，圖中的箭頭上省略了一些系數(shù))第一步為求四個(gè)子序列的四點(diǎn)DFT第二步求原序列的16點(diǎn)DFTMatlab

b程序此部分程序?yàn)樾蛄械亩x，分組以及第一次蝶形運(yùn)算的代碼Matlab

b程序此部分為第二次蝶形運(yùn)算和畫(huà)圖的代碼實(shí)驗(yàn)結(jié)果用自己的學(xué)號(hào)加6個(gè)零構(gòu)成的序列的圖形如下：

實(shí)驗(yàn)結(jié)果

所得序列的值如下：X[k]=[

19,

-7.39-0.12i,

12.07-5.66i,

-3.06+5.13i,5-8i,

3.06+5.96i,

-2.07-5.66i,7.39+4.71i,-1,7.39-4.71i,-2.07+5.66i,3.06-06-55.96i,5+8i,

-3

3.06

5.13i,12.07+5.66i,

-7

7.39+0.12i]。實(shí)驗(yàn)結(jié)果序列實(shí)部圖形：實(shí)驗(yàn)結(jié)果序列虛部圖形：實(shí)驗(yàn)結(jié)果

為了檢驗(yàn)結(jié)果的正確性，將X[k]與用matlab中的fft函數(shù)算出的結(jié)果相減在matlab中輸入：X-fft(x)得到的結(jié)果如下所示：說(shuō)明結(jié)果沒(méi)有問(wèn)題總結(jié)本實(shí)驗(yàn)的要點(diǎn)在于：對(duì)基41.

基4-FFT

FFT算法的理解2.流圖的繪制

其實(shí)畫(huà)圖很容易3.Matlab程序的編寫(xiě)

使用循環(huán)語(yǔ)句可以大大減少工作量2一、實(shí)驗(yàn)?zāi)康暮鸵驠FT是快速計(jì)算DFT的一類(lèi)算法的總稱(chēng)。通過(guò)序列分解，用短序列的DFT代替長(zhǎng)序列的DFT，使得計(jì)算量大大下降?；?-FFT是混合基FFT的一個(gè)特例。通過(guò)編寫(xiě)基4-FFT算法程序，加深對(duì)FFT思路、算法結(jié)構(gòu)的理解。∑∑∑∑∑x

(l)W+W∑x

(l)W+W∑x

(l)W+W∑x(3)

)(l)W4lk

m?13我們可以得到

X(k)的計(jì)算式為：4NNN4m?1?1

l=04m?1?1

l=04m?1?1

l=04m?1?1

l=0X(k)

=x(4l)WNlk

+x(4l

+1)W

(4l+1)k

+x(4l

+

2)W

(4l+2)k

+x(4l

+3)W

(4l+3)k(0)

lk(1)

lk(2)

lk(

kN

2kN

3kN4m?14m?14m?14m?1?1

l=04m?1?1

l=04m?1?1

l=04m?1?1

l=0=k

(1)

2km?1

+WN

X

((k))

m?1

+WN

X

(2)((k))

m?1

+WNkX

(3)((k))((k))4

4

4

4m?14(0)k

(1)

2kX(k)

=

X

其中X

(0)(k)

=

DFT4m?1{x(0)(n)}X

(1)(k)

=

DFT4m?1{x(1)(n)}X

(2)(k)

=

DFT4m?1{x(2)(n)}X

(3)(k)

=

DFT4m?1{x(3)(n)}

0≤

k

≤

N

?1=

4m

?15基4DIT-FFT中的蝶形運(yùn)算：16點(diǎn)基4時(shí)域抽取的流圖6三、MATLAB程序7(1)4點(diǎn)FFT的matlab編程（蝶形運(yùn)算，將被16點(diǎn)FFT調(diào)用）：function

xk=FFT4(x);xk=zeros(1,4);

%xk初始化為0for

m=0:3

%對(duì)每一組m進(jìn)行DFT

for

n=0:3

%對(duì)每一個(gè)固定xk(m)進(jìn)行累加

xk(m+1)=xk(m+1)+x(n+1)*(-i)^(m*n);

%用累加函數(shù)表示DFT，注意matlab中矩陣從n

=1開(kāi)始

endend(2)16點(diǎn)FFT的matlab編程（調(diào)用4點(diǎn)FFT）function

[xk,k]=FFT16(x);k=0:15;N=length(x);if

N>16;

sprintf('The

length

of

x(n)

cannot

match

the

Function');

%若輸入位數(shù)大于16位，函數(shù)報(bào)錯(cuò)else

if

N<16;

for

m=N:16;

%若輸入位數(shù)小于16位，則對(duì)未定義位補(bǔ)零

x(m)=0;

end

endxk1=zeros(1,16);

%初始化xk8for

m=0:3;

%第一級(jí)x1=[x(1+m),x(5+m),x(9+m),x(13+m)];

%將所需蝶算x(n)寄存x2=FFT4(x1);

%四點(diǎn)FFT產(chǎn)生x(l)(m)for

l=0:3;xk1(1+l*4+m)=x2(l+1)*exp(-j*pi/8*m*l);

%對(duì)第一級(jí)運(yùn)算后的數(shù)乘以旋轉(zhuǎn)因子后寄存

end

endfor

m=0:3;

%對(duì)第二級(jí)基4-16位FFT外循環(huán)

x1=[xk1(1+m*4),xk1(2+m*4),xk1(3+m*4),xk1(4+m*4)];

%將所需蝶算xk1（n）寄存

x2=FFT4(x1);

%進(jìn)行蝶形運(yùn)算，相當(dāng)于FFT

for

l=0:3;xk1(1+m*4+l)=x2(l+1);

%將所的寄存解輸出到xk1中end

endfor

m=0:3;for

l=0:3xk(m*4+l+1)=xk1(l*4+m+1);

%將最終xk1(k)解進(jìn)行4進(jìn)制倒序輸出X(0),X(4),

X(8),

X(12),

X(1),

X(5),

X(9),

X(13),

X(2),

X(6),

X(10),

X(14),

X(3),

X(7),

X(11),X(15)

轉(zhuǎn)為X(0),X(1),

X(2),

X(3),

X(4),

X(5),

X(6),

X(7),

X(8),

X(9),

X(10),

X(11),

X(12),

X(13),

X(14),X(15)end

%所得xk(k)為最終DFT變換解

endend9（3）主程序代碼（調(diào)用fft16（x））x=[3

1

1

0

1

0

0

9

5

6

0

0

0

0

0

0];y=FFT16(x);subplot(2,2,1);stem(abs(y),'filled');title('X(k)mag');subplot(2,2,2);stem(angle(y),'filled');title('X(k)arg');subplot(2,2,3);stem(real(y),'filled');title('X(k)real');

10x=3110101956000000y=

Columns

1

through

527.0000-14.9343

-

3.9449i

18.3137

+

1.4142i

-7.3576

-

4.1097i

7.0000

+

2.0000iColumns

6

through

10

3.3576

-

3.2813i

-4.3137

+

1.4142i

10.9343

+

0.8835i

-5.000010.9343

-

0.8835i

Columns

11

through

15

-4.3137

-

1.4142i

3.3576

+

3.2813i