『壹』 基於FFT的演算法優化 要c語言完整程序（利用旋轉因子的性質），有的請留言，答謝！！！(有核心代碼，望指教
實現(C描述)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include "complex.h"
// --------------------------------------------------------------------------
#define N 8 //64
#define M 3 //6 //2^m=N
#define PI 3.1415926
// --------------------------------------------------------------------------
float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};
float x_i[N]; //N=8
/*
float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,
0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,
0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64
float x_r[N]={1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,};
float x_i[N];
*/
FILE *fp;
// ----------------------------------- func -----------------------------------
/**
* 初始化輸出虛部
*/
static void fft_init( void )
{
int i;
for(i=0; i<N; i++) x_i[i] = 0.0;
}
/**
* 反轉演算法.將時域信號重新排序.
* 這個演算法有改進的空間
*/
static void bitrev( void )
{
int p=1, q, i;
int bit_rev[ N ]; //
float xx_r[ N ]; //
bit_rev[ 0 ] = 0;
while( p < N )
{
for(q=0; q<p; q++)
{
bit_rev[ q ] = bit_rev[ q ] * 2;
bit_rev[ q + p ] = bit_rev[ q ] + 1;
}
p *= 2;
}
for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];
for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];
}
/* ------------ add by sshc625 ------------ */
static void bitrev2( void )
{
return ;
}
/* */
void display( void )
{
printf("\n\n");
int i;
for(i=0; i<N; i++)
printf("%f\t%f\n", x_r[i], x_i[i]);
}
/**
*
*/
void fft1( void )
{ fp = fopen("log1.txt", "a+");
int L, i, b, j, p, k, tx1, tx2;
float TR, TI, temp; // 臨時變數
float tw1, tw2;
/* 深M. 對層進行循環. L為當前層, 總層數為M. */
for(L=1; L<=M; L++)
{
fprintf(fp,"----------Layer=%d----------\n", L);
/* b的意義非常重大,b表示當前層的顆粒具有的輸入樣本點數 */
b = 1;
i = L - 1;
while(i > 0)
{
b *= 2;
i--;
}
// -------------- 是否外層對顆粒循環, 內層對樣本點循環邏輯性更強一些呢! --------------
/*
* outter對參與DFT的樣本點進行循環
* L=1, 循環了1次(4個顆粒, 每個顆粒2個樣本點)
* L=2, 循環了2次(2個顆粒, 每個顆粒4個樣本點)
* L=3, 循環了4次(1個顆粒, 每個顆粒8個樣本點)
*/
for(j=0; j<b; j++)
{
/* 求旋轉因子tw1 */
p = 1;
i = M - L; // M是為總層數, L為當前層.
while(i > 0)
{
p = p*2;
i--;
}
p = p * j;
tx1 = p % N;
tx2 = tx1 + 3*N/4;
tx2 = tx2 % N;
// tw1是cos部分, 實部; tw2是sin部分, 虛數部分.
tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];
tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];
/*
* inner對顆粒進行循環
* L=1, 循環了4次(4個顆粒, 每個顆粒2個輸入)
* L=2, 循環了2次(2個顆粒, 每個顆粒4個輸入)
* L=3, 循環了1次(1個顆粒, 每個顆粒8個輸入)
*/
for(k=j; k<N; k=k+2*b)
{
TR = x_r[k]; // TR就是A, x_r[k+b]就是B.
TI = x_i[k];
temp = x_r[k+b];
/*
* 如果復習一下 (a+j*b)(c+j*d)兩個復數相乘後的實部虛部分別是什麼
* 就能理解為什麼會如下運算了, 只有在L=1時候輸入才是實數, 之後層的
* 輸入都是復數, 為了讓所有的層的輸入都是復數, 我們只好讓L=1時候的
* 輸入虛部為0
* x_i[k+b]*tw2是兩個虛數相乘
*/
fprintf(fp, "tw1=%f, tw2=%f\n", tw1, tw2);
x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;
x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;
x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;
x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k, x_r[k], x_i[k]);
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k+b, x_r[k+b], x_i[k+b]);
} //
} //
} //
}
/**
* ------------ add by sshc625 ------------
* 該實現的流程為
* for( Layer )
* for( Granule )
* for( Sample )
*
*
*
*
*/
void fft2( void )
{ fp = fopen("log2.txt", "a+");
int cur_layer, gr_num, i, k, p;
float tmp_real, tmp_imag, temp; // 臨時變數, 記錄實部
float tw1, tw2;// 旋轉因子,tw1為旋轉因子的實部cos部分, tw2為旋轉因子的虛部sin部分.
int step; // 步進
int sample_num; // 顆粒的樣本總數(各層不同, 因為各層顆粒的輸入不同)
/* 對層循環 */
for(cur_layer=1; cur_layer<=M; cur_layer++)
{
/* 求當前層擁有多少個顆粒(gr_num) */
gr_num = 1;
i = M - cur_layer;
while(i > 0)
{
i--;
gr_num *= 2;
}
/* 每個顆粒的輸入樣本數N' */
sample_num = (int)pow(2, cur_layer);
/* 步進. 步進是N'/2 */
step = sample_num/2;
/* */
k = 0;
/* 對顆粒進行循環 */
for(i=0; i<gr_num; i++)
{
/*
* 對樣本點進行循環, 注意上限和步進
*/
for(p=0; p<sample_num/2; p++)
{
// 旋轉因子, 需要優化...
tw1 = cos(2*PI*p/pow(2, cur_layer));
tw2 = -sin(2*PI*p/pow(2, cur_layer));
tmp_real = x_r[k+p];
tmp_imag = x_i[k+p];
temp = x_r[k+p+step];
/*(tw1+jtw2)(x_r[k]+jx_i[k])
*
* real : tw1*x_r[k] - tw2*x_i[k]
* imag : tw1*x_i[k] + tw2*x_r[k]
* 我想不抽象出一個
* typedef struct {
* double real; // 實部
* double imag; // 虛部
* } complex; 以及針對complex的操作
* 來簡化復數運算是否是因為效率上的考慮!
*/
/* 蝶形演算法 */
x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
/* X[k] = A(k)+WB(k)
* X[k+N/2] = A(k)-WB(k) 的性質可以優化這里*/
// 旋轉因子, 需要優化...
tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);
}
/* 開跳!:) */
k += 2*step;
}
}
}
/*
* 後記:
* 究竟是顆粒在外層循環還是樣本輸入在外層, 好象也差不多, 復雜度完全一樣.
* 但以我資質愚鈍花費了不少時間才弄明白這數十行代碼.
* 從中我發現一個於我非常有幫助的教訓, 很久以前我寫過一部分演算法, 其中絕大多數都是遞歸.
* 將數據量減少, 減少再減少, 用歸納的方式來找出數據量加大代碼的規律
* 比如FFT
* 1. 先寫死LayerI的代碼; 然後再把LayerI的輸出作為LayerII的輸入, 又寫死代碼; ......
* 大約3層就可以統計出規律來. 這和遞歸也是一樣, 先寫死一兩層, 自然就出來了!
* 2. 有的功能可以寫偽代碼, 不急於求出結果, 降低復雜性, 把邏輯結果定出來後再添加.
* 比如旋轉因子就可以寫死, 就寫1.0. 流程出來後再寫旋轉因子.
* 寥寥數語, 我可真是流了不少汗! Happy!
*/
void dft( void )
{
int i, n, k, tx1, tx2;
float tw1,tw2;
float xx_r[N],xx_i[N];
/*
* clear any data in Real and Imaginary result arrays prior to DFT
*/
for(k=0; k<=N-1; k++)
xx_r[k] = xx_i[k] = x_i[k] = 0.0;
// caculate the DFT
for(k=0; k<=(N-1); k++)
{
for(n=0; n<=(N-1); n++)
{
tx1 = (n*k);
tx2 = tx1+(3*N)/4;
tx1 = tx1%(N);
tx2 = tx2%(N);
if(tx1 >= (N/2))
tw1 = -twiddle[tx1-(N/2)];
else
tw1 = twiddle[tx1];
if(tx2 >= (N/2))
tw2 = -twiddle[tx2-(N/2)];
else
tw2 = twiddle[tx2];
xx_r[k] = xx_r[k]+x_r[n]*tw1;
xx_i[k] = xx_i[k]+x_r[n]*tw2;
}
xx_i[k] = -xx_i[k];
}
// display
for(i=0; i<N; i++)
printf("%f\t%f\n", xx_r[i], xx_i[i]);
}
// ---------------------------------------------------------------------------
int main( void )
{
fft_init( );
bitrev( );
// bitrev2( );
//fft1( );
fft2( );
display( );
system( "pause" );
// dft();
return 1;
}
本文來自CSDN博客，轉載請標明出處：http://blog.csdn.net/sshcx/archive/2007/06/14/1651616.aspx
『貳』 求FFT的c語言程序
快速傅里葉變換 要用C++ 才行吧 你可以用MATLAB來實現更方便點啊
此FFT 是用VC6.0編寫，由FFT.CPP；STDAFX.H和STDAFX.CPP三個文件組成，編譯成功。程序可以用文件輸入和輸出為文件。文件格式為TXT文件。測試結果如下：
輸入文件：8.TXT 或手動輸入
8 //N
1
2
3
4
5
6
7
8
輸出結果為：或保存為TXT文件。（8OUT.TXT）
8
(36,0)
(-4,9.65685)
(-4,4)
(-4,1.65685)
(-4,0)
(-4,-1.65685)
(-4,-4)
(-4,-9.65685)
下面為FFT.CPP文件：
// FFT.cpp : 定義控制台應用程序的入口點。
#include "stdafx.h"
#include <iostream>
#include <complex>
#include <bitset>
#include <vector>
#include <conio.h>
#include <string>
#include <fstream>
using namespace std;
bool inputData(unsigned long &, vector<complex<double> >&); //手工輸入數據
void FFT(unsigned long &, vector<complex<double> >&); //FFT變換
void display(unsigned long &, vector<complex<double> >&); //顯示結果
bool readDataFromFile(unsigned long &, vector<complex<double> >&); //從文件中讀取數據
bool saveResultToFile(unsigned long &, vector<complex<double> >&); //保存結果至文件中
const double PI = 3.1415926;
int _tmain(int argc, _TCHAR* argv[])
{
vector<complex<double> > vecList; //有限長序列
unsigned long ulN = 0; //N
char chChoose = ' '; //功能選擇
//功能循環
while(chChoose != 'Q' && chChoose != 'q')
{
//顯示選擇項
cout << "\nPlease chose a function" << endl;
cout << "\t1.Input data manually, press 'M':" << endl;
cout << "\t2.Read data from file, press 'F':" << endl;
cout << "\t3.Quit, press 'Q'" << endl;
cout << "Please chose:";
//輸入選擇
chChoose = getch();
//判斷
switch(chChoose)
{
case 'm': //手工輸入數據
case 'M':
if(inputData(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
case 'f': //從文檔讀取數據
case 'F':
if(readDataFromFile(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
}
}
return 0;
}
bool Is2Power(unsigned long ul) //判斷是否是2的整數次冪
{
if(ul < 2)
return false;
while( ul > 1 )
{
if( ul % 2 )
return false;
ul /= 2;
}
return true;
}
bool inputData(unsigned long & ulN, vector<complex<double> >& vecList)
{
//題目
cout<< "\n\n\n==============================Input Data===============================" << endl;
//輸入N
cout<< "\nInput N:";
cin>>ulN;
if(!Is2Power(ulN)) //驗證N的有效性
{
cout<< "N is invalid (N must like 2, 4, 8, .....), please retry." << endl;
return false;
}
//輸入各元素
vecList.clear(); //清空原有序列
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
cout << "Input x(" << i << "):";
cin >> c;
vecList.push_back(c);
}
return true;
}
bool readDataFromFile(unsigned long & ulN, vector<complex<double> >& vecList) //從文件中讀取數據
{
//題目
cout<< "\n\n\n===============Read Data From File==============" << endl;
//輸入文件名
string strfilename;
cout << "Input filename:" ;
cin >> strfilename;
//打開文件
cout << "open file " << strfilename << "......." <<endl;
ifstream loadfile;
loadfile.open(strfilename.c_str());
if(!loadfile)
{
cout << "\tfailed" << endl;
return false;
}
else
{
cout << "\tsucceed" << endl;
}
vecList.clear();
//讀取N
loadfile >> ulN;
if(!loadfile)
{
cout << "can't get N" << endl;
return false;
}
else
{
cout << "N = " << ulN << endl;
}
//讀取元素
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
loadfile >> c;
if(!loadfile)
{
cout << "can't get enough infomation" << endl;
return false;
}
else
cout << "x(" << i << ") = " << c << endl;
vecList.push_back(c);
}
//關閉文件
loadfile.close();
return true;
}
bool saveResultToFile(unsigned long & ulN, vector<complex<double> >& vecList) //保存結果至文件中
{
//詢問是否需要將結果保存至文件
char chChoose = ' ';
cout << "Do you want to save the result to file? (y/n):";
chChoose = _getch();
if(chChoose != 'y' && chChoose != 'Y')
{
return true;
}
//輸入文件名
string strfilename;
cout << "\nInput file name:" ;
cin >> strfilename;
cout << "Save result to file " << strfilename << "......" << endl;
//打開文件
ofstream savefile(strfilename.c_str());
if(!savefile)
{
cout << "can't open file" << endl;
return false;
}
//寫入N
savefile << ulN << endl;
//寫入元素
for(vector<complex<double> >::iterator i = vecList.begin(); i < vecList.end(); i++)
{
savefile << *i << endl;
}
//寫入完畢
cout << "save succeed." << endl;
//關閉文件
savefile.close();
return true;
}
void FFT(unsigned long & ulN, vector<complex<double> >& vecList)
{
//得到冪數
unsigned long ulPower = 0; //冪數
unsigned long ulN1 = ulN - 1;
while(ulN1 > 0)
{
ulPower++;
ulN1 /= 2;
}
//反序
bitset<sizeof(unsigned long) * 8> bsIndex; //二進制容器
unsigned long ulIndex; //反轉後的序號
unsigned long ulK;
for(unsigned long p = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset<sizeof(unsigned long) * 8>(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;
ulK *= 2;
}
if(ulIndex > p)
{
complex<double> c = vecList[p];
vecList[p] = vecList[ulIndex];
vecList[ulIndex] = c;
}
}
//計算旋轉因子
vector<complex<double> > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecW.push_back(complex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW[m];
}
//計算FFT
unsigned long ulGroupLength = 1; //段的長度
unsigned long ulHalfLength = 0; //段長度的一半
unsigned long ulGroupCount = 0; //段的數量
complex<double> cw; //WH(x)
complex<double> c1; //G(x) + WH(x)
complex<double> c2; //G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength *= 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength];
c1 = vecList[j + k] + cw;
c2 = vecList[j + k] - cw;
vecList[j + k] = c1;
vecList[j + k + ulHalfLength] = c2;
}
}
}
}
void display(unsigned long & ulN, vector<complex<double> >& vecList)
{
cout << "\n\n===========================Display The Result=========================" << endl;
for(unsigned long d = 0; d < ulN;d++)
{
cout << "X(" << d << ")\t\t\t = " << vecList[d] << endl;
}
}
下面為STDAFX.H文件：
// stdafx.h : 標准系統包含文件的包含文件，
// 或是常用但不常更改的項目特定的包含文件
#pragma once
#include <iostream>
#include <tchar.h>
// TODO: 在此處引用程序要求的附加頭文件
下面為STDAFX.CPP文件：
// stdafx.cpp : 只包括標准包含文件的源文件
// FFT.pch 將成為預編譯頭
// stdafx.obj 將包含預編譯類型信息
#include "stdafx.h"
// TODO: 在 STDAFX.H 中
//引用任何所需的附加頭文件，而不是在此文件中引用
『叄』 16點DFT的FFT演算法
FFT（快速傅里葉變換）是DFT的一種特殊情況，就是當運算點的個數是2的整數次冪的時候進行的運算（不夠用0補齊）。
FFT計算原理及流程圖：
原理：FFT的計算要求點數必須為2的整數次冪，如果點數不夠用0補齊。例如計算{2，3，5，8，4}的16點FFT，需要補11個0後進行計算。FFT計算運用蝶形運算，在蝶形運算中變化規律由W(N, p)推導，其中N為FFT計算點數，J為下角標的值。
L = 1時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;
L = 2時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;
L = 3時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;
所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1
又因為2^L = 2^M*2^(L-M) = N*2^(L-M)，這里N為2的整數次冪，即N=2^M，
W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))
所以，p = J*2^(M-L)，此處J = 0, 1, ..., 2^(L-1)-1，當J遍歷結束但計算點數不夠N時，J=J+2^L，後繼續遍歷，直到計算點數為N時不再循環。
流程圖：
/*======================================================================
*方法名：fft
*方法功能：計算數組的FFT，運用蝶形運算
*
*變數名稱：
*yVector-原始數據
*length-原始數據長度
*N-FFT計算點數
*fftYreal-FFT後的實部
*fftYImg-FFT後的虛部
*
*返回值：是否成功的標志，若成功返回true，否則返回false
*=====================================================================*/
+(BOOL)fft:(floatfloat*)yVectorandOriginalLength:(NSInteger)lengthandFFTCount:(NSInteger)NandFFTReal:(floatfloat*)fftYRealandFFTYImg:(floatfloat*)fftYImg
{
//確保計算時時2的整數冪點數計算
NSIntegerN1=[selfnextNumOfPow2:N];
//定義FFT運算是否成功的標志
BOOLisFFTOK=false;
//判斷計算點數是否為2的整數次冪
if(N!=N1)
{
//不是2的整數次冪，直接計算DFT
isFFTOK=[selfdft:yVectorandOriginalLength:lengthandFFTCount:NandFFTReal:fftYRealandFFTYImg:fftYImg];
//返回成功標志
returnisFFTOK;
}
//如果計算點數位2的整數次冪，用FFT計算，如下
//定義變數
floatyVectorN[N1];//N點運算的原始數據
NSIntegerpowOfN=log2(N1);//N=2^powOfN，用於標記最大運算級數(公式中表示為：M)
NSIntegerlevel=1;//運算級數（第幾次運算），最大為powOfN，初值為第一級運算(公式中表示為：L)
NSIntegerlineNum;//行號，倒序排列後的蝶形運算行號(公式中表示為：k)
floatinverseOrderY[N1];//yVector倒序x
NSIntegerdistanceLine=1;//行間距，第level級運算每個蝶形的兩個節點距離為distanceLine=2^(L-1)(公式中表示為：B)
NSIntegerp;//旋轉因子的階數，旋轉因子表示為W(N,p)，p=J*2^(M-L)
NSIntegerJ;//旋轉因子的階數，旋轉因子表示為W(2^L,J)，J=0,1,2,...,2^(L-1)-1=distanceLine-1
floatrealTemp,imgTemp,twiddleReal,twiddleImg,twiddleTheta,twiddleTemp=PI_x_2/N1;
NSIntegerN_4=N1/4;
//判斷點數是否夠FFT運算點數
if(length<=N1)
{
//如果N至少為length，先把yVector全部賦值
for(NSIntegeri=0;i<length;i++)
{
yVectorN[i]=yVector[i];
}
if(length<N1)
{
//如果N>length後面補零
for(NSIntegeri=length;i<N1;i++)
{
yVectorN[i]=0.0;
}
}
}
else
{
//如果N<length截取相應長度的數據進行運算
for(NSIntegeri=0;i<N1;i++)
{
yVectorN[i]=yVector[i];
}
}
//調用倒序方法
[selfinverseOrder:yVectorNandN:N1andInverseOrderVector:inverseOrderY];
//初始值
for(NSIntegeri=0;i<N1;i++)
{
fftYReal[i]=inverseOrderY[i];
fftYImg[i]=0.0;
}
//三層循環
//第三層（最里）：完成相同旋轉因子的蝶形運算
//第二層（中間）：完成旋轉因子的變化，步進為2^level
//第一層（最外）：完成M次迭代過程，即計算出x(k)=A0(k),A1(k),...,Am(k)=X(k)
//第一層循環
while(level<=powOfN)
{
distanceLine=powf(2,level-1);//初始條件distanceLine=2^(level-1)
J=0;
NSIntegerpow2_Level=distanceLine*2;//2^level
NSIntegerpow2_NSubL=N1/pow2_Level;//2^(M-L)
//第二層循環
while(J<distanceLine)
{
p=J*pow2_NSubL;
lineNum=J;
NSIntegerstepCount=0;//J運算的步進計數
//求旋轉因子
if(p==0)
{
twiddleReal=1.0;
twiddleImg=0.0;
}
elseif(p==N_4)
{
twiddleReal=0.0;
twiddleImg=-1.0;
}
else
{
//計算尤拉公式中的θ
twiddleTheta=twiddleTemp*p;
//計算復數的實部與虛部
twiddleReal=cos(twiddleTheta);
twiddleImg=-11*sin(twiddleTheta);
}
//第三層循環
while(lineNum<N1)
{
//計算下角標
NSIntegerfootNum=lineNum+distanceLine;
/*---------------------------------------
*用復數運算計算每級中各行的蝶形運算結果
*X(k)=X(k)+X(k+B)*W(N,p)
*X(k+B)=X(k)-X(k+B)*W(N,p)
*---------------------------------------*/
realTemp=fftYReal[footNum]*twiddleReal-fftYImg[footNum]*twiddleImg;
imgTemp=fftYReal[footNum]*twiddleImg+fftYImg[footNum]*twiddleReal;
//將計算後的實部和虛部分別存放在返回數組中
fftYReal[footNum]=fftYReal[lineNum]-realTemp;
fftYImg[footNum]=fftYImg[lineNum]-imgTemp;
fftYReal[lineNum]=fftYReal[lineNum]+realTemp;
fftYImg[lineNum]=fftYImg[lineNum]+imgTemp;
stepCount+=pow2_Level;
//行號改變
lineNum=J+stepCount;
}
//旋轉因子的階數變換，達到旋轉因子改變的效果
J++;
}
//運算級數加一
level++;
}
isFFTOK=true;
returnisFFTOK;
}
『肆』 求FFT的c語言程序
分類: 教育/科學 >> 學習幫助
問題描述:
追20分
解析:
快速傅里葉變換 要用C++ 才行吧 你可以用MATLAB來實現更方便點啊
此FFT 是用VC6.0編寫，由FFT.CPP；STDAFX.H和STDAFX.CPP三個文件組成，編譯成功。程序可以用文件輸入和輸出為文件。文件格式為TXT文件。測試結果如下：
輸入文件：8.TXT 或手動輸入
8 N
1
2
3
4
5
6
7
8
輸出結果為：或保存為TXT文件。（8OUT.TXT）
8
(36,0)
(-4,9.65685)
(-4,4)
(-4,1.65685)
(-4,0)
(-4,-1.65685)
(-4,-4)
(-4,-9.65685)
下面為FFT.CPP文件：
FFT.cpp : 定義控制台應用程序的入口點。
#include "stdafx.h"
#include <iostream>
#include <plex>
#include <bitset>
#include <vector>
#include <conio.h>
#include <string>
#include <fstream>
using namespace std;
bool inputData(unsigned long &, vector<plex<double> >&); 手工輸入數據
void FFT(unsigned long &, vector<plex<double> >&); FFT變換
void display(unsigned long &, vector<plex<double> >&); 顯示結果
bool readDataFromFile(unsigned long &, vector<plex<double> >&); 從文件中讀取數據
bool saveResultToFile(unsigned long &, vector<plex<double> >&); 保存結果至文件中
const double PI = 3.1415926;
int _tmain(int argc, _TCHAR* argv[])
{
vector<plex<double> > vecList; 有限長序列
unsigned long ulN = 0; N
char chChoose = ' '; 功能選擇
功能循環
while(chChoose != 'Q' && chChoose != 'q')
{
顯示選擇項
cout << "\nPlease chose a function" << endl;
cout << "\t1.Input data manually, press 'M':" << endl;
cout << "\t2.Read data from file, press 'F':" << endl;
cout << "\t3.Quit, press 'Q'" << endl;
cout << "Please chose:";
輸入選擇
chChoose = getch();
判斷
switch(chChoose)
{
case 'm': 手工輸入數據
case 'M':
if(inputData(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
case 'f': 從文檔讀取數據
case 'F':
if(readDataFromFile(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
}
}
return 0;
}
bool Is2Power(unsigned long ul) 判斷是否是2的整數次冪
{
if(ul < 2)
return false;
while( ul > 1 )
{
if( ul % 2 )
return false;
ul /= 2;
}
return true;
}
bool inputData(unsigned long & ulN, vector<plex<double> >& vecList)
{
題目
cout<< "\n\n\n==============================Input Data===============================" << endl;
輸入N
cout<< "\nInput N:";
cin>>ulN;
if(!Is2Power(ulN)) 驗證N的有效性
{
cout<< "N is invalid (N must like 2, 4, 8, .....), please retry." << endl;
return false;
}
輸入各元素
vecList.clear(); 清空原有序列
plex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
cout << "Input x(" << i << "):";
cin >> c;
vecList.push_back(c);
}
return true;
}
bool readDataFromFile(unsigned long & ulN, vector<plex<double> >& vecList) 從文件中讀取數據
{
題目
cout<< "\n\n\n===============Read Data From File==============" << endl;
輸入文件名
string strfilename;
cout << "Input filename:" ;
cin >> strfilename;
打開文件
cout << "open file " << strfilename << "......." <<endl;
ifstream loadfile;
loadfile.open(strfilename.c_str());
if(!loadfile)
{
cout << "\tfailed" << endl;
return false;
}
else
{
cout << "\tsucceed" << endl;
}
vecList.clear();
讀取N
loadfile >> ulN;
if(!loadfile)
{
cout << "can't get N" << endl;
return false;
}
else
{
cout << "N = " << ulN << endl;
}
讀取元素
plex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
loadfile >> c;
if(!loadfile)
{
cout << "can't get enough infomation" << endl;
return false;
}
else
cout << "x(" << i << ") = " << c << endl;
vecList.push_back(c);
}
關閉文件
loadfile.close();
return true;
}
bool saveResultToFile(unsigned long & ulN, vector<plex<double> >& vecList) 保存結果至文件中
{
詢問是否需要將結果保存至文件
char chChoose = ' ';
cout << "Do you want to save the result to file? (y/n):";
chChoose = _getch();
if(chChoose != 'y' && chChoose != 'Y')
{
return true;
}
輸入文件名
string strfilename;
cout << "\nInput file name:" ;
cin >> strfilename;
cout << "Save result to file " << strfilename << "......" << endl;
打開文件
ofstream savefile(strfilename.c_str());
if(!savefile)
{
cout << "can't open file" << endl;
return false;
}
寫入N
savefile << ulN << endl;
寫入元素
for(vector<plex<double> >::iterator i = vecList.begin(); i < vecList.end(); i++)
{
savefile << *i << endl;
}
寫入完畢
cout << "save succeed." << endl;
關閉文件
savefile.close();
return true;
}
void FFT(unsigned long & ulN, vector<plex<double> >& vecList)
{
得到冪數
unsigned long ulPower = 0; 冪數
unsigned long ulN1 = ulN - 1;
while(ulN1 > 0)
{
ulPower++;
ulN1 /= 2;
}
反序
bitset<sizeof(unsigned long) * 8> bsIndex; 二進制容器
unsigned long ulIndex; 反轉後的序號
unsigned long ulK;
for(unsigned long p = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset<sizeof(unsigned long) * 8>(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;
ulK *= 2;
}
if(ulIndex > p)
{
plex<double> c = vecList[p];
vecList[p] = vecList[ulIndex];
vecList[ulIndex] = c;
}
}
計算旋轉因子
vector<plex<double> > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecW.push_back(plex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW[m];
}
計算FFT
unsigned long ulGroupLength = 1; 段的長度
unsigned long ulHalfLength = 0; 段長度的一半
unsigned long ulGroupCount = 0; 段的數量
plex<double> cw; WH(x)
plex<double> c1; G(x) + WH(x)
plex<double> c2; G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength *= 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength];
c1 = vecList[j + k] + cw;
c2 = vecList[j + k] - cw;
vecList[j + k] = c1;
vecList[j + k + ulHalfLength] = c2;
}
}
}
}
void display(unsigned long & ulN, vector<plex<double> >& vecList)
{
cout << "\n\n===========================Display The Result=========================" << endl;
for(unsigned long d = 0; d < ulN;d++)
{
cout << "X(" << d << ")\t\t\t = " << vecList[d] << endl;
}
}
下面為STDAFX.H文件：
stdafx.h : 標准系統包含文件的包含文件，
或是常用但不常更改的項目特定的包含文件
#pragma once
#include <iostream>
#include <tchar.h>
TODO: 在此處引用程序要求的附加頭文件
下面為STDAFX.CPP文件：
stdafx.cpp : 只包括標准包含文件的源文件
FFT.pch 將成為預編譯頭
stdafx.obj 將包含預編譯類型信息
#include "stdafx.h"
TODO: 在 STDAFX.H 中
引用任何所需的附加頭文件，而不是在此文件中引用