一个关于128点的快速傅立叶的C语言程序

一个关于128点的快速傅立叶的C语言程序

这是我写的1024点的快速傅里叶变换程序，下面有验证，你把数组
double A[2049]={0};
double B[1100]={0};
double powerA[1025]={0};

改成 A[256]={0};
B[130]={0};
power[129]={0};就行了，
void FFT(double data[], int nn, int isign)

的程序可以针对任何点数，只要是2的n次方
具体程序如下：
#include <iostream.h>
#include "math.h"
#include<stdio.h>
#include<string.h>
#include <stdlib.h>
#include <fstream.h>
#include <afx.h>

void FFT(double data[], int nn, int isign)
{
//复数的快速傅里叶变换
int n,j,i,m,mmax,istep;
double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp;
n = 2 * nn;
j = 1;
for (i = 1; i<=n ; i=i+2) //这个循环进行的是码位倒置。
{
if( j > i)
{
tempr = data[j];
tempi = data[j + 1];
data[j] = data[i];
data[j + 1] = data[i + 1];
data[i] = tempr;
data[i + 1] = tempi;
}
m = n / 2;
while (m >= 2 && j > m)
{
j = j - m;
m = m / 2;
}
j = j + m;
}
mmax = 2;
while( n > mmax )
{
istep = 2 * mmax; //这里表示一次的数字的变化。也体现了级数，若第一级时，也就是书是的第0级，其为两个虚数，所以对应数组应该增加4，这样就可以进入下一组运算
theta = -6.28318530717959 / (isign * mmax);
wpr = -2.0 * sin(0.5 * theta)*sin(0.5 * theta);
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for( m = 1; m<=mmax; m=m+2)
{
for (i = m; i<=n; i=i+istep)
{
j = i + mmax;
tempr=double(wr)*data[j]-double(wi)*data[j+1];//这两句表示蝶形因子的下一个数乘以W因子所得的实部和虚部。
tempi=double(wr)*data[j+1]+double(wi)*data[j];
data[j] = data[i] - tempr; //蝶形单元计算后下面单元的实部，下面为虚部，注意其变换之后的数组序号与书上蝶形单元是一致的
data[j + 1] = data[i + 1] - tempi;
data[i] = data[i] + tempr;
data[i + 1] = data[i + 1] + tempi;
}
wtemp = wr;
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}

}

void main()
{
//本程序已经和MATLAB运算结果对比，准确无误，需要注意的的是，计算中数组都是从1开始取得，丢弃了A[0]等数据
double A[2049]={0};
double B[1100]={0};
double powerA[1025]={0};
char line[50];
char dataA[20], dataB[20];
int ij;
char ch1[3]="\t";
char ch2[3]="\n";
int strl1,strl2;
CString str1,str2;
ij=1;
//********************************读入文件data1024.txt中的数据, 其中的数据格式见该文件
FILE *fp = fopen("data1024.txt","r");
if(!fp)
{
cout<<"Open file is failing!"<<endl;
return;
}

while(!feof(fp)) //feof(fp)有两个返回值:如果遇到文件结束，函数feof（fp）的值为1，否则为0。
{
memset(line,0,50); //清空为0
memset(dataA,0,20);
memset(dataB,0,20);

fgets(line,50,fp); //函数的功能是从fp所指文件中读入n-1个字符放入line为起始地址的空间内

sscanf(line, "%s%s", dataA, dataB); //我同时读入了两列值，但你要求1024个，那么我就只用了第一列的1024个值
//dataA读入第一列，dataB读入第二列
B[ij]=atof(dataA); //将字符型的dataA值转化为float型
ij++;
}

for (int mm=1;mm<1025;mm++)//A[2*mm-1]是实部，A[2*mm]是虚部，当只要输入实数时，那么保证虚部A[mm*2]为零即可
{
A[2*mm-1]=B[mm];
A[2*mm]=0;

}
//*******************************************正式计算FFT
FFT(A,1024,1);
//********************************************写入数据到workout.txt文件中

for (int k=1;k<2049;k=k+2)
{
powerA[(k+1)/2]=sqrt(pow(A[k],2.0)+pow(A[k+1],2.0));//求功率谱

FILE *pFile=fopen("workout.txt","a+"); //?a+只能在文件最后补充，光标在结尾。没有则创建
memset(ch1,0,15);

str1.Format("%.4f",powerA[(k+1)/2]);
if (A[k+1]>=0)
str2.Format("%d\t%6.4f%s%6.4f %s",(k+1)/2,A[k],"+",A[k+1],"i");//保存fft计算的频谱，是复数频谱
else
str2.Format("%d\t%6.4f%6.4f %s",(k+1)/2,A[k],A[k+1],"i");

strl1=strlen(str1);
strl2=strlen(str2);
// 用 法:fwrite(buffer,size,count,fp);
// buffer：是一个指针，对fwrite来说，是要输出数据的地址。
// size：要写入的字节数；
// count:要进行写入size字节的数据项的个数；
// fp:目标文件指针。
fwrite(str2,1,strl2,pFile);
fwrite(ch1,1,3,pFile);
fwrite(ch1,1,3,pFile);
fwrite(str1,1,strl1,pFile);
fwrite(ch2,1,3,pFile);
fclose(pFile);
}
cout<<"计算完毕，到fft_test\workout.txt查看结果"<<endl;
}

这是我写的1024点的快速傅里叶变换程序，下面有验证，你把数组 double A[2049]={0}; double B[1100]={0}; double powerA[1025]={0}; 改成 A[256]={0}; B[130]={0}; power[129]={0};就行了， void FFT(double data[], int nn, int isign) 的程序可以针对任何点数，只要是2的n次方 具体程序如下： #include <iostream.h> #include "math.h" #include<stdio.h> #include<string.h> #include <stdlib.h> #include <fstream.h> #include <afx.h> void FFT(double data[], int nn, int isign) { //复数的快速傅里叶变换 int n,j,i,m,mmax,istep; double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp; n = 2 * nn; j = 1; for (i = 1; i<=n ; i=i+2) //这个循环进行的是码位倒置。 { if( j > i) { tempr = data[j]; tempi = data[j + 1]; data[j] = data[i]; data[j + 1] = data[i + 1]; data[i] = tempr; data[i + 1] = tempi; } m = n / 2; while (m >= 2 && j > m) { j = j - m; m = m / 2; } j = j + m; } mmax = 2; while( n > mmax ) { istep = 2 * mmax; //这里表示一次的数字的变化。也体现了级数，若第一级时，也就是书是的第0级，其为两个虚数，所以对应数组应该增加4，这样就可以进入下一组运算 theta = -6.28318530717959 / (isign * mmax); wpr = -2.0 * sin(0.5 * theta)*sin(0.5 * theta); wpi = sin(theta); wr = 1.0; wi = 0.0; for( m = 1; m<=mmax; m=m+2) { for (i = m; i<=n; i=i+istep) { j = i + mmax; tempr=double(wr)*data[j]-double(wi)*data[j+1];//这两句表示蝶形因子的下一个数乘以W因子所得的实部和虚部。 tempi=double(wr)*data[j+1]+double(wi)*data[j]; data[j] = data[i] - tempr; //蝶形单元计算后下面单元的实部，下面为虚部，注意其变换之后的数组序号与书上蝶形单元是一致的 data[j + 1] = data[i + 1] - tempi; data[i] = data[i] + tempr; data[i + 1] = data[i + 1] + tempi; } wtemp = wr; wr = wr * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } mmax = istep; } } void main() { //本程序已经和MATLAB运算结果对比，准确无误，需要注意的的是，计算中数组都是从1开始取得，丢弃了A[0]等数据 double A[2049]={0}; double B[1100]={0}; double powerA[1025]={0}; char line[50]; char dataA[20], dataB[20]; int ij; char ch1[3]="\t"; char ch2[3]="\n"; int strl1,strl2; CString str1,str2; ij=1; //********************************读入文件data1024.txt中的数据, 其中的数据格式见该文件 FILE *fp = fopen("data1024.txt","r"); if(!fp) { cout<<"Open file is failing!"<<endl; return; } while(!feof(fp)) //feof(fp)有两个返回值:如果遇到文件结束，函数feof（fp）的值为1，否则为0。 { memset(line,0,50); //清空为0 memset(dataA,0,20); memset(dataB,0,20); fgets(line,50,fp); //函数的功能是从fp所指文件中读入n-1个字符放入line为起始地址的空间内 sscanf(line, "%s%s", dataA, dataB); //我同时读入了两列值，但你要求1024个，那么我就只用了第一列的1024个值 //dataA读入第一列，dataB读入第二列 B[ij]=atof(dataA); //将字符型的dataA值转化为float型 ij++; } for (int mm=1;mm<1025;mm++)//A[2*mm-1]是实部，A[2*mm]是虚部，当只要输入实数时，那么保证虚部A[mm*2]为零即可 { A[2*mm-1]=B[mm]; A[2*mm]=0; } //*******************************************正式计算FFT FFT(A,1024,1); //********************************************写入数据到workout.txt文件中 for (int k=1;k<2049;k=k+2) { powerA[(k+1)/2]=sqrt(pow(A[k],2.0)+pow(A[k+1],2.0));//求功率谱 FILE *pFile=fopen("workout.txt","a+"); //?a+只能在文件最后补充，光标在结尾。没有则创建 memset(ch1,0,15); str1.Format("%.4f",powerA[(k+1)/2]); if (A[k+1]>=0) str2.Format("%d\t%6.4f%s%6.4f %s",(k+1)/2,A[k],"+",A[k+1],"i");//保存fft计算的频谱，是复数频谱 else str2.Format("%d\t%6.4f%6.4f %s",(k+1)/2,A[k],A[k+1],"i"); strl1=strlen(str1); strl2=strlen(str2); // 用 法:fwrite(buffer,size,count,fp); // buffer：是一个指针，对fwrite来说，是要输出数据的地址。 // size：要写入的字节数； // count:要进行写入size字节的数据项的个数； // fp:目标文件指针。 fwrite(str2,1,strl2,pFile); fwrite(ch1,1,3,pFile); fwrite(ch1,1,3,pFile); fwrite(str1,1,strl1,pFile); fwrite(ch2,1,3,pFile); fclose(pFile); } cout<<"计算完毕，到fft_test\workout.txt查看结果"<<endl; }