welch功率谱的c语言程序

#includemath.h¡¡¡¡#includestdlib.h¡¡¡¡#includefft.c¡¡¡¡¡¡¡¡voidpmpse(doublex[],intn,intm,intnfft,intwin,doublefs,doubler[],doublefreq[],doublesxx[],intsdb)¡¡¡¡{¡¡¡¡inti,j,k,s;¡¡¡¡intm2,nrd,kmax,ns1,nsectp;¡¡¡¡intnfft21,NumOfSections,NumUsed;¡¡¡¡doubleu,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);¡¡¡¡}