快速傅利葉變換FFT

fft的作用就不多說了，搞訊號處理的人都會用上。

fft的由來：傅利葉變換ft、離散傅利葉變換dft、快速傅利葉變換fft。

學習資料：

1、陳後金的《數字訊號處理》，裡面深入淺出，該有的公式都有，程式設計思想也有。

2、一篇系統講述傅利葉變換的帖子：

3、學生對fft的理解：

4、工程人員對fft的簡單明瞭的總結：

5、英文版的闡述：

6、維基：庫利－圖基快速傅利葉變換演算法

7、一篇碉堡了的博文：

fft就是dft的快速演算法，最經典的形式莫過於cooley-tukey演算法，除此之外有多種變體，直到今天依然有人研究。

cooley-tukey演算法就是基2-fft，其基本思想是分治，把原始資料的奇偶項分成兩部分，並一直**直到剩下乙個，通過這樣做，將dft的複雜度 o(n^2)降為o( n*log2(n) )

這樣做的原理是公式推導發現前半部分和後半部分的表示式是共軛對稱的，只需要做前半部分的運算。

以下是將原始資料進行**的圖示：

8點的原始資料按 0、2、4、6 和 1、3、5、7 的奇偶順序劃分為兩部分，前半部分作為余弦分量，後半部分作為正弦分量，對應相加，這又被稱為蝶形運算，如下圖：

上面的兩部分，各自又可以繼續分解，變成：

到這裡，8點的dft就被分解為3次的蝶形運算組合。如果按照dft的公式

計算，將要進行n^n次迴圈，k和n均為[0,,n-1]的迴圈變數，而採用fft則只需要進行n*log2(n)次迴圈，外層迴圈為log2(n)次分解，內層迴圈為n個資料點運算

fft的對稱性

實訊號的頻譜時共軛對稱的，對於實序列同樣成立。由dft定義可以驗證。

dft定義：x(k)=sigma(x(n)*exp(2*pi*k*n/n)),式中n為序列長度，求和限為n=0,1,...,n-1；結果k的取值範圍為k=0,1,...,n-1。(fft為dft的快速演算法)

這裡序列序號從0到n-1,有時是從1到n，注意區別。

由此可知，x(0)=sigma(x(n))。如果x(n)為實序列，x(0)就一定是實的，代表的是序列的直流(零頻)分量。n取1和n-1,2和n-2,...時，由exp(.)指數函式的性質可以分析得到dft序列的對稱性。共軛對稱的乙個直接結果就是零頻兩側對稱的部分幅度相等【dft直接得到的是[0,2*pi]區間，根據數字頻率的2*pi週期性，可把[pi,2*pi]等價到[-pi,0]區間】。

關於dft的對稱性，數字訊號處理的書籍均有介紹，可參考：

1.丁玉美。《數字訊號處理》，西電出版

2.奧本海姆，《離散時間訊號處理》（第二版），西安交大出版

dft根據該公式編寫迴圈：

用尤拉公式進行複數運算：

為了方便直接在matlab上進行資料的生成和顯式，vc負責dft，**如下：

inputdata.m，產生輸入資料

pi = 3.141526;  
%% 產生原始資料
n = 4096;      %取樣點數
fs = 5000;    %取樣頻率
ts = 1/fs;     %取樣週期
t = 0:ts:(n-1)*ts;  %取樣範圍
xr = sin(2*pi*100*t)+sin(2*pi*750*t);       %100hz、750hz的正弦波
xr = xr + randn(size(t)); %加雜訊
plot(t,xr);
axis([t(1),t(100),1.1*min(xr),1.1*max(xr)]);
xi = zeros(1,n);
x(1:n) = double(xr);
x(n+1:2*n) = double(xi);
% %matlab 的 fft函式
% x = fft(x,n);
% r = abs(x);
% f = fs*linspace(0,1,n); 
% plot(f,r)
% grid on
% maxx = find(r==max(r));
% text(f(maxx(1)),max(r),['max=',num2str(f(maxx(1)))]);
fid = fopen('inputdata.dat','wb');
if(fid)
fwrite(fid,x,'float64');
else
msgbox('儲存失敗');
endfclose(fid);

dft.cpp，控制台程式

// dft.cpp : 定義控制台應用程式的入口點。
//#include "stdafx.h"
#include "math.h"
#include "iostream"
#include using namespace std;
void dft(double *x, double *x, int n);
int _tmain(int argc, _tchar* argv)
void complexmul(double re0, double im0, double re1, double im1, double *pre, double *pim)
// dft
// x原始資料，x變換資料
// n為資料的取樣點數
void dft(double *x, double *x, int n)
// 將時域點寫入x1  
memcpy(x1, td, sizeof(complex) * count);  
// 採用蝶形演算法進行快速付立葉變換  
for(k = 0; k < r; k++)//k為蝶形運算的級數  
{  for(j = 0; j < 1 << k; j++)  
{  bfsize = 1 << (r-k);//做蝶形運算兩點間距離  
for(i = 0; i < bfsize / 2; i++)   
{  p = j * bfsize;  
x2[i + p] = x1[i + p] + x1[i + p + bfsize / 2];  
x2[i + p + bfsize / 2] = (x1[i + p] - x1[i + p + bfsize / 2])   
* w[i * (1二維fft
第三方庫
效率

快速傅利葉變換FFT

快速傅利葉變換 FFT

快速傅利葉變換 FFT

FFT（快速傅利葉變換）

快速傅利葉變換FFT

快速傅利葉變換 FFT

快速傅利葉變換 FFT

FFT（快速傅利葉變換）

相關推薦