我通过NI-6221采集了一个1KHz的正弦波信号,现在要通过FFT变换后计算其THDN,但结果总不对,不知是什么问题.
请大家指点下,我现在贴出相关代码:
//采样频率
#define SCANFREQ 210000
//采样点数
#define SCANCOUNT 4096
#define PI
3.1415926535897932384626433832795028841971
//定义圆周率值
/*********************************************/
// FFT 变量声明
/*********************************************/
//系数
#define
FFTOUTCOEF
(SCANCOUNT / 2) / 1.4142
//??此处为什么要除以1.4142??
//中间数据的实部
double xr[SCANCOUNT];
//中间数据的虚部
double xi[SCANCOUNT];
//旋转因子的实部
double rCf[SCANCOUNT];
//旋转因子的虚部
double iCf[SCANCOUNT];
double v_iOutaver;
/*********************************************/
// FFT 函数声明
/*********************************************/
void Calculcf();
// 计算旋转因子
void ChangeOrder(double *,double *,int);
//数据作蝶形排列
void Fft(double *,int );
void FftConvert(double *,int );
/*********************************************/
// 计算旋转因子
/*********************************************/
void Calculcf()
{
long int itemp;
for(itemp=0;itemp<SCANCOUNT;itemp++)
{
rCf[itemp]=cos(itemp*PI*2/SCANCOUNT);
iCf[itemp]=sin(itemp*PI*2/SCANCOUNT);
}
}
/*********************************************/
// 数据作蝶形排列
/*********************************************/
void ChangeOrder(double xr[], double xi[], int N)
{
int LH,N1,I,J,K;
double T;
LH= N/2;J=LH;N1=N-2;
for(I=1;I<=N1;I++)
{
if(I<J)
{
T=xr[I]; xr[I]=xr[J];xr[J]=T;
T=xi[I]; xi[I]=xi[J];xi[J]=T;
}
K=LH;
while(J>=K)
{
J=J-K;
K=(int) (K/2+0.5);
}
J=J+K;
}
}
/*********************************************/
// FFT函数实现
/*********************************************/
void Fft(double xr[],int npoint)
{
int L,B,J,P,k,N,M;
int iloop;
double rPartKB,iPartKB;
N=1;
M=0;
for(iloop=0;iloop<npoint;iloop++) xi[iloop]=0;
do{
N=2*N;
M++;
}while(N<npoint);
ChangeOrder(xr,xi,N);
for(L=1;L<=M;L++)
{
B=(int)(pow(2,(L-1))+0.5);
for(J=0;J<=B-1;J++)
{
P=J*((int)(pow(2,(M-L))+0.5));
for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5))
{
rPartKB=xr[k+B]*rCf[P]-xi[k+B]*iCf[P];
iPartKB=xi[k+B]*rCf[P]+xr[k+B]*iCf[P];
xr[k+B]=xr[k]-rPartKB;
xi[k+B]=xi[k]-iPartKB;
xr[k]=xr[k]+rPartKB;
xi[k]=xi[k]+iPartKB;
}
}
}
}
/*********************************************/
// FFT函数调用
/*********************************************/
void FftConvert(double xrt[],int npoint )
{
int iloop;
Calculcf();
for(iloop=0;iloop<npoint;iloop++) xr[iloop]=xrt[iloop];
Fft(xr,npoint);
for(iloop=0;iloop<npoint;iloop++)
xr[iloop]=sqrt(xr[iloop]*xr[iloop]+xi[iloop]*xi[iloop]) / FFTOUTCOEF;
//此处为什么要除以FFTOUTCOEF???
}
/*********************************************/
// 计算THD+N
/*********************************************/
double WINAPI PXI_THD(double *VarrayThd) //数组VarrayThd存放的是6221采集的数据
{
static double THDN;
v_iOutaver=0;
for(ii=0;ii<SCANCOUNT;ii++)
{
v_iOutaver += VarrayThd[ii];
}
v_iOutaver/=SCANCOUNT;
for(ii=0; ii<SCANCOUNT; ii++)
{
VarrayThd[ii] -= v_iOutaver;
}
FftConvert(VarrayThd,SCANCOUNT);
/*********************************************/
// 计算失真度，10个谐波
/*********************************************/
THDN=0;
for( ii=1; ii<=HarmonicNum; ii++)
{
THDN=THDN+xr[20+20*ii]*xr[20+20*ii]; //????为什么要取步长为20???
}
THDN=sqrt(THDN)/xr[20];
//????xr[20]应该为基波,为什么基波不是xr[1]呢,xr[0]为直流分量
return THDN;
} |
楼主
标题置顶
标题高亮
