一起聊聊基于相位導數校正的高分辨率時頻分析算法 - 二階同步壓縮變換

步驚云_32

發布于 2025-8-20 07:28

瀏覽

0收藏

實現了二階同步壓縮變換算法，主要用于高分辨率時頻分析。核心功能包括：

1. 信號生成：創建具有時變特性的多分量測試信號（含高斯幅度調制和復雜頻率調制）

2. 噪聲添加：按指定信噪比添加高斯白噪聲

3. 時頻分析：

短時傅里葉變換（STFT）

標準同步壓縮變換（FSST）

二階同步壓縮變換（FSST2）

連續小波變換（CWT）

小波同步壓縮變換（WSST）

小波二階同步壓縮變換（WSST2）

4. 可視化：比較不同時頻分析方法的能量聚集性和分辨率

算法創新點在于通過引入二階相位導數校正，顯著提高了非平穩信號瞬時頻率估計的精度，特別適合分析頻率快速變化的信號。

一起聊聊基于相位導數校正的高分辨率時頻分析算法 - 二階同步壓縮變換-AI.x社區

算法步驟詳解

信號預處理：

生成測試信號（多分量非平穩信號）

添加高斯白噪聲至指定SNR

對信號進行對稱填充

STFT路徑：
a. 計算短時傅里葉變換
b. 估計瞬時頻率（一階導數）
c. 估計群延遲（時間偏移）
d. 計算相位二階導數
e. 執行二階同步壓縮：

通過相位導數校正瞬時頻率估計

將能量重分配到校正后的頻率位置

小波路徑：
a. 計算連續小波變換
b. 尺度轉頻率
c. 估計小波域瞬時頻率
d. 計算小波域相位導數
e. 執行二階同步壓縮：

校正小波域頻率估計

跨尺度能量重分配

結果可視化：

繪制6種時頻分析方法結果

比較能量聚集性和時頻分辨率

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft

# =========================================================
# 輔助函數定義
# =========================================================

def amgauss(N, t0=None, T=None):
    """
    生成高斯幅度調制信號
    Args:
        N: 點數
        t0: 時間中心 (默認N/2)
        T: 時間展寬 (默認2*sqrt(N))
    Returns:
        y: 高斯幅度調制信號
    """
    if t0 is None:
        t0 = N / 2
    if T is None:
        T = 2 * np.sqrt(N)

    tmt0 = np.arange(1, N+1) - t0
    y = np.exp(-(tmt0 / T)**2 * np.pi)
    return y

def cmor(Fb, Fc, xi):
    """
    復Morlet小波函數及其導數
    Args:
        Fb: 帶寬參數
        Fc: 中心頻率
        xi: 頻率軸向量
    Returns:
        psih: 小波函數
        dpsih: 一階導數
        ddpsih: 二階導數
    """
    psih = np.sqrt(Fb) * np.exp(-Fb**2 * np.pi * (xi - Fc)**2)
    dpsih = -2 * Fb**2 * np.pi * (xi - Fc) * psih
    ddpsih = 2 * np.pi * Fb**2 * (2 * np.pi * Fb**2 * Fc**2 - 4 * np.pi * Fb**2 * Fc * xi + 
                                  2 * np.pi * Fb**2 * xi**2 - 1) * psih
    return psih, dpsih, ddpsih

def mypad(s):
    """
    對信號進行對稱填充
    Args:
        s: 輸入信號
    Returns:
        N: 填充后長度
        x: 填充后信號
        n1: 左側填充點數
    """
    n = len(s)
    N = 2**(1 + int(np.round(np.log2(n + np.finfo(float).eps))))
    n1 = int(np.floor((N - n) / 2))

    pad_width = (n1, N - n - n1)
    x = np.pad(s, pad_width, mode='symmetric')

    return N, x, n1

def sigmerge(x1, x2, ratio):
    """
    按指定信噪比合并兩個信號
    Args:
        x1: 主信號
        x2: 噪聲信號
        ratio: 信噪比 (dB)
    Returns:
        sig: 合并后信號
    """
    x1 = np.array(x1).flatten()
    x2 = np.array(x2).flatten()

    if ratio == np.inf:
        return x1

    Ex1 = np.mean(np.abs(x1)**2)
    Ex2 = np.mean(np.abs(x2)**2)

    h = np.sqrt(Ex1 / (Ex2 * 10**(ratio / 10)))
    sig = x1 + h * x2
    return sig

一起聊聊基于相位導數校正的高分辨率時頻分析算法 - 二階同步壓縮變換-AI.x社區