一起聊聊基于相位導(dǎo)數(shù)校正的高分辨率時頻分析算法 - 二階同步壓縮變換
實現(xiàn)了二階同步壓縮變換算法，主要用于高分辨率時頻分析。核心功能包括：
1. 信號生成：創(chuàng)建具有時變特性的多分量測試信號（含高斯幅度調(diào)制和復(fù)雜頻率調(diào)制）
2. 噪聲添加：按指定信噪比添加高斯白噪聲
3. 時頻分析：
短時傅里葉變換（STFT）
標(biāo)準(zhǔn)同步壓縮變換（FSST）
二階同步壓縮變換（FSST2）
連續(xù)小波變換（CWT）
小波同步壓縮變換（WSST）
小波二階同步壓縮變換（WSST2）
4. 可視化：比較不同時頻分析方法的能量聚集性和分辨率
算法創(chuàng)新點在于通過引入二階相位導(dǎo)數(shù)校正，顯著提高了非平穩(wěn)信號瞬時頻率估計的精度，特別適合分析頻率快速變化的信號。
算法步驟詳解
信號預(yù)處理：
生成測試信號（多分量非平穩(wěn)信號）
添加高斯白噪聲至指定SNR
對信號進(jìn)行對稱填充
STFT路徑：
a. 計算短時傅里葉變換
b. 估計瞬時頻率（一階導(dǎo)數(shù)）
c. 估計群延遲（時間偏移）
d. 計算相位二階導(dǎo)數(shù)
e. 執(zhí)行二階同步壓縮：
通過相位導(dǎo)數(shù)校正瞬時頻率估計
將能量重分配到校正后的頻率位置
小波路徑：
a. 計算連續(xù)小波變換
b. 尺度轉(zhuǎn)頻率
c. 估計小波域瞬時頻率
d. 計算小波域相位導(dǎo)數(shù)
e. 執(zhí)行二階同步壓縮：
校正小波域頻率估計
跨尺度能量重分配
結(jié)果可視化：
繪制6種時頻分析方法結(jié)果
比較能量聚集性和時頻分辨率
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft
# =========================================================
# 輔助函數(shù)定義
# =========================================================
def amgauss(N, t0=None, T=None):
"""
生成高斯幅度調(diào)制信號
Args:
N: 點數(shù)
t0: 時間中心 (默認(rèn)N/2)
T: 時間展寬 (默認(rèn)2*sqrt(N))
Returns:
y: 高斯幅度調(diào)制信號
"""
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):
"""
復(fù)Morlet小波函數(shù)及其導(dǎo)數(shù)
Args:
Fb: 帶寬參數(shù)
Fc: 中心頻率
xi: 頻率軸向量
Returns:
psih: 小波函數(shù)
dpsih: 一階導(dǎo)數(shù)
ddpsih: 二階導(dǎo)數(shù)
"""
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):
"""
對信號進(jìn)行對稱填充
Args:
s: 輸入信號
Returns:
N: 填充后長度
x: 填充后信號
n1: 左側(cè)填充點數(shù)
"""
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
本文轉(zhuǎn)載自??????高斯的手稿??
贊
收藏
回復(fù)
分享
微博
QQ
微信
舉報
回復(fù)
相關(guān)推薦
