welch功率谱的c语言程序 .txt

Wetch功率谱计算,采用c语言编制!

#include "math.h"
　　#include "stdlib.h"
　　#include "fft.c"
　　
　　void pmpse(double 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+1)*sizeof(double));
　　 ximag = (double *)malloc((nfft+1)*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;
　　 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;
　　 }//k循环结束
　　 norm = 2.0*u*NumOfSections;
　　 for(i=0;i　　 sxx[i]=sxx[i]/norm;
　　 xreal[i]=sxx[i];
　　 ximag[i]=0.0;
　　 j=nfft-i;
　　 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]=20.0*log10(sxx[i]);
　　 }
　　 }
　　 free(xa);
　　 free(xreal);
　　 free(ximag);
　　 free(window); 
　　}