python 两信号互相关 两个信号的互相关函数

目录

• 互相关函数的定义
• 互相关函数的计算
• 存在具体项目参数时

互相关函数的定义

描述两个不同的信号在不同时期上的相关性的函数，主要应用：

• 混有周期成分数据（信号）的频率（周期）提取，例如两列数据在其中一列数据滞后三期时相关性最高，则该类数据的周期为3。

互相关函数的计算

公式：

$C_{1,2}(t)\approx\int_{0}^{t_{c}}v_{1}(\tau)v_{2}(t+\tau)d\tau\tag{1}$

实际计算举例，本实验中主要针对采用傅里叶变换的循环互相关，因此后续以循环互相关为例进行讲解！

循环互相关的思想如下：

注意到在计算时要用到超出原始数据索引范围的数据，其数据补充方式并不是“补零”而是“周期延拓”。

循环互相关实际则通过傅里叶变换进行计算，将上述两列数据分别变换到频率域，从而得到频率域同等长度的两列数据$[10. +0.j -2.5+3.4409548j -2.5+0.81229924j]，[10. +0.j -2.5+3.4409548j -2.5+0.81229924j]$，则输出互相关矩阵为第一个元素为：

$(10*10+0j*0j)+(10*0j-0j*10)=(100+0j)$，同理第二个元素为$（-2.5*-2.5+3.4409548j *3.4409548j ）+（-2.5*3.4409548j+2.5*3.4409548j)$

存在具体项目参数时

计算方式一：

for c in range(len(irover)):
                # ================= Non center Station Index =================
                Atr=st[irover[c]]
                # ================= Location =================
                # ================= Do FFT at non center =================
                freq_rover, time_rover, Sxx_rover_ = sg.spectrogram(Atr.data  ,
                                                                    nperseg=wd, # 窗口大小，对应fftlen
                                                                    noverlap=0, # noverlap对应overlaprate，但是noverlap是具体值，overlaprate是重叠率
                                                                    nfft=nft, # 实际就是窗口大小，表示是否对数据进行填充。
                                                                    fs=Fs, # 频率，暂时不清楚频率在计算中的作用
                                                                    scaling='spectrum', # 输出设置为功率谱
                                                                    mode='complex')
                
                # ======== Find Data in range output frequency ===========
                Sxx_rover_=Sxx_rover_[idx]
                
                # ================= Measure AutoCorrelation Ratio ================
                rat_autocorr=np.zeros((len(Sxx_pusat_[:,0]),len(Sxx_pusat_[0,:])), dtype=np.complex)
                for nwd in range(len(Sxx_pusat_[0,:])):
                    Abs_Sxx_pusat=(abs(Sxx_pusat_[:,nwd]))
                    Smooth_Abs_Sxx_pusat=sg.convolve(Abs_Sxx_pusat,smooth,mode='same')
                    
                    Abs_Sxx_rover=(abs(Sxx_rover_[:,nwd]))
                    Smooth_Abs_Sxx_rover=sg.convolve(Abs_Sxx_rover,smooth,mode='same')
                    
                    rat_autocorr[:,nwd]=((Sxx_pusat_[:,nwd]*np.conj(Sxx_rover_[:,nwd]))/
                                (Smooth_Abs_Sxx_pusat*Smooth_Abs_Sxx_rover))
                    
                time_autocorr=np.mean(rat_autocorr,axis=1)
                autocorr[:,c]=time_autocorr

计算方式二：

for(long i=0;i<d.nsta;i++){
            if(((k*d.steplen)>=d.startend[i*2])&&((k*d.steplen+d.fftlen)<d.startend[i*2+1]))
            {
                ifCC[i] = 1;
            }
            else
            {
                ifCC[i] = 0;   
            }
            
            if (ifCC[i]){
                for(long j=0;j<d.fftlen;j++){
                    in[i][j] = d.data[j+k*d.steplen+d.npts*i];
                }
                fftwf_execute(p[i]);
                if(d.ifspecwhitenning){
                    for (long j=0;j<d.nf;j++){
                        absv = sqrtf(out[i][j*d.fstride][0]*out[i][j*d.fstride][0]+out[i][j*d.fstride][1]*out[i][j*d.fstride][1]);
                        absv = max(absv,1e-8f);
                        if(absv!=0){
                            datar[j+d.nf*i] = out[i][j*d.fstride][0]/absv;
                            datai[j+d.nf*i] = out[i][j*d.fstride][1]/absv;
                        }
                    }
                }
                else{
                    for(long j=0;j<d.nf;j++){
                        datar[j+d.nf*i] = out[i][j*d.fstride][0];
                        datai[j+d.nf*i] = out[i][j*d.fstride][1];
                    }
                }
            }
        }
#ifdef useOMP
#pragma omp parallel
#pragma omp for
#endif
        for(long i=0;i<npairs; i++){
            if(ifCC[d.Pairs[i*2]]&ifCC[d.Pairs[i*2+1]]){
                for(long j=0;j<d.nf;j++){
                    d.ncfsr[j+i*d.nf] += datar[j+d.Pairs[i*2]*d.nf]*datar[j+d.Pairs[i*2+1]*d.nf] + datai[j+d.Pairs[i*2]*d.nf]*datai[j+d.Pairs[i*2+1]*d.nf];
                    d.ncfsi[j+i*d.nf] += datar[j+d.Pairs[i*2]*d.nf]*datai[j+d.Pairs[i*2+1]*d.nf] - datai[j+d.Pairs[i*2]*d.nf]*datar[j+d.Pairs[i*2+1]*d.nf];
                }
                CCnumbers[i] = CCnumbers[i] + 1;
            }
        }
    }

相同之处：均采用傅里叶变换求解
不同之处：方法二窗口的指定是在傅里叶变换函数包的外面，且其将数据集切分为每步再传入傅里叶变换函数。方法一在傅里叶变换中指定数据窗口大小和每次移动步长，暂时不知道这两种方式是否具有差异。