# Demo entry 6362996

FFT

Submitted by anonymous on May 11, 2017 at 13:34
Language: C. Code size: 7.3 kB.

```#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define PI 3.1415926
#define N 128
#define M 7
#define A0  255.0		//定义输入波形的幅值

/*进行FFT*/
void fft(void)
{
/*输入函数：实部-dataR ,虚部-dataI*/
/*输出函数：实部-output_real ,虚部-ourput_imag*/

int i,j,k;
int p,J,L,B;
float SinIn[N];
float output_real[N];
float output_imag[N];
float dataR[N],dataI[N];
float dataA[N];
float Tr,Ti,temp;

/*将下面即将完成的fft函数的实部数据存入fft_real.txt文件中*/
/*将下面即将完成的fft函数的虚部数据存入fft_imag.txt文件中*/
FILE *fpt;
fpt = fopen("D:\word.txt","w");
if(!fpt)
{
printf("can't open word.txt file\n");

}

FILE *fp=fopen("D:\debugu.txt","w");
if(!fp)
{
printf("can't open debug.txt file\n");

}

/*输入信号波形设置*/
for(i = 0; i < N; i++)
{
SinIn[i] = A0 * (sin(2*PI*i/25)+sin(2*PI* i * 0.4 ));
dataR[i] = SinIn[i];
dataI[i] = 0;
dataA[i] = 0;
}

j = N / 2;
/*输入序列倒序*/
/*第0个数(二进制数都为0)和最后一个第N-1个数(二进制数都为1)不需倒序*/
for(i = 1; i < N - 2; i++) {
if(i < j) {
temp = dataR[i];
dataR[i] = dataR[j];
dataR[j] = temp;

//因为波形虚数部分都为0,所以不用交换
//temp = dataI[i];
//dataI[i] = dataI[j];
//dataI[j] = temp;
}
k = N / 2;
while(1) {
if(j < k) {
j = j + k;
break;
}
else {
j = j - k;
k = k / 2;
}
}
}
//进行FFT
//Xr[J] = Xr'(J) + Tr
//Xi[J] = Xi'(J) + Ti
//Xr[J+B] = Xr'(J) - Tr
//Xi[J+B] = Xi'(J) - Ti
//(其中 Xr'为上一级的Xr, Xi'为上一级的Xi)
//其中Tr = Xr'(J+B)cos(2.0*PI*p/N) + Xi'(J+B)sin(2.0*PI*p/N)
//    Ti = Xi'(J+B)cos(2.0*PI*p/N) - Xr'(J+B)sin(2.0*PI*p/N)

for(L=1; L<=M;L++)	//FFT蝶形级数L从1--M
{
//计算每个蝶形的两个输入数据相距 B = 2^(L-1);
B = 1;
i = L-1;
while(i) {
B = B * 2;
i--;
}
//或采用运算,即B = 2^(L-1);    B = (int)pow(2,L-1);

//第L级蝶形有pow(2,L-1),即2的L-1次方个蝶形运算和pow(2,L-1)个旋转因子p
for(J=0;J<=B-1;J++)		//J = 0,1,2,...,2^(L-1) - 1
{
//计算p = J*2^(M-L)
p = 1;		i = M - L;
while(i)   {p = p * 2;        i--;  }
p = J * p;
//第L级蝶形中同一个旋转因子包含多少个蝶形运算
//每个蝶形的两个输入数据相距B = 2^(L-1)个点
//同一旋转因子对应着间隔为2^L点的2^(M-L)个蝶形  (2^L = 2*B)
for(k = J; k <= N-1; k = k+2*B) {
Tr = dataR[k];
Ti = dataI[k];
temp = dataR[k+B];
dataR[k] = dataR[k] + dataR[k+B]*cos(2.0*PI*p/N)
+ dataI[k+B]*sin(2.0*PI*p/N);
dataI[k] = dataI[k] + dataI[k+B]*cos(2.0*PI*p/N)
- dataR[k+B]*sin(2.0*PI*p/N);
dataR[k+B] = Tr - dataR[k+B]*cos(2.0*PI*p/N)
- dataI[k+B]*sin(2.0*PI*p/N);
dataI[k+B] = Ti - dataI[k+B]*cos(2.0*PI*p/N)
+ temp*sin(2.0*PI*p/N);
}
}
}

for ( i=0;i<N/2;i++ )
{
//计算出输出幅值的模，具有对称性
dataA[i]=(sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]));

}
/*将经FFT输出的数据保存在txt中*/
for(i=0;i<N;i++)
{
output_real[i]=dataR[i];
output_imag[i]=dataI[i];
}
for ( i=0;i<N;i++ )
{
fprintf(fpt,"%f\n",output_real[i]);
fprintf(fp,"%f\n",(output_imag[i]));
}
//关闭文件
fclose(fpt);
fclose(fp);

}

/*读取经matlab计算得到的un（n），并计算其功率谱，并输出进txt中*/
void power()
{

float function[128];
float P[128];

FILE *fp=fopen("D:\data.txt","r");
if(!fp)
{
printf("can't open file\n");

}
FILE *fp1 = fopen("D:\power.txt","w");
if(!fp1)
{
printf("can't open file\n");

}
for(int t=0;t<N;t++)
{
fscanf(fp,"%f\n",&function[t]);
P[t] = function[t]*function[t];
printf("%f\n",P[t]);
fprintf(fp1,"%f\n",P[t]);
}
printf("\n");
fclose(fp);
fclose(fp1);
}

/*进行IFFT*/
void ifft()
{
/*输入函数：实部-input_real ,虚部-input_imag*/
/*输出函数：实部-output_real ,虚部-ourput_imag*/
int i,j,k,t;
int p,J,L,B;
float SinIn[N];
float output_real[N];
float output_imag[N];
float input_real[N];
float input_imag[N];
float dataR[N],dataI[N];
float dataA[N];
float Tr,Ti,temp;

//读取输入实部的数据
FILE *fp=fopen("D:\word.txt","r");
if(!fp)
{
printf("can't open file\n");

}
//读取输入虚部的数据
FILE *fpde=fopen("D:\debugu.txt","r");
if(!fpde)
{
printf("can't open debugu file\n");

}

//即将输出的实部到此文件中
FILE *fpt=fopen("D:\ifft.txt","w");
if(!fpt)
{
printf("can't open file\n");

}
//即将输出的虚部到此文件中
FILE *fpwf=fopen("D:\ifftfushu.txt","w");
if(!fpwf)
{
printf("can't open ifftfushu file\n");

}

//读取实部数据到input_real[]
for(t=0;t<N;t++)
{
fscanf(fp,"%f\n",&input_real[t]);
}
fclose(fp);
//读取虚部数据到input_imag[]
for(t=0;t<N;t++)
{
fscanf(fpde,"%f\n",&input_imag[t]);
}

//输入波形
for(i = 0; i < N; i++)
{
dataR[i] = input_real[i];
dataI[i] = input_imag[i];
dataA[i] = 0;
}
j = N / 2;

//输入序列倒序
//第0个数(二进制数都为0)和最后一个第N-1个数(二进制数都为1)不需倒序
for(i = 1; i < N - 2; i++) {
if(i < j) {
temp = dataR[i];
dataR[i] = dataR[j];
dataR[j] = temp;

temp = dataI[i];
dataI[i] = dataI[j];
dataI[j] = temp;
}
k = N / 2;
while(1) {
if(j < k) {
j = j + k;
break;
}
else {
j = j - k;
k = k / 2;
}
}
}
//进行IFFT
//Xr[J] = Xr'(J) + Tr
//Xi[J] = Xi'(J) + Ti
//Xr[J+B] = Xr'(J) - Tr
//Xi[J+B] = Xi'(J) - Ti
//(其中 Xr'为上一级的Xr, Xi'为上一级的Xi)
//其中Tr = Xr'(J+B)cos(2.0*PI*p/N) + Xi'(J+B)sin(2.0*PI*p/N)
//    Ti = Xi'(J+B)cos(2.0*PI*p/N) - Xr'(J+B)sin(2.0*PI*p/N)
//FFT蝶形级数L从1--M
for(L=1; L<=M;L++)
{
//计算每个蝶形的两个输入数据相距 B = 2^(L-1);
B = 1;
i = L-1;
while(i) {
B = B * 2;
i--;
}
//或采用运算,即B = 2^(L-1);    B = (int)pow(2,L-1);
//第L级蝶形有pow(2,L-1),即2的L-1次方个蝶形运算和pow(2,L-1)个旋转因子p
for(J=0;J<=B-1;J++)		//J = 0,1,2,...,2^(L-1) - 1
{
//计算p = -J*2^(M-L)
//p旋转因子与FFT相比应取负数
p = -1;		i = M - L;
while(i)   {p = p * 2;        i--;  }
p = J * p;
//第L级蝶形中同一个旋转因子包含多少个蝶形运算
//每个蝶形的两个输入数据相距B = 2^(L-1)个点
//同一旋转因子对应着间隔为2^L点的2^(M-L)个蝶形  (2^L = 2*B)
for(k = J; k <= N-1; k = k+2*B) {
Tr = dataR[k];
Ti = dataI[k];
temp = dataR[k+B];
dataR[k] = dataR[k] + dataR[k+B]*cos(2.0*PI*p/N)
+ dataI[k+B]*sin(2.0*PI*p/N);
dataI[k] = dataI[k] + dataI[k+B]*cos(2.0*PI*p/N)
- dataR[k+B]*sin(2.0*PI*p/N);
dataR[k+B] = Tr - dataR[k+B]*cos(2.0*PI*p/N)
- dataI[k+B]*sin(2.0*PI*p/N);
dataI[k+B] = Ti - dataI[k+B]*cos(2.0*PI*p/N)
+ temp*sin(2.0*PI*p/N);
}
}
}

for ( i=0;i<N/2;i++ )
{
//幅值的模与FFT相比还要除以N
dataA[i]=(sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]))/N;
}

for(i=0;i<N;i++)
{
//此处与FFT相比还要除以N
output_real[i]=dataR[i]/N;
output_imag[i]=dataI[i]/N;
}
for ( i=0;i<N;i++ )
{
printf("%f\n",output_real[i]);
fprintf(fpt,"%f\n",output_real[i]);
fprintf(fpwf,"%f\n",output_imag[i]);

}
//关闭文件
fclose(fpt);
fclose(fpde);
fclose(fpwf);
}

void main()
{
printf("fft：\n\n\n");
fft();    //fft变换
printf("功率谱：\n\n\n");
power();  //进行功率谱计算
printf("ifft：\n\n\n");
ifft();   //Ifft变换
}
```

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.