160_dft.cpp

本程序是使用在电力谐波分析过程中160点ＤＦＴ算法

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define SAMP 160
#define frequency 47.5             //为方便仿真频率先固定
#define MN (1.0*SAMP*50/47.5)   
float tmp[21];
float ratio[20];
struct compx 
{	
	float real;
	float imag;
};
struct compx sum[SAMP+1];
float base;                 /*-基波-*/
float component_O[20];	  	/* -各次谐波分量-*/
float component_O_rt[20];
float ftmp;
float angle[20];
int N = SAMP;
#define RN SAMP*2
float AD_RESULT[RN];
float CZ_RESULT[RN];

/*float amend_mt[20]=                   //谐波修正系数
    {
		1.0,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,
        0.9752*1.08,0.9752*1.015,0.9752*1.022,0.9752*1.025,0.9752*1.030,
        0.9752*1.032,0.9752*1.042,0.9752*1.055,0.9752*1.065,0.9752*1.083
	};
*/



struct compx EE(struct compx b1,struct compx b2)
{
	struct compx b3;	
	b3.real=(b1.real*b2.real-b1.imag*b2.imag);
	b3.imag=(b1.real*b2.imag+b1.imag*b2.real);
	return(b3);
}

float pwcz(float x)	
{
	float y;
	float e1,e2,e3;
	unsigned int k;

	if((unsigned int)(x*10)%10<=5)
	{	
		k=(unsigned int)x;
		printf("k=%d\t",k);
	}
	else
	{
		k=(unsigned int)x+1;
		printf("k1=%d\t",k);
	}
	e1=(x-k)*(x-(k+1))/2.0;
	e2=(x-(k-1))*(x-(k+1))/-1.0;
	e3=(x-(k-1))*(x-k)/2.0;
	y=e1*AD_RESULT[k-1]+e2*AD_RESULT[k]+e3*AD_RESULT[k+1];
	printf("y=%f\n", y);
	return y;
}

static void draw_data(float *buf,int len,int mode)
{
    FILE *stream;
	float u;
    float temp;
	stream = fopen( "data.txt", "w" );
	int n = 0;
	while(len> 0) 
	{
		temp = *buf;
		if( temp < 0)
		{
			temp = temp/4 + 60;
		}
		else
		{
			temp = temp/4 + 60;
		}

		fprintf(stream, "%d %f\t", n++, *buf);
		for (u = 0; u < temp-1; u++)
		{
			fprintf( stream, " ");
		}
		fprintf( stream, "*\n");
		buf += 1;
		len -= 1;
	};
	fclose( stream );
}



void main()
{
	for(int n=1;n<=2*N;n++)		// 正弦波加谐波计算测试 
	{ 
		AD_RESULT[n]=220*sin(2*3.1415926*(n-1)/MN)
						+0*sin(2*2*3.1415926*(n-1)/MN)							//2次
//							+22*sin(3*2*3.1415926*(n-1)/MN)						//3次
//								+22*sin(4*2*3.1415926*(n-1)/MN)					//4次
//									+22*sin(5*2*3.1415926*(n-1)/MN)             //5次
//										+22*sin(6*2*3.1415926*(n-1)/MN)         //6次
//										+22*sin(7*2*3.1415926*(n-1)/MN)             //5次
//										+22*sin(8*2*3.1415926*(n-1)/MN)             //5次
//										+22*sin(9*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(10*2*3.1415926*(n-1)/MN)         //10次
//										+6.6*sin(11*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(12*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(13*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(14*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(15*2*3.1415926*(n-1)/MN)             //5次
//										+6.6*sin(16*2*3.1415926*(n-1)/MN)        //16次
										+11*sin(17*2*3.1415926*(n-1)/MN)
										+11*sin(18*2*3.1415926*(n-1)/MN)
										+11*sin(21*2*3.1415926*(n-1)/MN);
//										+22*sin(19*2*3.1415926*(n-1)/MN);   //19次	
//		printf("AD_RESULT[%d]=%f\n", n, AD_RESULT[n]);


	}
	for(n=1;n<=SAMP;n++)		// 正弦波加谐波计算测试 
	{
		printf("AD_RESULT[%d]=%f\t", n, AD_RESULT[n]);
		CZ_RESULT[n]=pwcz(n*50.0/frequency);
	}

	float tmp2;
	float temp;

	for(int m=1; m<22;m++)
	{
		for(n=1; n<=SAMP; n++)
		{
			sum[m].real+=CZ_RESULT[n]*cos(2*3.1415926*m*n/SAMP);
		}
		for(n=1; n<=SAMP; n++)
		{
			sum[m].imag+=CZ_RESULT[n]*sin(2*3.1415926*m*n/SAMP);
		}

		tmp2=2*sqrt((sum[m].real*sum[m].real+sum[m].imag*sum[m].imag))/SAMP;
		temp = tmp2;
		component_O[m]=temp;
		ftmp = component_O[m] * component_O[m];
		component_O_rt[m] = sqrt(ftmp);
		ratio[m] = component_O_rt[m]*100/component_O_rt[1];
//		printf("component_O_rt[%d]=%f  ratio[%d] = %f  angle[%d]=%f\n", n, component_O_rt[n], n, ratio[n], n, angle[n]);
		printf("component_O_rt[%d]=%f  ratio[%d] = %f  \n", m, component_O_rt[m], m, ratio[m]);
	} 	
	draw_data(AD_RESULT, 2*SAMP, 0);
	draw_data(CZ_RESULT, SAMP, 0);
}