gst.rar_GST_S变换_信号变换_信号的s变换_地震信号

//S 变换 //参数说明: //s[][] :二维数组，用于储存时频分析的结果 //pr,pi,fr,fi:参考fft中的说明 //len : 数组长度 //sfr,sfi ：用于存储时频分析的实部和虚部 //k : len = 2**k #ifndef ST_H #define ST_H #include <cmath> #include <iostream> #include "fft.h" using std::cout; using std::cin; using std::endl; void st(float s[],float *pr,float *pi,int len ,int k,float *fr,float *fi,float sfr[],float sfi[]) { const float pai=3.1415926; const float gf = 0.05; const float gb = 4.0; const float T = 1.0/len; const float L2 = pow((1/16.0/T),2.0); float wpr[len],wpi[len],wfr[len],wfi[len]; float hfr[2*len],hfi[2*len]; float xfr[len],xfi[len]; //用于复数相乘 float tt[len],g[len]; float sum; for(int i=0;i<len/2;++i) tt[i] = -1.0*i - (gb-gf)/2.0*sqrt(L2/gf/gb); for(int j=0;j<len/2;++j) tt[j+len/2] = len/2-j - (gb-gf)/2.0*sqrt(L2/gf/gb); float X[len]; for(int i=0;i<len;++i) X[i] = ((gb+gf)/(2*gb*gf)*tt[i] + (gb-gf)/(2*gb*gf)*sqrt(tt[i]*tt[i]+L2)); kkfft(pr,pi,len,k,fr,fi,0,0); for(int i=0;i<len/2;++i) { hfr[i] = fr[i]; hfi[i] = fi[i]; hfr[i+len] = hfr[i]; hfi[i+len] = hfi[i]; } for(int m=0;m<len;++m) { if(m==0) { wfr[m] = 1.0; wfi[m] = 1.0; } else { wfr[m] = 0.0; wfi[m] = 0.0; } } //S变换开始********************************* for(int n=0;n<len;++n){ //频率 //=====将FFT后的频谱 H1 进行平移，再与 w 相乘===== for(int m=0;m<len;++m) { fr[m] = hfr[m+n]; fi[m] = hfi[m+n]; } for(int i=0;i<len;++i) { xfr[i] = fr[i]*wfr[i]-fi[i]*wfi[i]; xfi[i] = fr[i]*wfi[i]+fi[i]*wfr[i]; } //对xfr[],xfi[]进行IFFT kkfft(xfr,xfi,len,k,pr,pi,1,1); for(int i=0;i<len;++i){ if(n>len/2-1) { s[n*len+i] = 0.0; sfr[n*len+i]= 0.0; sfi[n*len+i]= 0.0; } else { s[n*len+i] = xfr[i]; sfr[n*len+i]= pr[i]; sfi[n*len+i]= pi[i]; } } //===== 计算G(m,n)=exp(-2*pai**2*m**2/n**2) ===== for(int m=0;m<len;++m) { wpr[m] = abs(n+1)/(sqrt(2.0*pai))*2.0/(gb+gf)*exp(-pow((n+1),2.0)*pow(X[m],2.0)/(2.0*pow(float(len),2.0))); } sum=0.0; for(int i=0;i<len;++i) sum = sum +wpr[i]; for(int i=0;i<len;++i){ wpr[i] = wpr[i]/(sum); wpi[i] = 0.0; } kkfft(wpr,wpi,len,k,wfr,wfi,0,0); } //S变换结束****************************************** } #endif