#include "math.h"
#include "stdio.h"
float max=0,sum=0;
struct compx
{
float real;
float imag;
} compx ;
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);
}
void FFT(struct compx *xin,int N) //*xin为需要FFT变换的数组,N为转换的点数
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
long int le,B,ip;
float pi;
struct compx w,t;
LH=N>>1;
f=N;
for(m=1;(f=f>>1)!=1;m++){;}
nm=N-2;
j=N>>1;
for(i=1;i<=nm;i++) // 变址运算
{
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}
for(L=1;L<=m;L++) //实现蝶形运算
{
le=(long int)pow(2,L);
B=le>>1;
pi=3.14159;
for(j=0;j<=B-1;j++)
{
p=pow(2,m-L)*j;
ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);
for(i=j;i<=N-1;i=i+le)
{
ip=i+B;
t=EE(xin[ip],w);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
}
}
}
/********************************************************************
*功能: 失真度的计算
*
*****************************************************************/
float THDcount(unsigned int *result)
{
unsigned int i=0,n=2,m,t=0;
float THD=0;
max=result[1];
for(i=2;i<256;i++)
{
if(result[i]>max)
{
max=result[i];
t=i;
}
}
m=t*n;
while(m<256)
{
sum=result[m]*result[m]+sum;
n++;
m=t*n;
}
sum=sum/(max*max);
THD=100*sqrt(sum);
return THD;
}