a.cpp

可以实现FFT信号频谱分析 及fir 数字滤波

// a.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#ifndef M_PI
#define	M_PI		3.14159265358979323846  /* pi */
#endif

/*
* This is the function you will use the most often.
* Given an array of floats, this will compute the power
* spectrum by doing a Real FFT and then computing the
* sum of the squares of the real and imaginary parts.
* Note that the output array is half the length of the
* input array, and that NumSamples must be a power of two.
*/

void PowerSpectrum(int NumSamples, float *In, float *Out);

/*
* Computes an FFT when the input data is real but you still
* want complex data as output.  The output arrays are half
* the length of the input, and NumSamples must be a power of
* two.
*/

void RealFFT(int NumSamples,
			 float *RealIn, float *RealOut, float *ImagOut);

/*
* Computes a FFT of complex input and returns complex output.
* Currently this is the only function here that supports the
* inverse transform as well.
*/

void FFT(int NumSamples,
		 bool InverseTransform,
		 float *RealIn, float *ImagIn, float *RealOut, float *ImagOut);

/*
* Applies a windowing function to the data in place
*
* 0: Rectangular (no window)
* 1: Bartlett    (triangular)
* 2: Hamming
* 3: Hanning
*/

void WindowFunc(int whichFunction, int NumSamples, float *data);

/*
* Returns the name of the windowing function (for UI display)
*/

const char *WindowFuncName(int whichFunction);

/*
* Returns the number of windowing functions supported
*/

int NumWindowFuncs();
#if 1
int _tmain(int argc, _TCHAR* argv[])
{
	int NumSamples=25000;
	int M=1;	
	float *Real2   =new float[NumSamples];
	float *RealIn  =new float[NumSamples];
	float *RealIn1 =new float[NumSamples];
	float *RealIn2 =new float[NumSamples];
	float *RealRs  =new float[NumSamples];
	float *RealRc  =new float[NumSamples];
	float *RealOut =new float[NumSamples];
	
	float *err  =new float[NumSamples];
	float *ws=new float[M];
	float *wc=new float[M];
	float mu=0.0005f;
	float a1=6.0f;
	float a2=24.0f;
	float ms=5.0f;
	for(int i=0;i<M;i++)
	{ 
		ws[i]=(float)(rand()%M-rand()%M)/100.0f;
		wc[i]=(float)(rand()%M-rand()%M)/100.0f;
		RealRs[i]=sin(2.0f*M_PI*0.001f*(float)i*ms);
		RealRc[i]=cos(2.0f*M_PI*0.001f*(float)i*ms);
		Real2[i]=sin(2.0f*M_PI*0.0025f*(float)i*ms);
		RealIn1[i]=a1*sin(2.0f*M_PI*0.001f*(float)i*ms);
		RealIn2[i]=a2*sin(2.0f*M_PI*0.00085f*(float)i*ms);
		RealIn[i]=RealIn1[i]+RealIn2[i]+4.0f*sin(2.0f*M_PI*0.0025f*(float)i*ms)+2.0f*(float)(rand()%20-10)/10.0f;
		err[i]=0;
		RealOut[i]=0;
	}
	for(int i=M;i<NumSamples;i++)
	{
		RealRs[i]=sin(2.0f*M_PI*0.001f*(float)i*ms);
		RealRc[i]=cos(2.0f*M_PI*0.001f*(float)i*ms);
		RealIn1[i]=a1*sin(2.0f*M_PI*0.001f*(float)i*ms);
		RealIn2[i]=a2*sin(2.0f*M_PI*0.00085f*(float)i*ms);
		RealIn[i]=RealIn1[i]+RealIn2[i]+4.0f*sin(2.0f*M_PI*0.0025f*(float)i*ms)+2.0f*(float)(rand()%20-10)/10.0f;
			
		RealOut[i]=0.0f;
		for(int j=0;j<M;j++)
			RealOut[i]+=RealRs[i-M+j]*ws[j]+RealRc[i-M+j]*wc[j];
		err[i]=RealIn[i]-RealOut[i];
		for(int j=0;j<M;j++)
		{
			ws[j]=ws[j]+mu*(RealRs[i-M+j])*err[i];
			wc[j]=wc[j]+mu*(RealRc[i-M+j])*err[i];
		}
	}

	FILE* f=fopen("a.csv","w+");
	if(f)
	{
		fprintf(f,"%s\n","sin1hz,sin0.85,1+0.85,out,err,1hzsin,1hzcos");
		for(int i=0;i<NumSamples;i++)
		{
			fprintf(f,"%.05f,",RealIn1[i]);
			fprintf(f,"%.05f,",RealIn2[i]);
			fprintf(f,"%.05f,",RealIn[i]);
			fprintf(f,"%.05f,",RealOut[i]);
			fprintf(f,"%.05f,",err[i]);
			fprintf(f,"%.05f,",RealRs[i]);
			fprintf(f,"%.05f,",RealRc[i]);

			fwrite("\n",1,1,f);
		}
		fclose(f);
	}
	delete []RealIn;
	delete []RealRs;
	delete []RealRc;
	delete []RealOut;
	delete []ws;
	delete []wc;
	return 0;
}

