谐波分析源程序

#includestdio.h#includemath.h#includestdlib.h#defineN256#definedt0.078125//采样间隔0.078125ms#definet_min0//一个周期采样时间下限#definet_max20//一个周期采样时间上限#definePI3.14159265358979323846264338327//圆周率#definesize_x256//输入序列X[N]的大小，在本程序中仅限2的次幂typedefstruct//定义复数类型{doublereal;doubleimg;}complex;complexx[N],*W;//x[N]:FFT输入序列,*W:FFT变换核intU_t(doublet)//返回每个tms对应的U值{t=t-(long)(t)/20*20;if(t=4&&t=6){return100;}elseif(t=14&&t=16){return-100;}else{return0;}}voidadd(complexa,complexb,complex*c)//复数加法{c-real=a.real+b.real;c-img=a.img+b.img;}voidmul(complexa,complexb,complex*c)//复数乘法{c-real=a.real*b.real-a.img*b.img;c-img=a.real*b.img+a.img*b.real;}voidsub(complexa,complexb,complex*c)//复数减法{c-real=a.real-b.real;c-img=a.img-b.img;}voidchange()//变址计算，将x(n)码位倒置{complextemp;unsignedshorti=0,j=0,k=0;doublet;for(i=0;isize_x;i++){k=i;j=0;t=(log(size_x)/log(2));while((t--)0){j=j1;j|=(k&1);k=k1;}if(ji){temp=x[i];x[i]=x[j];x[j]=temp;}}}voidinitW()//初始化变换核{inti;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;isize_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}voidfft()//快速傅里叶变换{inti=0,j=0,k=0,l=0;complexup,down,product;change();for(i=0;ilog(size_x)/log(2);i++){/*一级蝶形运算*/l=1i;for(j=0;jsize_x;j+=2*l){/*一组蝶形运算*/for(k=0;kl;k++){/*一个蝶形运算*/mul(x[j+k+l],W[size_x*k/2/l],&product);add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}intmain(){doublet,sum=0,square_sum=0,average,effective;//t:时间，average：平均值，effective：有效值longcount=0;for(t=t_min;tt_max;t=t+dt){count++;sum+=U_t(t);square_sum+=U_t(t)*U_t(t);}average=sum/count;effective=sqrt(square_sum/count);printf(平均值：%f，有效值：%f\n***************************************************\n*******************谐波有效值**********************\n***************************************************\n,average,effective);for(t=t_min,count=0;countsize_x;t=t+dt,count++){x[count].real=U_t(t);x[count].img=0;}initW();fft();for(count=1;count=15;count++){printf(第%d次谐波的有效值------%f\n,count,sqrt(x[count].real*x[count].real+x[count].img*x[count].img));}}