相信網上現在有很多關于FFT的教程,我曾經也參閱了很多網上的教程,感覺都不怎么通俗易懂。在基本上的研究FFT,并且通過編程的形式實現之后。我決定寫一篇通俗易懂的關于FFT的講解。因此我在接下來的敘述中盡量非常通俗細致的講解。
本人最早知道傅里葉變換的時候是沉迷于音樂的頻譜跳動無法自拔,當時就很想做一個音樂頻譜顯示器。搜閱了很多資料之后,才了解到傅里葉變換,和FFT。當然這都是以前的事情了,經過了系統的學習+2個星期的研究,自制了一個FFT的算法,不敢說速度上的優勢,但是個人認為是一種通俗易懂的實現方法。經過實際的VC++模擬實驗、和STM32跑的也很成功。
首先,要會FFT,就必須要對DFT有所了解,因為兩者之間本質上是一樣的。在此之前,先列出離散傅里葉變換對(DFT):
,k=0,1,…N-1
n=0,1…N-1
其中:
但是FFT之所以稱之為快速傅里葉變換,就利用了以下的幾個性質(重中之重!)
周期性:
對稱性:
可約性:
先把這仨公式放到這,接下來會用到。
根據這幾個特性,就可以將一個長的DFT運算分解為若干短序列的DFT運算的組合,從而減少運算量。在這里,為了方便理解,我就用了一個按時間抽取的快速傅里葉變換(DIT-FFT)的方法。
首先,將一個序列x(n)一分為二
對于,k=0,1,…N-1? ?設N是2的整次冪 也就是N=2^M
將x(n)按照奇偶分組
x(2r)=x1(r)
x(2r+1)=x2(r) r=0,1,…..(N/2)-1
那么將DFT也分為兩組來預算
(第一項是偶第二項是奇)
這個時候,我們上邊提的三個性質中的可約性就起到作用了:
回顧一下:
那么這個式子就可以化為:
(式1-1)
則X1(k)和X2(k)都是長度為N/2的序列 x1(k)和x2(k)的N/2點的離散傅里葉變換
其中:
K=0,1,2…N/2-1
至此,一個N點的DFT就被分解為2個N/2的DFT。但是X1(k),和X2(k)只有N/2個點,也就是說式子(1-1)只是DFT前半部分。要求DFT的后半部分可以利用其對稱性求出后半部分為:
(式1-2)
那么式(1-1)和(1-2)就可以專用一個蝶形信號流圖符號來表示。如圖:
以N=8為例,可以用下圖表示:
通過這樣的分解,每一個N/2點DFT只需(N^2)/4次復數相乘計算,明顯的節省了運算量。
以此類推,繼續將已經得出的X1(k)和X2(k)按照奇偶分解,還是和上邊一樣的方法。那么就得出了百度上都可以找到的一大推的這個圖片了。(笑)
對于這張圖片我想強調的一點就是第二階蝶形運算的時候,WN0之后不應該是WN1嗎,為什么是W2N了,這個問題之前困擾了我一段時間,但是不要忘了,前者的 W0N的展開是W0N/2因為此時N已經按照奇偶分開了,所以是N/2而W2N/2也就是W2N是根據可約性得出的,在這里不能混淆.。
對于運算效率就不用多提了
以上就是FFT算法的理論內容了,接下來就是用C語言對這個算法的實現了,對于FFT算法C語言的實現,網上的方法層出不窮,介于本人比較懶(懶得看別人的程序),再加上自給自足豐衣足食的原則,我自己也寫了一個個人認為比較通俗易懂的程序,并且為了幫助讀者理解,我特意盡量減少了庫函數的使用,一些基本的函數都是自己寫的(難免有很多BUG),但是作為FFT算法已經夠用了,目前這個程序只能處理2^N的數據,理論上來講如果不夠2^N的話可以在后面數列補0這種操作為了簡約我也就沒有實現,但是主要的功能是具備的,下面是代碼:
/*
時間:2018年11月24日
功能:將input里的數據進行快速傅里葉變換 并且輸出
*/
#include
#include
#define PI 3.1415926
#define FFT_LENGTH 8
double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};
struct complex1{ //定義一個復數結構體
double real; //實部
double image; //虛部
};
//將input的實數結果存放為復數
struct complex1 result_dat[8];
/* 虛數的乘法
*/
struct complex1 con_complex(struct complex1 a,struct complex1 b){
struct complex1 temp;
temp.real=(a.real*b.real)-(a.image*b.image);
temp.image=(a.image*b.real)+(a.real*b.image);
return temp;
}
/*
簡單的a的b次方
*/
int mypow(int a,int b){
int i,sum=a;
if(b==0)return 1;
for(i=1;i
sum*=a;
}
return sum;
}
/*
簡單的求以2為底的正整數
*/
int log2(int n){
unsigned i=1;
int sum=1;
for(i;;i++){
sum*=2;
if(sum>=n)break;
}
return i;
}
/*
簡單的交換數據的函數
*/
void swap(struct complex1 *a,struct complex1 *b){
struct complex1 temp;
temp=*a;
*a=*b;
*b=temp;
}
/*
dat為輸入數據的數組
N為抽樣次數 也代表周期 必須是2^N次方
*/
void fft(struct complex1 dat[],unsigned char N){
/*最終 dat_buf計算出 當前蝶形運算奇數項與W 乘積
dat_org存放上一個偶數項的值
*/
struct complex1 dat_buf,dat_org;
/* L為幾級蝶形運算 也代表了2進制的位數
n為當前級蝶形的需要次數 n最初為N/2 每級蝶形運算后都要/2
i j為倒位時要用到的自增符號 同時 i也用到了L碟級數 j是計算當前碟級的計算次數
re_i i_copy均是倒位時用到的變量
k為當前碟級 cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k
*/
unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
//經過觀察,發現每級蝶形運算需要N/2次運算,共運算N/2*log2N 次
unsigned char fft_counter=0;
//在此要進行補2 N必須是2^n 在此略
//蝶形級數 (L級)
L=log2(N);
//計算每級蝶形計算的次數(這里只是一個初始值) 之后每次要/2
//n=N/2;
//對dat的順序進行倒位
for(i=1;i
i_copy=i;
re_i=0;
for(j=L-1;j>0;j--){
//判斷i的副本最低位的數字 并且移動到最高位 次高位 ..
//re_i為交換的數 每次它的數字是不能移動的 并且循環之后要清0
re_i|=((i_copy&0x01)<
i_copy>>=1;
}
swap(&dat[i],&dat[re_i]);
}
//進行fft計算
for(i=0;i
fft_flag=1;
fft_counter=0;
for(j=0;j
if(fft_counter==mypow(2,i)){ //控制隔幾次,運算幾次,
fft_flag=0;
}else if(fft_counter==0){ //休止結束,繼續運算
fft_flag=1;
}
//當不判斷這個語句的時候 fft_flag保持 這樣就可以持續運算了
if(fft_flag){
dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
//計算 當前蝶形運算奇數項與W 乘積
dat_org.real=dat[j].real;
dat_org.image=dat[j].image; //暫存
dat[j].real=dat_org.real+dat_buf.real;
dat[j].image=dat_org.image+dat_buf.image;
//實部加實部 虛部加虛部
dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
//實部減實部 虛部減虛部
k++;
fft_counter++;
}else{
fft_counter--; //運算幾次,就休止幾次
k=0;
}
}
}
}
void main(){
int i;
//先將輸入信號轉換成復數
for(i=0;i
result_dat[i].image=0;
//輸入信號是二維的,暫時不存在復數
result_dat[i].real=input[i];
//result_dat[i].real=10;
//輸入信號都為實數
}
fft(result_dat,FFT_LENGTH);
for(i=0;i
input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
//取模
printf("%lf\n",input[i]);
}
while(1);
}
這個程序中input這個數組是輸入信號,在這里只模擬抽樣了8次,輸出的數據也是input,如果想看其它序列的話,可以改變FFT_LENGTH的值以及input里的內容,程序輸出的是實部和虛部的模,如果單純想看實部或者虛部的幅度的話,請自行修改程序~這就是本次淺談FFT以及FFT算法的基本實現的全部內容了參考書籍:數字信號處理
-
FFT
+關注
關注
15文章
434瀏覽量
59366 -
STM32
+關注
關注
2270文章
10895瀏覽量
355743 -
DFT
+關注
關注
2文章
231瀏覽量
22712
原文標題:高手分享:一篇通俗易懂的關于FFT的講解
文章出處:【微信號:mcu168,微信公眾號:硬件攻城獅】歡迎添加關注!文章轉載請注明出處。
發布評論請先 登錄
相關推薦
評論