#else
int _tmain(int argc, _TCHAR* argv[])
{
//ȡpda
	FILE *data=fopen("E:\\fft\\a\\x.txt","r+");
	if(data==0)
		return 1; 
	int NumSamples=6000;
	int M=1;	
	float *RealIn  =new float[NumSamples];
	float *RealIn0s =new float[NumSamples];
	float *RealIn0c =new float[NumSamples];
	float *RealIn1s =new float[NumSamples];
	float *RealIn1c =new float[NumSamples];
	float *RealIn2s =new float[NumSamples];
	float *RealIn2c =new float[NumSamples];
	float *RealIn3s =new float[NumSamples];
	float *RealIn3c =new float[NumSamples];
	float *RealIn4s =new float[NumSamples];
	float *RealIn4c =new float[NumSamples];
	float *RealIn5s =new float[NumSamples];
	float *RealIn5c =new float[NumSamples];

	float *RealOut  =new float[NumSamples];
	float *RealOut2  =new float[NumSamples];
	float *RealOut3  =new float[NumSamples];
	
	float *err  =new float[NumSamples];
	float *w0s=new float[M];
	float *w0c=new float[M];
	float *w1s=new float[M];
	float *w1c=new float[M];
	float *w2s=new float[M];
	float *w2c=new float[M];
	float *w3s=new float[M];
	float *w3c=new float[M];
	float *w4s=new float[M];
	float *w4c=new float[M];
	float *w5s=new float[M];
	float *w5c=new float[M];

	float mu=0.0003f;
	float theta=0.0f;
	char buffer[265];

	fscanf(data,"%s\t\t%s\n",buffer,buffer);
	for(int i=0;i<M;i++)
	{ 
		float ff,fs;
		fscanf(data,"%f\t\t%f\n",&fs,&ff);
		w0s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w0c[i]=(float)(rand()%M-rand()%M)/100.0f;
		w1s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w1c[i]=(float)(rand()%M-rand()%M)/100.0f;
		w2s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w2c[i]=(float)(rand()%M-rand()%M)/100.0f;
		w3s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w3c[i]=(float)(rand()%M-rand()%M)/100.0f;
		w4s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w4c[i]=(float)(rand()%M-rand()%M)/100.0f;
		w5s[i]=(float)(rand()%M-rand()%M)/100.0f;
		w5c[i]=(float)(rand()%M-rand()%M)/100.0f;
		theta+=fs/60.0f/1250.0f*2.0f*M_PI;
		RealIn0s[i] =1.0f*sin(theta*0.85f);
		RealIn0c[i] =1.0f*cos(theta*0.85f);

		RealIn1s[i] =1.0f*sin(theta);
		RealIn1c[i] =1.0f*cos(theta);
		RealIn2s[i]=0.5f*sin(2.0f*theta);
		RealIn2c[i]=0.5f*cos(2.0f*theta);
		RealIn3s[i]=0.5f*sin(3.0f*theta);
		RealIn3c[i]=0.5f*cos(3.0f*theta);
		RealIn4s[i]=0.5f*sin(4.0f*theta);
		RealIn4c[i]=0.5f*cos(4.0f*theta);
		RealIn5s[i]=0.5f*sin(5.0f*theta);
		RealIn5c[i]=0.5f*cos(5.0f*theta);


		RealIn[i]=ff;
		err[i]=0;
		RealOut[i]=0;
		RealOut2[i]=0;
	}

	for(int i=M;i<NumSamples;i++)
	{
		float ff,fs;
		fscanf(data,"%f\t\t%f\n",&fs,&ff);
		theta+=fs/60.0f/1250.0f*2.0f*M_PI;
		RealIn0s[i] =1.0f*sin(theta*0.85f);
		RealIn0c[i] =1.0f*cos(theta*0.85f);

		RealIn1s[i] =1.0f*sin(theta);
		RealIn1c[i] =1.0f*cos(theta);
		RealIn2s[i]=0.5f*sin(2.0f*theta);
		RealIn2c[i]=0.5f*cos(2.0f*theta);
		RealIn3s[i]=0.5f*sin(3.0f*theta);
		RealIn3c[i]=0.5f*cos(3.0f*theta);
		RealIn4s[i]=0.5f*sin(4.0f*theta);
		RealIn4c[i]=0.5f*cos(4.0f*theta);
		RealIn5s[i]=0.5f*sin(5.0f*theta);
		RealIn5c[i]=0.5f*cos(5.0f*theta);

		RealIn[i]=ff;
		RealOut[i]=0.0f;
		for(int j=0;j<M;j++)	
		{
			RealOut[i]+=RealIn1s[i-M+j]*w1s[j]+RealIn1c[i-M+j]*w1c[j];
			RealOut[i]+=RealIn2s[i-M+j]*w2s[j]+RealIn2c[i-M+j]*w2c[j];
			RealOut[i]+=RealIn3s[i-M+j]*w3s[j]+RealIn3c[i-M+j]*w3c[j];
			RealOut[i]+=RealIn4s[i-M+j]*w4s[j]+RealIn4c[i-M+j]*w4c[j];
			RealOut[i]+=RealIn5s[i-M+j]*w5s[j]+RealIn5c[i-M+j]*w5c[j];
			//RealOut[i]+=RealIn0s[i-M+j]*w0s[j]+RealIn0c[i-M+j]*w0c[j];
		 
		}
		err[i]=RealIn[i]-RealOut[i];
		for(int j=0;j<M;j++)
		{
			w1s[j] =w1s[j] +mu*(RealIn1s[i-M+j])*err[i];
			w1c[j ]=w1c[j] +mu*(RealIn1c[i-M+j])*err[i];
			w2s[j] =w2s[j] +mu*(RealIn2s[i-M+j])*err[i];
			w2c[j] =w2c[j] +mu*(RealIn2c[i-M+j])*err[i];
			w3s[j] =w3s[j] +mu*(RealIn3s[i-M+j])*err[i];
			w3c[j] =w3c[j] +mu*(RealIn3c[i-M+j])*err[i];
			w4s[j] =w4s[j] +mu*(RealIn4s[i-M+j])*err[i];
			w4c[j] =w4c[j] +mu*(RealIn4c[i-M+j])*err[i];
			w5s[j] =w5s[j] +mu*(RealIn5s[i-M+j])*err[i];
			w5c[j] =w5c[j] +mu*(RealIn5c[i-M+j])*err[i];
			w0s[j] =w0s[j] +mu*(RealIn0s[i-M+j])*err[i];
			w0c[j] =w0c[j] +mu*(RealIn0c[i-M+j])*err[i];
		}
	}
	fclose(data);
	FILE* f=fopen("a.csv","w+");
	if(f)
	{
		for(int i=0;i<NumSamples;i++)
		{
			fprintf(f,"%.05f,",RealIn[i]);
			fprintf(f,"%.05f,",RealOut[i]);
			fprintf(f,"%.05f,",err[i]);
			fprintf(f,"%.05f,",RealIn1s[i]);
			fprintf(f,"%.05f,",RealIn1c[i]);

			fwrite("\n",1,1,f);
		}
		fclose(f);
	}
	delete []RealIn;
	delete []RealIn0s;
	delete []RealIn0c;
	delete []RealIn1s;
	delete []RealIn1c;
	delete []RealIn2s;
	delete []RealIn2c;
	delete []RealIn3s;
	delete []RealIn3c;
	delete []RealIn4s;
	delete []RealIn4c;
	delete []RealIn5s;
	delete []RealIn5c;

	delete []RealOut;
	delete []RealOut2;
	delete []RealOut3;
	delete []w0s;
	delete []w0c;
	delete []w1s;
	delete []w1c;
	delete []w2s;
	delete []w2c;
	delete []w3s;
	delete []w3c;
	delete []w4s;
	delete []w4c;
	delete []w5s;
	delete []w5c;
	return 0;
}

#endif
/**********************************************************************

FFT.cpp

Dominic Mazzoni

September 2000

This file contains a few FFT routines, including a real-FFT
routine that is almost twice as fast as a normal complex FFT,
and a power spectrum routine when you know you don't care
about phase information.

Some of this code was based on a free implementation of an FFT
by Don Cross, available on the web at:

http://www.intersrv.com/~dcross/fft.html

The basic algorithm for his code was based on Numerican Recipes
in Fortran.  I optimized his code further by reducing array
accesses, caching the bit reversal table, and eliminating
float-to-double conversions, and I added the routines to
calculate a real FFT and a real power spectrum.

**********************************************************************/




int **gFFTBitTable = NULL;
const int MaxFastBits = 16;

/* Declare Static functions */
static int IsPowerOfTwo(int x);
static int NumberOfBitsNeeded(int PowerOfTwo);