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.

Delete this entry (admin only).