‘壹’ 基于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 中
引用任何所需的附加头文件，而不是在此文件中引用