交叉谱分析，凝聚谱始终为1

渣渣顼

各位前辈，我在按照书上写交叉谱的时候，做到凝聚谱，发现我按书上公式来算分子分母相等。我觉得肯定是我理解有问题，麻烦大家帮我看看哪里有问题吧。代码如下
def comput_cross_spe(u, v):
    u_fft = np.fft.fftshift(np.fft.fft(u))
    v_fft = np.fft.fftshift(np.fft.fft(v))
    a1 = u_fft.real
    a2 = v_fft.real
    b1 = u_fft.imag
    b2 = v_fft.imag
    P = (a1*a2 + b1*b2)/4
    Q = (a1*b2 - a2*b1)/4
    P11 = (a1*a1 + b1*b1)/4
    P22 = (a2*a2 + b2*b2)/4
    R = (P**2 + Q**2)/(P11*P22)
    print(R)
    T = np.arctan(Q/P)
    return [P, Q, R, T]

渣渣顼

问题解决了，需要平滑才行

寒夜未明

渣渣顼 发表于 2020-7-30 15:44
问题解决了，需要平滑才行

楼主可以分享下平滑的方法吗

渣渣顼

寒夜未明 发表于 2021-8-6 11:23
楼主可以分享下平滑的方法吗

def comput_cross_spe(u, v):
    u_fft = np.fft.fftshift(np.fft.fft(u))
    v_fft = np.fft.fftshift(np.fft.fft(v))
    a1 = u_fft.real
    a2 = v_fft.real
    b1 = u_fft.imag
    b2 = v_fft.imag
    P = (a1 * a2 + b1 * b2) / 4
    Q = (a1 * b2 - a2 * b1) / 4
    P11 = smooth_operation((a1 * a1 + b1 * b1) / 4)
    P22 = smooth_operation((a2 * a2 + b2 * b2) / 4)
    P = smooth_operation(P)
    Q = smooth_operation(Q)
    R = (P ** 2 + Q ** 2) / (P11 * P22)
    T = np.arctan(Q / P)
    return [P, Q, R, T]

def smooth_operation(P):
    P_size = len(P)
    for i in range(P_size):
        if i == 0:
            P[0] = (P[0] + P[1]) * 0.5
        elif i == P_size - 1:
            P = (P + P[i - 1]) * 0.5
        else:
            P = 0.25 * P[i - 1] + 0.5 * P + 0.25 * P[i + 1]
    return P