main.c

快速傅立叶变化的fft算法的C语言实现。

#include "include.h"
#include <math.h>

#if FFT_SAMPLE

int hd_u[3][(FFT_N/2)-2];			/* 三相电压2-15次谐波畸变率*/
int hd_i[3][(FFT_N/2)-2];			/* 三相电流2-15次谐波畸变率*/
int thd_u[3];                /* 三相电压谐波总畸变率 */
int thd_i[3];                /* 三相电流谐波总畸变率 */
int thd_u_even[3];           /* 三相电压偶次谐波畸变率 */
int thd_u_odd[3];            /* 三相电压奇次谐波畸变率 */
int thd_i_even[3];           /* 三相电流偶次谐波畸变率 */
int thd_i_odd[3];            /* 三相电流奇次谐波畸变率 */
int k_i[3];                    	/* 电流K因子 */
int k_u[3];                   	 /* 电压K因子 */


/* 2005.9.16为了解决局部变量太多导致堆栈溢出的问题，
把fft模块定义的局部变量定义为全局变量 */
double ur[FFT_N/2], ui[FFT_N/2], ir[FFT_N/2], ii[FFT_N/2];  /*电压电流的实部和虚部*/
complex k[FFT_N];		/*fft算法模块输入输出值*/
double u_even_sum, u_odd_sum, u_total_sum, i_even_sum, i_odd_sum, i_total_sum;/*电压电流畸变率计算中间值*/
double k_sum_u, sum_u, k_sum_i, sum_i;	/*K因子计算中间值*/
#endif

const int demo_u_sample[64] = 
{
		0,527,609,823,945,1126,1249,1402,
		1511,1635,1720,1817,1878,1938,1964,1995,
	1990,1981,1944,1902,1826,1755,1650,1545,
	1413,1284,1131,980,805,647,453,301,
	0,-301,-453,-647,-805,-980,-1131,-1284,
	-1413,-1545,-1650,-1755,-1826,-1902,-1944,-1981,
	-1990,-1995,-1964,-1938,-1878,-1817,-1720,-1635,
	-1511,-1402,-1249,-1126,-945,-823,-609,-527		
};

const int demo_i_sample[64] = 
{
		0,851,821,1053,1114,1289,1362,1506,
		1572,1691,1744,1832,1862,1920,1924,1951,
		1928,1925,1874,1842,1766,1700,1600,1507,
		1380,1274,1122,1001,830,699,507,401,
		0,-401,-507,-699,-830,-1001,-1122,-1274,
		-1380,-1507,-1600,-1700,-1766,-1842,-1874,-1925,
		-1928,-1951,-1924,-1920,-1862,-1832,-1744,-1691,
		-1572,-1506,-1362,-1289,-1114,-1053,-821,-851	
};

void fft_calc();

void main()
{
	

	while(1)
	{
		fft_calc();	
	}

} 

void fft_calc()
{
	unsigned char i,j,l;

	for( i = 0; i < 3;i ++ )
	{
		l = 64 / FFT_N;
		for (j=0;j<FFT_N;j++)   /* 输入的电流电压实序列，组成一个复序列*/
  		{
  			k[j].r = demo_u_sample[j*l];
  			k[j].i = demo_i_sample[j*l];
  		}

		fft( FFT_N,&k );


		/* 根据傅立叶计算公式，计算结果除以N/2*/
		for (j=1;j<(FFT_N/2);j++)    
		{
			ur[j] = (k[j].r+k[FFT_N-j].r)/FFT_N;
			ui[j] = (k[j].i-k[FFT_N-j].i)/FFT_N;
			ir[j] = (k[j].i+k[FFT_N-j].i)/FFT_N;
			ii[j] = (k[FFT_N-j].r-k[j].r)/FFT_N;

		}
		

		/*   只需要对(N/2-1)次以下的谐波进行分析    */ 
		/*   幅值X(k)=sqrt(ak^2+bk^2)          */
		for( j = 1; j <(FFT_N/2); j ++ )
		{
			k[j].r=fSqrt(ur[j]*ur[j]+ui[j]*ui[j]);
			k[j].i=fSqrt(ir[j]*ir[j]+ii[j]*ii[j]);

			/*  计算谐波有效值时幅值除以根号2 */
			k[j].r = k[j].r * 0.7071068;
			k[j].i = k[j].i * 0.7071068;
		}


		/* ********计算谐波畸变率和K因子*************/
		u_even_sum = u_odd_sum = u_total_sum = 0;
		i_even_sum  = i_odd_sum  = i_total_sum = 0;
		for ( j = 2; j < (FFT_N/2); j ++ )
		{
			u_total_sum += k[j].r * k[j].r;
			i_total_sum += k[j].i * k[j].i ;

			if( (j-1) % 2 == 0)
			{
				u_odd_sum += k[j].r * k[j].r ;
				i_odd_sum += k[j].i * k[j].i ;
			}

			if( j % 2 == 0)
			{
				u_even_sum += k[j].r * k[j].r ;
				i_even_sum += k[j].i * k[j].i ;
			}

			hd_u[i][j-2] =  k[j].r * 10000L / k[1].r ; /*各次谐波畸变率*/
			hd_i[i][j-2] =  k[j].i * 10000L / k[1].i;
			
		}

		k_sum_u= sum_u = k_sum_i= sum_i  = 0;
		for ( j = 1; j < (FFT_N/2); j ++ )
		{
			k_sum_u += k[j].r * k[j].r * j * j;
			sum_u += k[j].r * k[j].r;
			k_sum_i += k[j].i * k[j].i * j * j;
			sum_i += k[j].i * k[j].i;
		}

		/* 谐波畸变率单位%，保留1位小数*/
		thd_u[i] = fSqrt( u_total_sum ) * 10000L / k[1].r;
		thd_u_odd[i] =  fSqrt( u_odd_sum ) * 10000L / k[1].r;
		thd_u_even[i] = fSqrt( u_even_sum ) * 10000L / k[1].r;
		thd_i[i] = fSqrt( i_total_sum ) * 10000L / k[1].i;
		thd_i_odd[i] = fSqrt( i_odd_sum ) * 10000L /k[1].i;
		thd_i_even[i] = fSqrt( i_even_sum ) * 10000L / k[1].i;

		/* 计算K因子，保留1位小数*/
		k_i[i] =  ( k_sum_i / sum_i ) * 100L;    
		k_u[i] = ( k_sum_u / sum_u )* 100L; 

		/**********谐波畸变率和K因子计算结束*********/

#if 0 	/* 不需要计算基波及谐波的幅值 */
		for ( j = 2; j < (FFT_N/2); j ++ )
		{
			/* 电压保留两位小数，电流保留3位小数*/
			/* 计算结果为一次侧的值*/

			k[j].r = (k[j].r * 19250L / 10000L) * volts_scale;
			k[j].i = (k[j].i * 10728L / 10000L) * amps_scale;
		}
#endif
		
	}


}