基于廣義高斯循環(huán)平穩(wěn)性最大化的自適應(yīng)盲反卷積脈沖提取算法（Python）
實(shí)現(xiàn)了一種基于最大廣義高斯循環(huán)平穩(wěn)性的盲反卷積算法，主要用于從噪聲背景中提取周期性脈沖信號。核心功能包括：
1.盲反卷積處理：
通過估計(jì)的逆濾波器恢復(fù)原始脈沖信號
無需先驗(yàn)知識即可分離混合信號中的周期性脈沖成分
2.關(guān)鍵算法特性：
自適應(yīng)形狀參數(shù)估計(jì)：使用廣義高斯分布模型（β參數(shù)）描述信號統(tǒng)計(jì)特性
循環(huán)平穩(wěn)性優(yōu)化：最大化輸出信號的循環(huán)平穩(wěn)性指標(biāo)（IGGCSGGS）
循環(huán)頻率自適應(yīng)：在指定帶寬內(nèi)搜索最優(yōu)循環(huán)頻率
3.主要處理流程：
信號預(yù)處理（去均值、白化）
迭代優(yōu)化FIR濾波器：
計(jì)算當(dāng)前輸出信號
估計(jì)廣義高斯形狀參數(shù)β
構(gòu)建加權(quán)協(xié)方差矩陣
求解廣義特征值問題更新濾波器
收斂判斷（基于循環(huán)平穩(wěn)性指標(biāo)變化）
4.應(yīng)用場景：
旋轉(zhuǎn)機(jī)械故障診斷（軸承/齒輪損傷檢測）
振動信號分析
周期性脈沖信號提取
強(qiáng)噪聲背景下的特征提取
詳細(xì)流程圖
流程圖說明：
1.初始化階段：
輸入觀測信號和算法參數(shù)
信號去均值、白化預(yù)處理
初始化單位向量作為FIR濾波器起點(diǎn)
2.迭代優(yōu)化核心：
信號處理：計(jì)算當(dāng)前濾波器輸出
參數(shù)估計(jì)：自適應(yīng)估計(jì)廣義高斯形狀參數(shù)β
循環(huán)分量提?。涸谥付◣拑?nèi)搜索最優(yōu)循環(huán)頻率
矩陣構(gòu)建：形成加權(quán)協(xié)方差矩陣
濾波器更新：通過廣義特征值問題求解新濾波器
3.收斂判斷：
基于循環(huán)平穩(wěn)性指標(biāo)(IGGCSGGS)的相對變化
滿足誤差閾值或達(dá)到最大迭代次數(shù)時(shí)停止
4.結(jié)果輸出：
最優(yōu)逆濾波器
提取的脈沖信號
自適應(yīng)確定的循環(huán)頻率
性能指標(biāo)值
估計(jì)的形狀參數(shù)β
def Periodic(x, alpha, fs, bd):
"""
提取信號中的循環(huán)分量
參數(shù):
x: 輸入信號
alpha: 目標(biāo)循環(huán)頻率
fs: 采樣率
bd: 帶寬
返回:
p: 循環(huán)分量
alpha1: 調(diào)整后的循環(huán)頻率
"""
L = len(x)
f = np.fft.fftfreq(L, 1/fs)
f = f[:L//2] # 只取正頻率部分
X = np.fft.fft(x) / L
AX = np.abs(X[:L//2])
t = np.arange(L) / fs
p = np.zeros(L)
alpha1 = []
# 對每個(gè)目標(biāo)頻率進(jìn)行搜索
for k in range(len(alpha)):
# 確定要搜索的目標(biāo)頻率
target_freq = alpha[k] if k == 0 else alpha1[0] * (k + 1)
# 在帶寬范圍內(nèi)搜索
idx_start = np.argmin(np.abs(f - (target_freq - bd/2)))
idx_end = np.argmin(np.abs(f - (target_freq + bd/2)))
if idx_end < idx_start:
idx_start, idx_end = idx_end, idx_start
search_range = np.arange(idx_start, min(idx_end + 1, len(f)))
if len(search_range) == 0:
# 如果沒有找到合適的頻率，使用原始alpha
alpha1.append(alpha[k])
continue
# 尋找最大幅值對應(yīng)的頻率
max_idx = search_range[np.argmax(AX[search_range])]
alpha1.append(f[max_idx])
# 重建循環(huán)分量
comp = 2 * np.real(X[max_idx] * np.exp(1j * 2 * np.pi * alpha1[k] * t))
p += comp
# 閾值處理
th = np.percentile(p, 75) # 上四分位數(shù)作為閾值
p[p < th] = 0
# 確保沒有NaN或inf
p = np.nan_to_num(p, nan=0.0, posinf=0.0, neginf=0.0)
return p, np.array(alpha1)
完整數(shù)據(jù)和詳細(xì)注釋代碼通過知乎學(xué)術(shù)咨詢獲得：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
