pmp matlab 代码 【welch功率谱估计】C语言实现的代码与matlab计算结果不同

其中FFT部分已经调试通过，只是功率谱值的计算结果与matlab出入，求高手告知。

/*

x-双精度实型一维数组，长度为n，输入信号x(i)

n-整形变量，输入信号的长度

m-整形变量，分段的长度

nfft-整形变量，估计功率谱所用FFT的长度，它必须是2的整数次幂且nfft>=m

win-整形变量，窗函数的类型。win=1，表示矩形窗；win=2，表示海明窗

fs-双精度实型变量，采样频率

r-双精度实型一维数组，长度为(nfft/2+1),用于存放相关函数r(i)的值

freq-双精度实型一维数组，长度为(nfft/2+1),用于存放与功率谱相对应的频率值

sxx-双精度实型一维数组，长度为(nfft/2+1),用于存放功率谱的值

sdb-整形变量，用于表示功率谱的类型，sdb=0，表示线形谱；sdb=1，表示以dB为单位的对数谱*/

#include "math.h"

#include "stdlib.h"

#include "fft.h"

void pmpse(int x[],int n,int m,int nfft,int win,double fs,double r[],double freq[],double sxx[],int sdb)

{

int i,j,k,s;

int m2,nrd,kmax,ns1,nsectp;

int nfft21,NumOfSections,NumUsed;

double u,fl,xsum,norm,twopi,rexmn,imxmn,xmean;

double *xa,*xreal,*ximag,*window;

xa=(double*)malloc(nfft*sizeof(double));

xreal=(double*)malloc(nfft*sizeof(double));

ximag=(double*)malloc(nfft*sizeof(double));

window=(double*)malloc(m*sizeof(double));

nfft21=nfft/2+1;

NumOfSections=(n-m/2)/(m/2);

NumUsed=NumOfSections*(m/2)+m/2;

s=0;

xsum=0.0;

ns1=NumOfSections+1;//分配的段数

m2=m/2;

for(k=0;k

{

for(i=0;i

{

xa[i]=x[s+i];

}

for(i=0;i

{

xsum=xsum+xa[i];

}

s+=m2;

}

xmean=xsum/NumUsed;

rexmn=xmean;

imxmn=xmean;

u=(double)m;//计算u，使用海明窗

if(win==2)

{

u=0.0;

twopi=8.0*atan(1.0);

fl=m-1.0;

for (i = 0;i

{

window[i]=0.54-0.46*cos(twopi*i/fl);

u+=window[i]*window[i];

}

}

s=0;

for(i=0;i

{sxx[i]=0.0;

}

m2=m/2;

for(i=0;i

{xa[i+m2]=x[s+i];}

s+=m2;

kmax=(NumOfSections+1)/2;

nsectp=(NumOfSections+1)/2;

nrd=m;

for(k=0;k

{

for(i=0;i

{

j=m2+i;

xreal[i]=xa[j];

ximag[i]=0.0;

}

if((k==(kmax-1))&&(nsectp!=NumOfSections))

{

for(i=m2;i

{xa[i]=0.0;}

nrd=m/2;

}

for(i=0;i

{xa[i]=x[s+i];}

for(i=0;i

{

j=m2+i;

xreal[j]=xa[i]-rexmn;

ximag[j]=xa[j]-imxmn;

xreal[i]=xreal[i]-rexmn;

ximag[i]=xa[i]-imxmn;

}

if((k==(kmax-1))&&(nsectp!=NumOfSections))

{

for(i=0;i

{ximag[i]=0.0;}

}

s=s+nrd;

if(win==2)//给每段加窗函数

{

for(i=0;i

{

xreal[i]*=window[i];

ximag[i]*=window[i];

}

}

if(m!=nfft)//不够的补零

{

for(i=m;i

{

xreal[i]=0.0;

ximag[i]=0.0;

}

}

fft(xreal,ximag,nfft,1);

for(i=1;i

{

j=nfft-i;

sxx[i]+=xreal[i]*xreal[i]+ximag[i]*ximag[i];

sxx[i]+=xreal[j]*xreal[j]+ximag[j]*ximag[j];

}

sxx[0]+=xreal[0]*xreal[0]*2.0;

sxx[0]+=ximag[0]*ximag[0]*2.0;

}

norm=2.0*u*NumOfSections;

for(i=0;i<=nfft21;i++)

{

sxx[i]=sxx[i]/norm;

xreal[i]=sxx[i];

ximag[i]=0.0;

j=nfft-i-1;

xreal[j]=xreal[i];

ximag[j]=ximag[i];

}

fft(xreal,ximag,nfft,-1);

for(i=0;i

{r[i]=xreal[i];}

for(i=0;i

{

freq[i]=i*fs/(double)nfft;

if(sdb==1)

{

if(sxx[i]==0.0)

sxx[i]=1.0e-15;

sxx[i]=10.0*log10(sxx[i]);

}

}

free(xa);

free(xreal);

free(ximag);

free(window);

}