交叉谱分析——Python

最新推荐文章于 2024-11-03 11:24:59 发布

小小蒲公英

最新推荐文章于 2024-11-03 11:24:59 发布

阅读量5.4k

点赞数

CC 4.0 BY-SA版权

分类专栏： Python 数据分析

本文链接：https://blog.youkuaiyun.com/weixin_39777626/article/details/89009030

Python 同时被 2 个专栏收录

120 篇文章

订阅专栏

数据分析

9 篇文章

订阅专栏

本文介绍了一个基于交叉谱原理的程序改编实现，通过定义`spectrum`函数来计算信号的功率谱、协谱、正交谱及相位谱等，并展示了如何通过循环调用来处理不同信号对的数据。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

程序改编至交叉谱程序
由于知识受限，不知对方所用语言，仅凭个人理解改编，如有错误之处望大家指出和谅解

函数定义

def spectrum(m,x1,x2):
    pi=3.1415926
    x1bar=x1.mean()
    x2bar=x2.mean()
    s1=math.sqrt(x1.var())
    s2=math.sqrt(x2.var())
    n=len(x1)
    result=pd.DataFrame()
    M=np.array(range(m+1))
    M.dtype = 'float32'
    B=M.copy()
    B[:]=1
    B[0],B[-1]=0.5,0.5
    R11=M.copy()
    R22=M.copy()
    R12=M.copy()
    R21=M.copy()
    P11=M.copy()
    P22=M.copy()
    P12=M.copy()
    Q12=M.copy()
    M=np.array(range(m+1))
    for later in range(m+1):
        x1_series1=x1[:n-later]
        x1_series2=x1[later:]
        x2_series1=x2[:n-later]
        x2_series2=x2[later:]
        
        R11[later]=((x1_series1-x1bar)/s1*(x1_series2-x1bar)/s1).sum()/(n-later)
        R22[later]=((x2_series1-x2bar)/s2*(x2_series2-x2bar)/s2).sum()/(n-later)
        R12[later]=((x1_series1-x1bar)/s1*(x2_series2-x2bar)/s2).sum()/(n-later)
        R21[later]=((x2_series1-x2bar)/s2*(x1_series2-x1bar)/s1).sum()/(n-later)
    for later in range(m+1):
        P11[later]=R11[0]+(R11*(1+np.cos(pi*M/m))*np.cos(later*pi*M/m)).sum()
        P22[later]=R22[0]+(R22*(1+np.cos(pi*M/m))*np.cos(later*pi*M/m)).sum()
        P12[later]=R12[0]+(0.5*(R12+R21)*(1+np.cos(pi*M/m))*np.cos(later*pi*M/m)).sum()
        Q12[later]=(0.5*(1+np.cos(pi*M/m))*np.sin(pi*later*M/m)*(R12-R21)).sum()

    result['周期']=np.around(2*m/(M+1),decimals=1)
    result['功率谱_X1']=B*P11/m
    result['功率谱_X2']=B*P22/m
    result['协谱']=B*P12/m
    result['正交谱']=B*Q12/m
    result['凝聚谱']=(P12**2+Q12**2)/(P11*P22)
    result['相位谱']=180/pi*np.arctan(Q12/P12)
    return result

调用

P12=pd.DataFrame()
Q12=pd.DataFrame()
R122=pd.DataFrame()
THITA=pd.DataFrame()
for i in range(16):
    for j in range(i+1,16):
        x=np_data1[:,i+1]
        y=np_data1[:,j+1]
        Result=spectrum(20,x,y)
        P12['{}_{}'.format(i+1,j+1)]=Result['协谱']
        Q12['{}_{}'.format(i+1,j+1)]=Result['正交谱']
        R122['{}_{}'.format(i+1,j+1)]=Result['凝聚谱']
        THITA['{}_{}'.format(i+1,j+1)]=Result['相位谱']
        Result.to_csv('F:\paper\dataset\PSD\data1\{}——{}.csv'.format(i+1,j+1),encoding='gbk')
P12.to_csv('F:\paper\dataset\PSD\data1\P12.csv',encoding='gbk')
Q12.to_csv('F:\paper\dataset\PSD\data1\Q12.csv',encoding='gbk')
R122.to_csv('F:\paper\dataset\PSD\data1\R122.csv',encoding='gbk')
THITA.to_csv('F:\paper\dataset\PSD\data1\THITA.csv',encoding='gbk')


P12=pd.DataFrame()
Q12=pd.DataFrame()
R122=pd.DataFrame()
THITA=pd.DataFrame()
for i in range(5):
    x=np_data2[:,i+1]
    y=np_data2[:,6]
    Result=spectrum(39800,x,y)
    P12['{}'.format(i+1)]=Result['协谱']
    Q12['{}'.format(i+1)]=Result['正交谱']
    R122['{}'.format(i+1)]=Result['凝聚谱']
    THITA['{}'.format(i+1)]=Result['相位谱']
    Result.to_csv('F:\paper\dataset\PSD\data2\{}.csv'.format(i+1),encoding='gbk')
P12.to_csv('F:\paper\dataset\PSD\data2\P12.csv',encoding='gbk')
Q12.to_csv('F:\paper\dataset\PSD\data2\Q12.csv',encoding='gbk')
R122.to_csv('F:\paper\dataset\PSD\data2\R122.csv',encoding='gbk')
THITA.to_csv('F:\paper\dataset\PSD\data2\THITA.csv',encoding='gbk')

结果展示
在这里插入图片描述