頻超分辨率的信號(hào)指紋高清修復(fù))
1. 項(xiàng)目概述從“模糊”到“銳利”的信號(hào)指紋提取如果你處理過(guò)一段音頻、一段振動(dòng)信號(hào)或者任何隨時(shí)間變化的非平穩(wěn)信號(hào)你大概率用過(guò)或者聽說(shuō)過(guò)短時(shí)傅里葉變換。它就像給信號(hào)戴上了一副“時(shí)間-頻率”眼鏡讓我們能同時(shí)看到信號(hào)在什么時(shí)候、有什么頻率成分。但用過(guò)的人都知道這副眼鏡有個(gè)“硬傷”——分辨率是固定的。一旦你選定了分析窗口的長(zhǎng)度時(shí)間分辨率和頻率分辨率就相互制約無(wú)法兼顧。這導(dǎo)致在時(shí)頻圖上一個(gè)尖銳的瞬時(shí)沖擊會(huì)變得“胖乎乎”一個(gè)純凈的單頻信號(hào)也會(huì)在時(shí)間軸上“拖泥帶水”我們看到的更像是一個(gè)模糊的“指紋”輪廓而非清晰的細(xì)節(jié)。這正是同步壓縮變換要解決的問(wèn)題。它不是一個(gè)全新的變換而是建立在STFT或小波變換結(jié)果之上的一種“后處理”技術(shù)。你可以把它想象成一個(gè)智能的“銳化”和“聚焦”算法。它通過(guò)分析STFT結(jié)果中每個(gè)點(diǎn)的“局部頻率”信息將那些能量模糊分布在周圍的點(diǎn)重新“壓縮”匯聚到其真實(shí)的瞬時(shí)頻率軌跡上。最終得到的時(shí)頻表示其頻率方向上的能量帶會(huì)變得非?！笆荨?、非?！颁J利”極大地提高了頻率分辨率同時(shí)幾乎不損失時(shí)間分辨率。這對(duì)于精確提取信號(hào)的瞬時(shí)頻率、分離緊密相鄰的頻率成分、識(shí)別微弱的瞬態(tài)沖擊具有革命性的意義。這個(gè)項(xiàng)目就是帶你從原理到代碼親手實(shí)現(xiàn)這個(gè)“信號(hào)指紋高清修復(fù)”的過(guò)程。無(wú)論你是從事機(jī)械故障診斷從振動(dòng)信號(hào)中定位軸承損傷頻率、語(yǔ)音信號(hào)處理分離共振峰、生物醫(yī)學(xué)工程分析心電、腦電信號(hào)的時(shí)變特性還是地球物理勘探只要你的信號(hào)是非平穩(wěn)的且你需要看清其頻率成分如何隨時(shí)間精細(xì)演化那么SST就是你工具箱里不可或缺的利器。接下來(lái)我將以一個(gè)包含兩個(gè)頻率非常接近的線性調(diào)頻信號(hào)為例帶你一步步拆解SST的原理并用Matlab代碼將其實(shí)現(xiàn)過(guò)程中我會(huì)分享那些官方文檔里不會(huì)寫的參數(shù)調(diào)優(yōu)心得和避坑指南。2. 核心原理拆解SST如何實(shí)現(xiàn)“時(shí)頻超分辨率”要理解同步壓縮變換我們必須先回到它的基礎(chǔ)——短時(shí)傅里葉變換。只有深刻理解了STFT的局限才能明白SST設(shè)計(jì)的精妙之處。2.1 短時(shí)傅里葉變換的“海森堡測(cè)不準(zhǔn)”困境STFT的核心思想很直觀用一個(gè)滑動(dòng)的、有限長(zhǎng)的窗函數(shù)比如漢明窗去截取信號(hào)對(duì)每一段加窗后的信號(hào)做傅里葉變換從而得到該時(shí)間段內(nèi)的頻譜。將所有時(shí)間段的頻譜排列起來(lái)就得到了時(shí)頻譜圖。用公式表示信號(hào)x(t)的STFT為STFT(t, ω) ∫ x(τ) g(τ - t) e^(-iω(τ-t)) dτ其中g(shù)(t)是窗函數(shù)。這里的關(guān)鍵在于窗函數(shù)g(t)的長(zhǎng)度。一個(gè)短的窗時(shí)間分辨率高能看清信號(hào)的快速變化但頻率分辨率低無(wú)法區(qū)分相近的頻率一個(gè)長(zhǎng)的窗頻率分辨率高但時(shí)間分辨率低會(huì)模糊掉瞬態(tài)事件。這就是時(shí)頻分析中的“測(cè)不準(zhǔn)原理”兩者無(wú)法同時(shí)達(dá)到最優(yōu)。在時(shí)頻譜圖上這表現(xiàn)為一個(gè)點(diǎn)擴(kuò)散函數(shù)。即使是一個(gè)理想的、頻率為ω0的純正弦信號(hào)其STFT在時(shí)頻平面上也不是一條無(wú)限細(xì)的直線而是一條沿著頻率軸有一定展寬的“能量帶”。這個(gè)展寬的寬度直接由窗函數(shù)的傅里葉變換?(ω)的寬度決定。換句話說(shuō)每個(gè)頻率成分的能量被“涂抹”開了模糊了真實(shí)的時(shí)頻結(jié)構(gòu)。2.2 同步壓縮變換的“再分配”哲學(xué)SST的核心思想源于D. Iatsenko等人提出的時(shí)頻再分配。它問(wèn)了一個(gè)關(guān)鍵問(wèn)題STFT時(shí)頻平面上某一點(diǎn)(t, ω)的能量真的應(yīng)該屬于頻率ω嗎對(duì)于純調(diào)頻信號(hào)STFT系數(shù)的相位包含了信號(hào)的瞬時(shí)頻率信息。SST通過(guò)計(jì)算一個(gè)稱為瞬時(shí)頻率估計(jì)的量來(lái)回答上述問(wèn)題。對(duì)于STFT這個(gè)估計(jì)值ω?(t, ω)可以通過(guò)STFT系數(shù)的相位隨時(shí)間的變化率即相位導(dǎo)數(shù)來(lái)計(jì)算ω?(t, ω) ω - Im{ (?STFT(t, ω)/?t) / STFT(t, ω) }這里Im表示取虛部。這個(gè)公式可能看起來(lái)有點(diǎn)復(fù)雜但其物理意義非常清晰它計(jì)算的是在時(shí)間t和頻率ω這個(gè)點(diǎn)上信號(hào)成分的局部振蕩頻率。如果這一點(diǎn)恰好位于信號(hào)的真實(shí)瞬時(shí)頻率軌跡上那么這個(gè)估計(jì)值就會(huì)接近真實(shí)頻率如果這一點(diǎn)只是由于窗函數(shù)展寬造成的能量泄漏那么這個(gè)估計(jì)值就會(huì)偏離當(dāng)前頻率ω。SST的“壓縮”動(dòng)作就基于此它遍歷STFT時(shí)頻平面的每一個(gè)點(diǎn)(t, ω)計(jì)算其瞬時(shí)頻率估計(jì)ω?(t, ω)然后將該點(diǎn)的能量|STFT(t, ω)|2或復(fù)數(shù)系數(shù)本身用于重構(gòu)從原來(lái)的位置(t, ω)“搬運(yùn)”或“壓縮”到新的位置(t, ω?(t, ω))上去。注意這里有一個(gè)非常重要的細(xì)節(jié)。我們搬運(yùn)的是復(fù)數(shù)系數(shù)STFT(t, ω)本身而不僅僅是能量。這是因?yàn)镾ST的一個(gè)巨大優(yōu)勢(shì)是完全可逆只要處理得當(dāng)可以從SST的結(jié)果中近乎完美地重構(gòu)原始信號(hào)。如果只搬運(yùn)能量模的平方就會(huì)丟失相位信息無(wú)法實(shí)現(xiàn)重構(gòu)。2.3 從連續(xù)公式到離散實(shí)現(xiàn)的關(guān)鍵步驟上面的公式是連續(xù)域的。在數(shù)字世界我們的信號(hào)是離散的STFT也是通過(guò)離散傅里葉變換計(jì)算的。因此實(shí)現(xiàn)SST需要解決幾個(gè)關(guān)鍵的離散化問(wèn)題相位導(dǎo)數(shù)的計(jì)算?STFT(t, ω)/?t需要離散近似。最常用且穩(wěn)定的方法是利用STFT在時(shí)間方向上的差分。假設(shè)我們的時(shí)間采樣索引是n那么可以用STFT[n1, k] - STFT[n-1, k]除以2Δt來(lái)近似時(shí)間導(dǎo)數(shù)中心差分法。這比前向或后向差分更精確。頻率軸的重新映射計(jì)算出的ω?[n, k]是一個(gè)連續(xù)的頻率值但我們的時(shí)頻圖輸出是一個(gè)離散的網(wǎng)格。我們需要將能量“分配”到離散的頻率倉(cāng)上。通常采用“投票”或“積累”的方式對(duì)于每個(gè)(n, k)找到ω?[n, k]對(duì)應(yīng)的最鄰近的頻率倉(cāng)索引k?然后將STFT[n, k]加到輸出矩陣的[n, k?]位置上。避免分母為零在計(jì)算ω?的公式中需要除以STFT(t, ω)。當(dāng)STFT系數(shù)非常小接近零時(shí)這會(huì)引入巨大的數(shù)值誤差。因此在實(shí)際計(jì)算中必須設(shè)定一個(gè)閾值只對(duì)那些幅度大于閾值的點(diǎn)進(jìn)行同步壓縮操作。低于閾值的點(diǎn)其能量通常被認(rèn)為是噪聲或數(shù)值誤差可以直接舍棄或保留在原位。窗函數(shù)的影響雖然SST能極大改善頻率聚焦性但其性能仍受初始STFT中窗函數(shù)選擇的間接影響。窗長(zhǎng)決定了初始時(shí)頻表示的“模糊”程度也影響了瞬時(shí)頻率估計(jì)的準(zhǔn)確性。通常需要選擇一個(gè)在時(shí)間和頻率上都有較好聚集性的窗如高斯窗。理解了這些我們就可以著手用Matlab搭建一個(gè)屬于自己的SST分析工具了。下面我將進(jìn)入最核心的實(shí)操環(huán)節(jié)。3. Matlab代碼實(shí)現(xiàn)一步步構(gòu)建SST分析儀我們將通過(guò)一個(gè)完整的Matlab腳本示例來(lái)演示如何生成測(cè)試信號(hào)計(jì)算STFT并實(shí)現(xiàn)同步壓縮變換。我會(huì)在代碼中插入大量注釋解釋每一步的目的和注意事項(xiàng)。3.1 測(cè)試信號(hào)生成與參數(shù)設(shè)置首先我們創(chuàng)建一個(gè)包含兩個(gè)成分的復(fù)雜信號(hào)以便直觀對(duì)比STFT和SST的效果。%% 1. 參數(shù)設(shè)置與測(cè)試信號(hào)生成 clear; close all; clc; % 信號(hào)參數(shù) fs 1000; % 采樣頻率 (Hz) T 2; % 信號(hào)時(shí)長(zhǎng) (秒) t 0:1/fs:T-1/fs; % 時(shí)間向量 N length(t); % 信號(hào)長(zhǎng)度 % 生成測(cè)試信號(hào)兩個(gè)線性調(diào)頻信號(hào) 一個(gè)瞬態(tài)沖擊 噪聲 % 成分1頻率從50Hz線性增加到150Hz f1 50 50*t/T; comp1 cos(2*pi * cumsum(f1)/fs); % 使用累積和來(lái)近似積分生成相位 % 成分2頻率從180Hz線性減少到80Hz與成分1在中間時(shí)段頻率接近 f2 180 - 100*t/T; comp2 0.8 * cos(2*pi * cumsum(f2)/fs); % 成分3在t1秒處的一個(gè)瞬態(tài)高斯包絡(luò)脈沖 transient exp(-100*(t-1).^2) .* cos(2*pi*250*t); % 成分4隨機(jī)噪聲 noise 0.1 * randn(size(t)); % 合成信號(hào) x comp1 comp2 transient noise; % 繪制原始信號(hào) figure(‘Position‘, [100, 100, 800, 400]); subplot(2,1,1); plot(t, x); xlabel(‘時(shí)間 (s)‘); ylabel(‘幅值‘); title(‘原始合成信號(hào)‘); grid on;實(shí)操心得1信號(hào)生成這里用cumsum(f)/fs來(lái)近似∫ f(t) dt對(duì)于線性調(diào)頻這類頻率變化平滑的信號(hào)在采樣率足夠高時(shí)是可行且簡(jiǎn)便的。對(duì)于精確的仿真可以考慮直接積分相位函數(shù)φ(t) 2π ∫ f(τ) dτ。3.2 短時(shí)傅里葉變換的實(shí)現(xiàn)接下來(lái)我們實(shí)現(xiàn)一個(gè)基礎(chǔ)的STFT函數(shù)。Matlab自帶的spectrogram函數(shù)雖然方便但為了更清晰地控制每一步并用于后續(xù)的SST我們選擇自己實(shí)現(xiàn)。%% 2. 短時(shí)傅里葉變換實(shí)現(xiàn) % STFT 參數(shù) win_len 128; % 窗長(zhǎng)度點(diǎn)數(shù)直接影響時(shí)頻分辨率權(quán)衡 hop 4; % 幀移點(diǎn)數(shù)hop越小時(shí)間軸越密計(jì)算量越大 win hamming(win_len, ‘periodic‘); % 使用漢明窗減少頻譜泄漏 nfft 256; % FFT點(diǎn)數(shù)通常 win_len用于頻率插值 % 計(jì)算STFT [STFT, f_stft, t_stft] my_stft(x, win, hop, nfft, fs); % 繪制STFT時(shí)頻譜能量譜密度 figure(‘Position‘, [100, 100, 1200, 500]); subplot(1,2,1); imagesc(t_stft, f_stft, 20*log10(abs(STFT))); % 轉(zhuǎn)換為dB尺度 axis xy; % 確保頻率軸方向正確 xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); title(‘傳統(tǒng)STFT時(shí)頻譜圖‘); colorbar; clim([-60, 0]); % 設(shè)置顏色范圍便于觀察這里調(diào)用了自定義函數(shù)my_stft。其實(shí)現(xiàn)如下重點(diǎn)在于邊界處理和矩陣運(yùn)算的效率function [STFT, f, t] my_stft(x, win, hop, nfft, fs) % 自定義STFT函數(shù)返回復(fù)數(shù)STFT矩陣、頻率向量和時(shí)間向量 L length(x); win_len length(win); % 計(jì)算幀數(shù) num_frames fix((L - win_len) / hop) 1; % 初始化STFT矩陣 (頻率倉(cāng) x 時(shí)間幀) STFT zeros(nfft, num_frames); % 逐幀處理 for i 0:num_frames-1 idx (i*hop) (1:win_len); segment x(idx) .* win; % 加窗 STFT(:, i1) fft(segment, nfft); % 做FFT end % 生成頻率和時(shí)間向量 f (0:nfft-1) * (fs / nfft); % 只取正頻率部分單邊譜如果需要的話 % STFT STFT(1:nfft/21, :); % f f(1:nfft/21); t (0:num_frames-1) * hop / fs; end注意事項(xiàng)1窗函數(shù)與重疊hamming窗的‘periodic‘選項(xiàng)適用于FFT能提供更好的頻譜特性。hop幀移通常設(shè)為窗長(zhǎng)的1/4到1/8在時(shí)間分辨率和計(jì)算量之間折衷。這里設(shè)為4時(shí)間分辨率非常高但計(jì)算量也大。實(shí)操心得2顯示動(dòng)態(tài)范圍時(shí)頻譜用dB尺度20*log10(abs(STFT))顯示是行業(yè)標(biāo)準(zhǔn)因?yàn)樗芡瑫r(shí)顯示很強(qiáng)和很弱的成分。clim用于統(tǒng)一顏色軸方便對(duì)比不同方法的結(jié)果。3.3 同步壓縮變換的核心算法實(shí)現(xiàn)這是整個(gè)項(xiàng)目的核心。我們將嚴(yán)格按照2.2和2.3節(jié)所述的原理來(lái)實(shí)現(xiàn)。%% 3. 同步壓縮變換核心實(shí)現(xiàn) function [SST, f_sst, t_sst] my_sst(STFT, t_stft, f_stft, fs, hop, thr) % 輸入 % STFT - 短時(shí)傅里葉變換結(jié)果矩陣 (頻率倉(cāng) x 時(shí)間幀) % t_stft, f_stft - STFT對(duì)應(yīng)的時(shí)間和頻率向量 % fs - 采樣率 % hop - STFT計(jì)算時(shí)的幀移點(diǎn)數(shù) % thr - 幅度閾值低于此值的STFT系數(shù)不參與壓縮 % 輸出 % SST - 同步壓縮變換結(jié)果矩陣 % f_sst, t_sst - 對(duì)應(yīng)的頻率和時(shí)間向量通常t_sst t_stft [n_freq, n_time] size(STFT); df f_stft(2) - f_stft(1); % 頻率分辨率 dt t_stft(2) - t_stft(1); % 時(shí)間分辨率理論上等于hop/fs % 初始化SST矩陣與STFT同尺寸 SST zeros(size(STFT)); % 為了避免復(fù)數(shù)運(yùn)算中的相位纏繞問(wèn)題我們使用STFT的導(dǎo)數(shù)來(lái)計(jì)算瞬時(shí)頻率 % 計(jì)算STFT對(duì)時(shí)間的偏導(dǎo)數(shù)采用中心差分 STFT_pad [zeros(n_freq,1), STFT, zeros(n_freq,1)]; % 在時(shí)間邊界填充零 dSTFT_dt (STFT_pad(:, 3:end) - STFT_pad(:, 1:end-2)) / (2*dt); % 中心差分 % 遍歷每個(gè)時(shí)頻點(diǎn) for ti 1:n_time for fi 1:n_freq STFT_coef STFT(fi, ti); coef_mag abs(STFT_coef); % 只處理幅度大于閾值的點(diǎn) if coef_mag thr % 計(jì)算瞬時(shí)頻率估計(jì) (公式 omega_hat omega - Im{(dSTFT/dt) / STFT}) if abs(STFT_coef) eps % 防止除以零 omega_inst f_stft(fi) - imag(dSTFT_dt(fi, ti) / STFT_coef) / (2*pi); % 上面除以2π是為了將角頻率(rad/s)轉(zhuǎn)換為普通頻率(Hz) else omega_inst f_stft(fi); end % 將瞬時(shí)頻率映射到最近的頻率倉(cāng)索引 k_hat round(omega_inst / df) 1; % 1 因?yàn)镸atlab索引從1開始 % 確保映射后的索引在有效范圍內(nèi) if k_hat 1 k_hat n_freq % 將當(dāng)前STFT系數(shù)累加到SST矩陣的對(duì)應(yīng)位置 % 注意這里是復(fù)數(shù)累加以保留重構(gòu)能力 SST(k_hat, ti) SST(k_hat, ti) STFT_coef; end end end end f_sst f_stft; t_sst t_stft; end在主腳本中調(diào)用這個(gè)函數(shù)% 設(shè)置SST參數(shù) thr max(abs(STFT(:))) * 0.01; % 閾值設(shè)為STFT最大幅值的1% % 計(jì)算SST [SST, f_sst, t_sst] my_sst(STFT, t_stft, f_stft, fs, hop, thr); % 繪制SST時(shí)頻譜 subplot(1,2,2); imagesc(t_sst, f_sst, 20*log10(abs(SST)eps)); % 加eps避免log10(0) axis xy; xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); title(‘同步壓縮變換時(shí)頻譜圖‘); colorbar; clim([-60, 0]); % 使用與STFT相同的顏色范圍3.4 結(jié)果對(duì)比分析與解讀運(yùn)行上述代碼后你會(huì)得到并排的兩幅時(shí)頻譜圖。對(duì)比它們你可以立即發(fā)現(xiàn)SST的魔力頻率聚焦性在STFT圖中兩個(gè)線性調(diào)頻信號(hào)是兩條較粗的、有一定寬度的“能量帶”。尤其是在時(shí)間中部約1秒處當(dāng)兩個(gè)信號(hào)的頻率非常接近時(shí)它們的能量帶會(huì)重疊、模糊在一起難以清晰分辨。而在SST圖中這兩條軌跡變成了極其銳利的細(xì)線即使它們靠得很近也能被清晰地區(qū)分開。這就是頻率分辨率的大幅提升。瞬態(tài)成分表征對(duì)于t1秒處的瞬時(shí)脈沖中心頻率250Hz在STFT圖中由于窗函數(shù)的限制它在時(shí)間軸上被“拉長(zhǎng)”了在頻率軸上也有一定的展寬看起來(lái)像一個(gè)“斑點(diǎn)”。在SST圖中這個(gè)脈沖在時(shí)間上依然被精確定位沒(méi)有因?yàn)閴嚎s而模糊時(shí)間信息同時(shí)在頻率上也變得更加集中更接近一個(gè)理想的時(shí)頻點(diǎn)。噪聲抑制觀察背景噪聲圖像中均勻分布的藍(lán)色背景。在SST圖中背景噪聲的強(qiáng)度似乎有所降低或變得更加“稀疏”。這是因?yàn)樵肼暤腟TFT系數(shù)相位是隨機(jī)的其計(jì)算出的瞬時(shí)頻率估計(jì)ω?也會(huì)非常隨機(jī)導(dǎo)致在再分配過(guò)程中能量被分散地映射到各個(gè)頻率倉(cāng)而不會(huì)像真實(shí)信號(hào)那樣集中到一條線上。因此在SST結(jié)果中信號(hào)的能量更加集中而噪聲的能量相對(duì)更加分散這在一定程度上提升了時(shí)頻譜的信噪比。注意事項(xiàng)2閾值選擇閾值thr的選擇至關(guān)重要。設(shè)得太高會(huì)丟失微弱信號(hào)設(shè)得太低會(huì)讓大量噪聲點(diǎn)參與壓縮不僅增加計(jì)算量還可能因噪聲點(diǎn)的隨機(jī)瞬時(shí)頻率估計(jì)而污染結(jié)果。通常建議設(shè)為STFT最大幅值的0.5%到5%之間需要根據(jù)具體信號(hào)的信噪比進(jìn)行微調(diào)。4. 關(guān)鍵參數(shù)影響與調(diào)優(yōu)指南SST的效果并非一勞永逸它嚴(yán)重依賴于初始STFT的參數(shù)設(shè)置。下面我們通過(guò)一個(gè)參數(shù)研究來(lái)理解這些影響。4.1 窗長(zhǎng)時(shí)頻分辨率的“起跑線”窗長(zhǎng)是STFT最核心的參數(shù)也是SST效果的基石。%% 4. 參數(shù)影響分析窗長(zhǎng) win_lens [64, 128, 256]; % 測(cè)試三種窗長(zhǎng) figure(‘Position‘, [100, 100, 1200, 900]); for i 1:length(win_lens) win_len win_lens(i); win hamming(win_len, ‘periodic‘); hop max(4, floor(win_len/16)); % 幀移隨窗長(zhǎng)適度增加 nfft 2^nextpow2(win_len*2); [STFT_temp, f_temp, t_temp] my_stft(x, win, hop, nfft, fs); thr_temp max(abs(STFT_temp(:))) * 0.01; [SST_temp, ~, ~] my_sst(STFT_temp, t_temp, f_temp, fs, hop, thr_temp); % 繪制STFT subplot(3, 2, (i-1)*21); imagesc(t_temp, f_temp, 20*log10(abs(STFT_temp))); axis xy; title([‘STFT - 窗長(zhǎng) ‘, num2str(win_len)]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); % 繪制SST subplot(3, 2, (i-1)*22); imagesc(t_temp, f_temp, 20*log10(abs(SST_temp)eps)); axis xy; title([‘SST - 窗長(zhǎng) ‘, num2str(win_len)]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); end結(jié)果分析短窗64點(diǎn)STFT的時(shí)間分辨率很高兩個(gè)調(diào)頻信號(hào)的軌跡在時(shí)間起止點(diǎn)清晰但頻率分辨率極差軌跡非常粗幾乎無(wú)法分辨中間接近的部分。SST試圖壓縮但“原料”太粗糙效果提升有限軌跡仍然較寬且可能出現(xiàn)斷點(diǎn)。中窗128點(diǎn)我們的初始選擇STFT的時(shí)頻權(quán)衡相對(duì)均衡。SST效果顯著軌跡銳利分離清晰。長(zhǎng)窗256點(diǎn)STFT的頻率分辨率很高兩條軌跡在頻率上本身已較清晰但時(shí)間分辨率下降瞬態(tài)脈沖被嚴(yán)重拉長(zhǎng)、模糊。SST能進(jìn)一步銳化頻率軌跡但無(wú)法修復(fù)損失的時(shí)間分辨率脈沖在SST中依然是被拉寬的。核心結(jié)論SST能顯著提升頻率分辨率但無(wú)法突破STFT初始時(shí)間分辨率的理論上限。它主要修復(fù)由窗函數(shù)引起的頻率方向上的能量擴(kuò)散。因此選擇窗長(zhǎng)的首要原則是確保STFT能捕捉到信號(hào)中最快的時(shí)間變化即所需的時(shí)間分辨率。在這個(gè)基礎(chǔ)上SST來(lái)優(yōu)化頻率分辨率。4.2 閾值信號(hào)與噪聲的“分水嶺”閾值決定了哪些STFT系數(shù)參與再分配。%% 5. 參數(shù)影響分析閾值 win_len 128; win hamming(win_len, ‘periodic‘); hop 4; nfft 256; [STFT_base, f_base, t_base] my_stft(x, win, hop, nfft, fs); thresholds [0.001, 0.01, 0.05]; % 相對(duì)于最大幅值的比例 figure(‘Position‘, [100, 100, 1200, 400]); for i 1:length(thresholds) thr max(abs(STFT_base(:))) * thresholds(i); [SST_temp, ~, ~] my_sst(STFT_base, t_base, f_base, fs, hop, thr); subplot(1, 3, i); imagesc(t_base, f_base, 20*log10(abs(SST_temp)eps)); axis xy; title([‘SST - 閾值 ‘, num2str(thresholds(i)*100), ‘%‘]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); colorbar; end結(jié)果分析低閾值0.1%幾乎所有點(diǎn)都參與壓縮包括大量噪聲點(diǎn)。結(jié)果圖中背景噪聲也呈現(xiàn)出一些虛假的、稀疏的“點(diǎn)狀”或“短線段”結(jié)構(gòu)這是因?yàn)樵肼暤碾S機(jī)相位導(dǎo)致了隨機(jī)的瞬時(shí)頻率估計(jì)。整體圖像可能看起來(lái)有點(diǎn)“臟”。適中閾值1%推薦起點(diǎn)大部分噪聲被過(guò)濾掉信號(hào)軌跡清晰銳利背景干凈。這是通常的起始選擇。高閾值5%只有能量最強(qiáng)的信號(hào)核心部分參與壓縮?？赡軐?dǎo)致微弱信號(hào)成分如我們信號(hào)中幅度為0.8的第二個(gè)成分的某些部分丟失軌跡出現(xiàn)不連續(xù)。同時(shí)瞬態(tài)脈沖的邊緣部分可能被舍棄。調(diào)優(yōu)建議從1%的閾值開始。如果發(fā)現(xiàn)微弱信號(hào)丟失適當(dāng)降低閾值如0.5%。如果背景噪聲干擾嚴(yán)重呈現(xiàn)虛假結(jié)構(gòu)則適當(dāng)提高閾值如2%?？梢杂^察SST結(jié)果中噪聲基底的特征來(lái)判斷。4.3 頻率軸插值與迭代SST基礎(chǔ)的SST將能量壓縮到離散的頻率網(wǎng)格上這可能導(dǎo)致“量化誤差”。更高級(jí)的實(shí)現(xiàn)可以采用以下技巧頻率軸插值在計(jì)算k_hat時(shí)不使用簡(jiǎn)單的round取整而是將能量按一定權(quán)重分配到相鄰的兩個(gè)頻率倉(cāng)上如線性插值這可以減輕因離散化造成的“柵欄效應(yīng)”使結(jié)果更平滑。這通常能帶來(lái)輕微的視覺改善。迭代SST將第一次SST的結(jié)果作為輸入再次進(jìn)行同步壓縮。理論上可以進(jìn)一步聚焦。但實(shí)踐中一次壓縮通常已能達(dá)到很好效果多次迭代可能引入偽影且計(jì)算成本翻倍。除非對(duì)時(shí)頻脊線提取有極高要求否則不建議常規(guī)使用。5. 常見問(wèn)題、排查技巧與進(jìn)階應(yīng)用在實(shí)際使用自編SST代碼時(shí)你可能會(huì)遇到以下典型問(wèn)題。5.1 時(shí)頻譜出現(xiàn)水平條紋或斷裂現(xiàn)象SST結(jié)果圖中本應(yīng)連續(xù)的信號(hào)軌跡出現(xiàn)明顯的水平斷裂帶或者在整個(gè)時(shí)間軸上出現(xiàn)均勻的水平條紋?？赡茉蚺c排查相位導(dǎo)數(shù)計(jì)算不準(zhǔn)確這是最常見的原因。確保計(jì)算dSTFT_dt時(shí)使用的是中心差分法并且時(shí)間步長(zhǎng)dt計(jì)算正確dt hop / fs。避免使用前向或后向差分它們?cè)谶吔缣幷`差大且整體精度低。邊界效應(yīng)我們的my_sst函數(shù)在計(jì)算時(shí)間導(dǎo)數(shù)時(shí)通過(guò)補(bǔ)零來(lái)近似中心差分但信號(hào)兩端的導(dǎo)數(shù)計(jì)算本身就不準(zhǔn)確。這會(huì)導(dǎo)致時(shí)間第一幀和最后一幀的瞬時(shí)頻率估計(jì)錯(cuò)誤從而產(chǎn)生邊界處的畸變。一種改進(jìn)方法是使用更復(fù)雜的邊界處理或者簡(jiǎn)單地在分析時(shí)舍棄頭尾幾幀。閾值過(guò)高過(guò)高的閾值會(huì)剔除掉構(gòu)成連續(xù)軌跡所必需的、幅度稍低的點(diǎn)導(dǎo)致軌跡斷裂。嘗試降低閾值。5.2 重構(gòu)信號(hào)誤差大現(xiàn)象使用SST的系數(shù)進(jìn)行信號(hào)重構(gòu)逆變換時(shí)重構(gòu)信號(hào)與原始信號(hào)差異顯著。可能原因與排查能量歸一化問(wèn)題SST是一個(gè)線性重分配過(guò)程但簡(jiǎn)單的“投票式”累加會(huì)改變系數(shù)的總能量。嚴(yán)格的可逆SST需要滿足保范數(shù)條件即在再分配過(guò)程中每個(gè)源點(diǎn)貢獻(xiàn)的能量權(quán)重需要精心設(shè)計(jì)使得整個(gè)變換是等距的。我們的基礎(chǔ)實(shí)現(xiàn)未做此處理因此逆變換不完美。若需精確重構(gòu)需查閱文獻(xiàn)實(shí)現(xiàn)“二階”或“可逆”SST。僅使用了SST的模如果只壓縮了能量abs(STFT)^2而丟棄了相位信息則絕對(duì)無(wú)法重構(gòu)。我們的代碼壓縮的是復(fù)數(shù)STFT保留了重構(gòu)的可能性。數(shù)值誤差累積相位導(dǎo)數(shù)的計(jì)算涉及除法對(duì)數(shù)值誤差敏感。確保使用雙精度計(jì)算并對(duì)極小分母進(jìn)行保護(hù)代碼中的eps檢查。5.3 對(duì)多分量信號(hào)中交叉軌跡的處理現(xiàn)象當(dāng)兩個(gè)信號(hào)的時(shí)頻軌跡在某個(gè)時(shí)間點(diǎn)交叉時(shí)SST結(jié)果在交叉點(diǎn)附近可能出現(xiàn)模糊或畸變。原因與對(duì)策這是SST以及大多數(shù)時(shí)頻后處理方法的固有挑戰(zhàn)。在交叉點(diǎn)信號(hào)的局部相位特性變得復(fù)雜瞬時(shí)頻率估計(jì)可能失效。對(duì)于交叉軌跡嘗試更短的窗短窗雖然初始頻率分辨率差但能更好地分離時(shí)間上快速變化的成分可能使交叉點(diǎn)的影響區(qū)域變小。使用方向性SST有研究提出在交叉區(qū)域根據(jù)信號(hào)分量方向進(jìn)行選擇性壓縮的算法但這非常復(fù)雜。接受局限對(duì)于高度非平穩(wěn)、分量交叉的信號(hào)需要認(rèn)識(shí)到時(shí)頻分析工具的局限性結(jié)合其他方法如經(jīng)驗(yàn)?zāi)B(tài)分解EMD先進(jìn)行信號(hào)分離再對(duì)單分量做SST。5.4 在強(qiáng)噪聲環(huán)境下的表現(xiàn)現(xiàn)象信號(hào)信噪比很低時(shí)SST可能無(wú)法清晰提取出信號(hào)軌跡甚至可能因噪聲而產(chǎn)生虛假結(jié)構(gòu)。優(yōu)化策略前置去噪在SST之前先對(duì)信號(hào)進(jìn)行濾波或小波去噪預(yù)處理。閾值調(diào)優(yōu)提高閾值只壓縮能量顯著高于噪聲基底的點(diǎn)。結(jié)合魯棒性估計(jì)使用更魯棒的瞬時(shí)頻率估計(jì)方法例如基于時(shí)頻分布如Wigner-Ville分布的重分配方法但計(jì)算量更大。多次平均如果條件允許對(duì)同一現(xiàn)象進(jìn)行多次測(cè)量對(duì)SST幅值譜進(jìn)行平均可以抑制隨機(jī)噪聲。5.5 計(jì)算效率優(yōu)化我們的雙循環(huán)實(shí)現(xiàn)直觀但較慢。對(duì)于長(zhǎng)信號(hào)或?qū)崟r(shí)處理可以考慮以下優(yōu)化向量化利用Matlab的矩陣運(yùn)算避免雙重循環(huán)?？梢酝瑫r(shí)計(jì)算所有點(diǎn)的瞬時(shí)頻率估計(jì)需處理除以零問(wèn)題并使用accumarray函數(shù)進(jìn)行高效的“投票”累加。這能帶來(lái)數(shù)量級(jí)的速度提升。使用C/MEX編碼將核心循環(huán)用C語(yǔ)言編寫并通過(guò)MEX接口調(diào)用適用于對(duì)性能要求極高的場(chǎng)合。利用GPUSST的并行性很好可以使用Matlab的Parallel Computing Toolbox或直接使用CUDA進(jìn)行GPU加速。最后分享一個(gè)我個(gè)人的深刻體會(huì)同步壓縮變換是一個(gè)極其強(qiáng)大的工具但它不是“銀彈”。它完美解決了單分量調(diào)頻信號(hào)在時(shí)頻譜上頻率擴(kuò)散的問(wèn)題。理解它的原理基于相位導(dǎo)數(shù)的再分配比單純調(diào)用一個(gè)函數(shù)更重要。這能幫助你在面對(duì)復(fù)雜信號(hào)時(shí)正確解讀SST的結(jié)果判斷哪些是真實(shí)的信號(hào)特征哪些可能是方法局限帶來(lái)的偽影。從STFT到SST就像是從一幅模糊的素描到一張清晰的高清照片而掌握拍攝參數(shù)設(shè)置和后期處理SST算法的技巧才能讓你成為真正的“信號(hào)攝影師”。 同步壓縮變換：從STFT到時(shí)頻超分辨率的信號(hào)指紋高清修復(fù)
1. 項(xiàng)目概述從“模糊”到“銳利”的信號(hào)指紋提取如果你處理過(guò)一段音頻、一段振動(dòng)信號(hào)或者任何隨時(shí)間變化的非平穩(wěn)信號(hào)你大概率用過(guò)或者聽說(shuō)過(guò)短時(shí)傅里葉變換。它就像給信號(hào)戴上了一副“時(shí)間-頻率”眼鏡讓我們能同時(shí)看到信號(hào)在什么時(shí)候、有什么頻率成分。但用過(guò)的人都知道這副眼鏡有個(gè)“硬傷”——分辨率是固定的。一旦你選定了分析窗口的長(zhǎng)度時(shí)間分辨率和頻率分辨率就相互制約無(wú)法兼顧。這導(dǎo)致在時(shí)頻圖上一個(gè)尖銳的瞬時(shí)沖擊會(huì)變得“胖乎乎”一個(gè)純凈的單頻信號(hào)也會(huì)在時(shí)間軸上“拖泥帶水”我們看到的更像是一個(gè)模糊的“指紋”輪廓而非清晰的細(xì)節(jié)。這正是同步壓縮變換要解決的問(wèn)題。它不是一個(gè)全新的變換而是建立在STFT或小波變換結(jié)果之上的一種“后處理”技術(shù)。你可以把它想象成一個(gè)智能的“銳化”和“聚焦”算法。它通過(guò)分析STFT結(jié)果中每個(gè)點(diǎn)的“局部頻率”信息將那些能量模糊分布在周圍的點(diǎn)重新“壓縮”匯聚到其真實(shí)的瞬時(shí)頻率軌跡上。最終得到的時(shí)頻表示其頻率方向上的能量帶會(huì)變得非?！笆荨?、非?！颁J利”極大地提高了頻率分辨率同時(shí)幾乎不損失時(shí)間分辨率。這對(duì)于精確提取信號(hào)的瞬時(shí)頻率、分離緊密相鄰的頻率成分、識(shí)別微弱的瞬態(tài)沖擊具有革命性的意義。這個(gè)項(xiàng)目就是帶你從原理到代碼親手實(shí)現(xiàn)這個(gè)“信號(hào)指紋高清修復(fù)”的過(guò)程。無(wú)論你是從事機(jī)械故障診斷從振動(dòng)信號(hào)中定位軸承損傷頻率、語(yǔ)音信號(hào)處理分離共振峰、生物醫(yī)學(xué)工程分析心電、腦電信號(hào)的時(shí)變特性還是地球物理勘探只要你的信號(hào)是非平穩(wěn)的且你需要看清其頻率成分如何隨時(shí)間精細(xì)演化那么SST就是你工具箱里不可或缺的利器。接下來(lái)我將以一個(gè)包含兩個(gè)頻率非常接近的線性調(diào)頻信號(hào)為例帶你一步步拆解SST的原理并用Matlab代碼將其實(shí)現(xiàn)過(guò)程中我會(huì)分享那些官方文檔里不會(huì)寫的參數(shù)調(diào)優(yōu)心得和避坑指南。2. 核心原理拆解SST如何實(shí)現(xiàn)“時(shí)頻超分辨率”要理解同步壓縮變換我們必須先回到它的基礎(chǔ)——短時(shí)傅里葉變換。只有深刻理解了STFT的局限才能明白SST設(shè)計(jì)的精妙之處。2.1 短時(shí)傅里葉變換的“海森堡測(cè)不準(zhǔn)”困境STFT的核心思想很直觀用一個(gè)滑動(dòng)的、有限長(zhǎng)的窗函數(shù)比如漢明窗去截取信號(hào)對(duì)每一段加窗后的信號(hào)做傅里葉變換從而得到該時(shí)間段內(nèi)的頻譜。將所有時(shí)間段的頻譜排列起來(lái)就得到了時(shí)頻譜圖。用公式表示信號(hào)x(t)的STFT為STFT(t, ω) ∫ x(τ) g(τ - t) e^(-iω(τ-t)) dτ其中g(shù)(t)是窗函數(shù)。這里的關(guān)鍵在于窗函數(shù)g(t)的長(zhǎng)度。一個(gè)短的窗時(shí)間分辨率高能看清信號(hào)的快速變化但頻率分辨率低無(wú)法區(qū)分相近的頻率一個(gè)長(zhǎng)的窗頻率分辨率高但時(shí)間分辨率低會(huì)模糊掉瞬態(tài)事件。這就是時(shí)頻分析中的“測(cè)不準(zhǔn)原理”兩者無(wú)法同時(shí)達(dá)到最優(yōu)。在時(shí)頻譜圖上這表現(xiàn)為一個(gè)點(diǎn)擴(kuò)散函數(shù)。即使是一個(gè)理想的、頻率為ω0的純正弦信號(hào)其STFT在時(shí)頻平面上也不是一條無(wú)限細(xì)的直線而是一條沿著頻率軸有一定展寬的“能量帶”。這個(gè)展寬的寬度直接由窗函數(shù)的傅里葉變換?(ω)的寬度決定。換句話說(shuō)每個(gè)頻率成分的能量被“涂抹”開了模糊了真實(shí)的時(shí)頻結(jié)構(gòu)。2.2 同步壓縮變換的“再分配”哲學(xué)SST的核心思想源于D. Iatsenko等人提出的時(shí)頻再分配。它問(wèn)了一個(gè)關(guān)鍵問(wèn)題STFT時(shí)頻平面上某一點(diǎn)(t, ω)的能量真的應(yīng)該屬于頻率ω嗎對(duì)于純調(diào)頻信號(hào)STFT系數(shù)的相位包含了信號(hào)的瞬時(shí)頻率信息。SST通過(guò)計(jì)算一個(gè)稱為瞬時(shí)頻率估計(jì)的量來(lái)回答上述問(wèn)題。對(duì)于STFT這個(gè)估計(jì)值ω?(t, ω)可以通過(guò)STFT系數(shù)的相位隨時(shí)間的變化率即相位導(dǎo)數(shù)來(lái)計(jì)算ω?(t, ω) ω - Im{ (?STFT(t, ω)/?t) / STFT(t, ω) }這里Im表示取虛部。這個(gè)公式可能看起來(lái)有點(diǎn)復(fù)雜但其物理意義非常清晰它計(jì)算的是在時(shí)間t和頻率ω這個(gè)點(diǎn)上信號(hào)成分的局部振蕩頻率。如果這一點(diǎn)恰好位于信號(hào)的真實(shí)瞬時(shí)頻率軌跡上那么這個(gè)估計(jì)值就會(huì)接近真實(shí)頻率如果這一點(diǎn)只是由于窗函數(shù)展寬造成的能量泄漏那么這個(gè)估計(jì)值就會(huì)偏離當(dāng)前頻率ω。SST的“壓縮”動(dòng)作就基于此它遍歷STFT時(shí)頻平面的每一個(gè)點(diǎn)(t, ω)計(jì)算其瞬時(shí)頻率估計(jì)ω?(t, ω)然后將該點(diǎn)的能量|STFT(t, ω)|2或復(fù)數(shù)系數(shù)本身用于重構(gòu)從原來(lái)的位置(t, ω)“搬運(yùn)”或“壓縮”到新的位置(t, ω?(t, ω))上去。注意這里有一個(gè)非常重要的細(xì)節(jié)。我們搬運(yùn)的是復(fù)數(shù)系數(shù)STFT(t, ω)本身而不僅僅是能量。這是因?yàn)镾ST的一個(gè)巨大優(yōu)勢(shì)是完全可逆只要處理得當(dāng)可以從SST的結(jié)果中近乎完美地重構(gòu)原始信號(hào)。如果只搬運(yùn)能量模的平方就會(huì)丟失相位信息無(wú)法實(shí)現(xiàn)重構(gòu)。2.3 從連續(xù)公式到離散實(shí)現(xiàn)的關(guān)鍵步驟上面的公式是連續(xù)域的。在數(shù)字世界我們的信號(hào)是離散的STFT也是通過(guò)離散傅里葉變換計(jì)算的。因此實(shí)現(xiàn)SST需要解決幾個(gè)關(guān)鍵的離散化問(wèn)題相位導(dǎo)數(shù)的計(jì)算?STFT(t, ω)/?t需要離散近似。最常用且穩(wěn)定的方法是利用STFT在時(shí)間方向上的差分。假設(shè)我們的時(shí)間采樣索引是n那么可以用STFT[n1, k] - STFT[n-1, k]除以2Δt來(lái)近似時(shí)間導(dǎo)數(shù)中心差分法。這比前向或后向差分更精確。頻率軸的重新映射計(jì)算出的ω?[n, k]是一個(gè)連續(xù)的頻率值但我們的時(shí)頻圖輸出是一個(gè)離散的網(wǎng)格。我們需要將能量“分配”到離散的頻率倉(cāng)上。通常采用“投票”或“積累”的方式對(duì)于每個(gè)(n, k)找到ω?[n, k]對(duì)應(yīng)的最鄰近的頻率倉(cāng)索引k?然后將STFT[n, k]加到輸出矩陣的[n, k?]位置上。避免分母為零在計(jì)算ω?的公式中需要除以STFT(t, ω)。當(dāng)STFT系數(shù)非常小接近零時(shí)這會(huì)引入巨大的數(shù)值誤差。因此在實(shí)際計(jì)算中必須設(shè)定一個(gè)閾值只對(duì)那些幅度大于閾值的點(diǎn)進(jìn)行同步壓縮操作。低于閾值的點(diǎn)其能量通常被認(rèn)為是噪聲或數(shù)值誤差可以直接舍棄或保留在原位。窗函數(shù)的影響雖然SST能極大改善頻率聚焦性但其性能仍受初始STFT中窗函數(shù)選擇的間接影響。窗長(zhǎng)決定了初始時(shí)頻表示的“模糊”程度也影響了瞬時(shí)頻率估計(jì)的準(zhǔn)確性。通常需要選擇一個(gè)在時(shí)間和頻率上都有較好聚集性的窗如高斯窗。理解了這些我們就可以著手用Matlab搭建一個(gè)屬于自己的SST分析工具了。下面我將進(jìn)入最核心的實(shí)操環(huán)節(jié)。3. Matlab代碼實(shí)現(xiàn)一步步構(gòu)建SST分析儀我們將通過(guò)一個(gè)完整的Matlab腳本示例來(lái)演示如何生成測(cè)試信號(hào)計(jì)算STFT并實(shí)現(xiàn)同步壓縮變換。我會(huì)在代碼中插入大量注釋解釋每一步的目的和注意事項(xiàng)。3.1 測(cè)試信號(hào)生成與參數(shù)設(shè)置首先我們創(chuàng)建一個(gè)包含兩個(gè)成分的復(fù)雜信號(hào)以便直觀對(duì)比STFT和SST的效果。%% 1. 參數(shù)設(shè)置與測(cè)試信號(hào)生成 clear; close all; clc; % 信號(hào)參數(shù) fs 1000; % 采樣頻率 (Hz) T 2; % 信號(hào)時(shí)長(zhǎng) (秒) t 0:1/fs:T-1/fs; % 時(shí)間向量 N length(t); % 信號(hào)長(zhǎng)度 % 生成測(cè)試信號(hào)兩個(gè)線性調(diào)頻信號(hào) 一個(gè)瞬態(tài)沖擊 噪聲 % 成分1頻率從50Hz線性增加到150Hz f1 50 50*t/T; comp1 cos(2*pi * cumsum(f1)/fs); % 使用累積和來(lái)近似積分生成相位 % 成分2頻率從180Hz線性減少到80Hz與成分1在中間時(shí)段頻率接近 f2 180 - 100*t/T; comp2 0.8 * cos(2*pi * cumsum(f2)/fs); % 成分3在t1秒處的一個(gè)瞬態(tài)高斯包絡(luò)脈沖 transient exp(-100*(t-1).^2) .* cos(2*pi*250*t); % 成分4隨機(jī)噪聲 noise 0.1 * randn(size(t)); % 合成信號(hào) x comp1 comp2 transient noise; % 繪制原始信號(hào) figure(‘Position‘, [100, 100, 800, 400]); subplot(2,1,1); plot(t, x); xlabel(‘時(shí)間 (s)‘); ylabel(‘幅值‘); title(‘原始合成信號(hào)‘); grid on;實(shí)操心得1信號(hào)生成這里用cumsum(f)/fs來(lái)近似∫ f(t) dt對(duì)于線性調(diào)頻這類頻率變化平滑的信號(hào)在采樣率足夠高時(shí)是可行且簡(jiǎn)便的。對(duì)于精確的仿真可以考慮直接積分相位函數(shù)φ(t) 2π ∫ f(τ) dτ。3.2 短時(shí)傅里葉變換的實(shí)現(xiàn)接下來(lái)我們實(shí)現(xiàn)一個(gè)基礎(chǔ)的STFT函數(shù)。Matlab自帶的spectrogram函數(shù)雖然方便但為了更清晰地控制每一步并用于后續(xù)的SST我們選擇自己實(shí)現(xiàn)。%% 2. 短時(shí)傅里葉變換實(shí)現(xiàn) % STFT 參數(shù) win_len 128; % 窗長(zhǎng)度點(diǎn)數(shù)直接影響時(shí)頻分辨率權(quán)衡 hop 4; % 幀移點(diǎn)數(shù)hop越小時(shí)間軸越密計(jì)算量越大 win hamming(win_len, ‘periodic‘); % 使用漢明窗減少頻譜泄漏 nfft 256; % FFT點(diǎn)數(shù)通常 win_len用于頻率插值 % 計(jì)算STFT [STFT, f_stft, t_stft] my_stft(x, win, hop, nfft, fs); % 繪制STFT時(shí)頻譜能量譜密度 figure(‘Position‘, [100, 100, 1200, 500]); subplot(1,2,1); imagesc(t_stft, f_stft, 20*log10(abs(STFT))); % 轉(zhuǎn)換為dB尺度 axis xy; % 確保頻率軸方向正確 xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); title(‘傳統(tǒng)STFT時(shí)頻譜圖‘); colorbar; clim([-60, 0]); % 設(shè)置顏色范圍便于觀察這里調(diào)用了自定義函數(shù)my_stft。其實(shí)現(xiàn)如下重點(diǎn)在于邊界處理和矩陣運(yùn)算的效率function [STFT, f, t] my_stft(x, win, hop, nfft, fs) % 自定義STFT函數(shù)返回復(fù)數(shù)STFT矩陣、頻率向量和時(shí)間向量 L length(x); win_len length(win); % 計(jì)算幀數(shù) num_frames fix((L - win_len) / hop) 1; % 初始化STFT矩陣 (頻率倉(cāng) x 時(shí)間幀) STFT zeros(nfft, num_frames); % 逐幀處理 for i 0:num_frames-1 idx (i*hop) (1:win_len); segment x(idx) .* win; % 加窗 STFT(:, i1) fft(segment, nfft); % 做FFT end % 生成頻率和時(shí)間向量 f (0:nfft-1) * (fs / nfft); % 只取正頻率部分單邊譜如果需要的話 % STFT STFT(1:nfft/21, :); % f f(1:nfft/21); t (0:num_frames-1) * hop / fs; end注意事項(xiàng)1窗函數(shù)與重疊hamming窗的‘periodic‘選項(xiàng)適用于FFT能提供更好的頻譜特性。hop幀移通常設(shè)為窗長(zhǎng)的1/4到1/8在時(shí)間分辨率和計(jì)算量之間折衷。這里設(shè)為4時(shí)間分辨率非常高但計(jì)算量也大。實(shí)操心得2顯示動(dòng)態(tài)范圍時(shí)頻譜用dB尺度20*log10(abs(STFT))顯示是行業(yè)標(biāo)準(zhǔn)因?yàn)樗芡瑫r(shí)顯示很強(qiáng)和很弱的成分。clim用于統(tǒng)一顏色軸方便對(duì)比不同方法的結(jié)果。3.3 同步壓縮變換的核心算法實(shí)現(xiàn)這是整個(gè)項(xiàng)目的核心。我們將嚴(yán)格按照2.2和2.3節(jié)所述的原理來(lái)實(shí)現(xiàn)。%% 3. 同步壓縮變換核心實(shí)現(xiàn) function [SST, f_sst, t_sst] my_sst(STFT, t_stft, f_stft, fs, hop, thr) % 輸入 % STFT - 短時(shí)傅里葉變換結(jié)果矩陣 (頻率倉(cāng) x 時(shí)間幀) % t_stft, f_stft - STFT對(duì)應(yīng)的時(shí)間和頻率向量 % fs - 采樣率 % hop - STFT計(jì)算時(shí)的幀移點(diǎn)數(shù) % thr - 幅度閾值低于此值的STFT系數(shù)不參與壓縮 % 輸出 % SST - 同步壓縮變換結(jié)果矩陣 % f_sst, t_sst - 對(duì)應(yīng)的頻率和時(shí)間向量通常t_sst t_stft [n_freq, n_time] size(STFT); df f_stft(2) - f_stft(1); % 頻率分辨率 dt t_stft(2) - t_stft(1); % 時(shí)間分辨率理論上等于hop/fs % 初始化SST矩陣與STFT同尺寸 SST zeros(size(STFT)); % 為了避免復(fù)數(shù)運(yùn)算中的相位纏繞問(wèn)題我們使用STFT的導(dǎo)數(shù)來(lái)計(jì)算瞬時(shí)頻率 % 計(jì)算STFT對(duì)時(shí)間的偏導(dǎo)數(shù)采用中心差分 STFT_pad [zeros(n_freq,1), STFT, zeros(n_freq,1)]; % 在時(shí)間邊界填充零 dSTFT_dt (STFT_pad(:, 3:end) - STFT_pad(:, 1:end-2)) / (2*dt); % 中心差分 % 遍歷每個(gè)時(shí)頻點(diǎn) for ti 1:n_time for fi 1:n_freq STFT_coef STFT(fi, ti); coef_mag abs(STFT_coef); % 只處理幅度大于閾值的點(diǎn) if coef_mag thr % 計(jì)算瞬時(shí)頻率估計(jì) (公式 omega_hat omega - Im{(dSTFT/dt) / STFT}) if abs(STFT_coef) eps % 防止除以零 omega_inst f_stft(fi) - imag(dSTFT_dt(fi, ti) / STFT_coef) / (2*pi); % 上面除以2π是為了將角頻率(rad/s)轉(zhuǎn)換為普通頻率(Hz) else omega_inst f_stft(fi); end % 將瞬時(shí)頻率映射到最近的頻率倉(cāng)索引 k_hat round(omega_inst / df) 1; % 1 因?yàn)镸atlab索引從1開始 % 確保映射后的索引在有效范圍內(nèi) if k_hat 1 k_hat n_freq % 將當(dāng)前STFT系數(shù)累加到SST矩陣的對(duì)應(yīng)位置 % 注意這里是復(fù)數(shù)累加以保留重構(gòu)能力 SST(k_hat, ti) SST(k_hat, ti) STFT_coef; end end end end f_sst f_stft; t_sst t_stft; end在主腳本中調(diào)用這個(gè)函數(shù)% 設(shè)置SST參數(shù) thr max(abs(STFT(:))) * 0.01; % 閾值設(shè)為STFT最大幅值的1% % 計(jì)算SST [SST, f_sst, t_sst] my_sst(STFT, t_stft, f_stft, fs, hop, thr); % 繪制SST時(shí)頻譜 subplot(1,2,2); imagesc(t_sst, f_sst, 20*log10(abs(SST)eps)); % 加eps避免log10(0) axis xy; xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); title(‘同步壓縮變換時(shí)頻譜圖‘); colorbar; clim([-60, 0]); % 使用與STFT相同的顏色范圍3.4 結(jié)果對(duì)比分析與解讀運(yùn)行上述代碼后你會(huì)得到并排的兩幅時(shí)頻譜圖。對(duì)比它們你可以立即發(fā)現(xiàn)SST的魔力頻率聚焦性在STFT圖中兩個(gè)線性調(diào)頻信號(hào)是兩條較粗的、有一定寬度的“能量帶”。尤其是在時(shí)間中部約1秒處當(dāng)兩個(gè)信號(hào)的頻率非常接近時(shí)它們的能量帶會(huì)重疊、模糊在一起難以清晰分辨。而在SST圖中這兩條軌跡變成了極其銳利的細(xì)線即使它們靠得很近也能被清晰地區(qū)分開。這就是頻率分辨率的大幅提升。瞬態(tài)成分表征對(duì)于t1秒處的瞬時(shí)脈沖中心頻率250Hz在STFT圖中由于窗函數(shù)的限制它在時(shí)間軸上被“拉長(zhǎng)”了在頻率軸上也有一定的展寬看起來(lái)像一個(gè)“斑點(diǎn)”。在SST圖中這個(gè)脈沖在時(shí)間上依然被精確定位沒(méi)有因?yàn)閴嚎s而模糊時(shí)間信息同時(shí)在頻率上也變得更加集中更接近一個(gè)理想的時(shí)頻點(diǎn)。噪聲抑制觀察背景噪聲圖像中均勻分布的藍(lán)色背景。在SST圖中背景噪聲的強(qiáng)度似乎有所降低或變得更加“稀疏”。這是因?yàn)樵肼暤腟TFT系數(shù)相位是隨機(jī)的其計(jì)算出的瞬時(shí)頻率估計(jì)ω?也會(huì)非常隨機(jī)導(dǎo)致在再分配過(guò)程中能量被分散地映射到各個(gè)頻率倉(cāng)而不會(huì)像真實(shí)信號(hào)那樣集中到一條線上。因此在SST結(jié)果中信號(hào)的能量更加集中而噪聲的能量相對(duì)更加分散這在一定程度上提升了時(shí)頻譜的信噪比。注意事項(xiàng)2閾值選擇閾值thr的選擇至關(guān)重要。設(shè)得太高會(huì)丟失微弱信號(hào)設(shè)得太低會(huì)讓大量噪聲點(diǎn)參與壓縮不僅增加計(jì)算量還可能因噪聲點(diǎn)的隨機(jī)瞬時(shí)頻率估計(jì)而污染結(jié)果。通常建議設(shè)為STFT最大幅值的0.5%到5%之間需要根據(jù)具體信號(hào)的信噪比進(jìn)行微調(diào)。4. 關(guān)鍵參數(shù)影響與調(diào)優(yōu)指南SST的效果并非一勞永逸它嚴(yán)重依賴于初始STFT的參數(shù)設(shè)置。下面我們通過(guò)一個(gè)參數(shù)研究來(lái)理解這些影響。4.1 窗長(zhǎng)時(shí)頻分辨率的“起跑線”窗長(zhǎng)是STFT最核心的參數(shù)也是SST效果的基石。%% 4. 參數(shù)影響分析窗長(zhǎng) win_lens [64, 128, 256]; % 測(cè)試三種窗長(zhǎng) figure(‘Position‘, [100, 100, 1200, 900]); for i 1:length(win_lens) win_len win_lens(i); win hamming(win_len, ‘periodic‘); hop max(4, floor(win_len/16)); % 幀移隨窗長(zhǎng)適度增加 nfft 2^nextpow2(win_len*2); [STFT_temp, f_temp, t_temp] my_stft(x, win, hop, nfft, fs); thr_temp max(abs(STFT_temp(:))) * 0.01; [SST_temp, ~, ~] my_sst(STFT_temp, t_temp, f_temp, fs, hop, thr_temp); % 繪制STFT subplot(3, 2, (i-1)*21); imagesc(t_temp, f_temp, 20*log10(abs(STFT_temp))); axis xy; title([‘STFT - 窗長(zhǎng) ‘, num2str(win_len)]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); % 繪制SST subplot(3, 2, (i-1)*22); imagesc(t_temp, f_temp, 20*log10(abs(SST_temp)eps)); axis xy; title([‘SST - 窗長(zhǎng) ‘, num2str(win_len)]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); end結(jié)果分析短窗64點(diǎn)STFT的時(shí)間分辨率很高兩個(gè)調(diào)頻信號(hào)的軌跡在時(shí)間起止點(diǎn)清晰但頻率分辨率極差軌跡非常粗幾乎無(wú)法分辨中間接近的部分。SST試圖壓縮但“原料”太粗糙效果提升有限軌跡仍然較寬且可能出現(xiàn)斷點(diǎn)。中窗128點(diǎn)我們的初始選擇STFT的時(shí)頻權(quán)衡相對(duì)均衡。SST效果顯著軌跡銳利分離清晰。長(zhǎng)窗256點(diǎn)STFT的頻率分辨率很高兩條軌跡在頻率上本身已較清晰但時(shí)間分辨率下降瞬態(tài)脈沖被嚴(yán)重拉長(zhǎng)、模糊。SST能進(jìn)一步銳化頻率軌跡但無(wú)法修復(fù)損失的時(shí)間分辨率脈沖在SST中依然是被拉寬的。核心結(jié)論SST能顯著提升頻率分辨率但無(wú)法突破STFT初始時(shí)間分辨率的理論上限。它主要修復(fù)由窗函數(shù)引起的頻率方向上的能量擴(kuò)散。因此選擇窗長(zhǎng)的首要原則是確保STFT能捕捉到信號(hào)中最快的時(shí)間變化即所需的時(shí)間分辨率。在這個(gè)基礎(chǔ)上SST來(lái)優(yōu)化頻率分辨率。4.2 閾值信號(hào)與噪聲的“分水嶺”閾值決定了哪些STFT系數(shù)參與再分配。%% 5. 參數(shù)影響分析閾值 win_len 128; win hamming(win_len, ‘periodic‘); hop 4; nfft 256; [STFT_base, f_base, t_base] my_stft(x, win, hop, nfft, fs); thresholds [0.001, 0.01, 0.05]; % 相對(duì)于最大幅值的比例 figure(‘Position‘, [100, 100, 1200, 400]); for i 1:length(thresholds) thr max(abs(STFT_base(:))) * thresholds(i); [SST_temp, ~, ~] my_sst(STFT_base, t_base, f_base, fs, hop, thr); subplot(1, 3, i); imagesc(t_base, f_base, 20*log10(abs(SST_temp)eps)); axis xy; title([‘SST - 閾值 ‘, num2str(thresholds(i)*100), ‘%‘]); xlabel(‘時(shí)間 (s)‘); ylabel(‘頻率 (Hz)‘); clim([-60, 0]); colorbar; end結(jié)果分析低閾值0.1%幾乎所有點(diǎn)都參與壓縮包括大量噪聲點(diǎn)。結(jié)果圖中背景噪聲也呈現(xiàn)出一些虛假的、稀疏的“點(diǎn)狀”或“短線段”結(jié)構(gòu)這是因?yàn)樵肼暤碾S機(jī)相位導(dǎo)致了隨機(jī)的瞬時(shí)頻率估計(jì)。整體圖像可能看起來(lái)有點(diǎn)“臟”。適中閾值1%推薦起點(diǎn)大部分噪聲被過(guò)濾掉信號(hào)軌跡清晰銳利背景干凈。這是通常的起始選擇。高閾值5%只有能量最強(qiáng)的信號(hào)核心部分參與壓縮?？赡軐?dǎo)致微弱信號(hào)成分如我們信號(hào)中幅度為0.8的第二個(gè)成分的某些部分丟失軌跡出現(xiàn)不連續(xù)。同時(shí)瞬態(tài)脈沖的邊緣部分可能被舍棄。調(diào)優(yōu)建議從1%的閾值開始。如果發(fā)現(xiàn)微弱信號(hào)丟失適當(dāng)降低閾值如0.5%。如果背景噪聲干擾嚴(yán)重呈現(xiàn)虛假結(jié)構(gòu)則適當(dāng)提高閾值如2%?？梢杂^察SST結(jié)果中噪聲基底的特征來(lái)判斷。4.3 頻率軸插值與迭代SST基礎(chǔ)的SST將能量壓縮到離散的頻率網(wǎng)格上這可能導(dǎo)致“量化誤差”。更高級(jí)的實(shí)現(xiàn)可以采用以下技巧頻率軸插值在計(jì)算k_hat時(shí)不使用簡(jiǎn)單的round取整而是將能量按一定權(quán)重分配到相鄰的兩個(gè)頻率倉(cāng)上如線性插值這可以減輕因離散化造成的“柵欄效應(yīng)”使結(jié)果更平滑。這通常能帶來(lái)輕微的視覺改善。迭代SST將第一次SST的結(jié)果作為輸入再次進(jìn)行同步壓縮。理論上可以進(jìn)一步聚焦。但實(shí)踐中一次壓縮通常已能達(dá)到很好效果多次迭代可能引入偽影且計(jì)算成本翻倍。除非對(duì)時(shí)頻脊線提取有極高要求否則不建議常規(guī)使用。5. 常見問(wèn)題、排查技巧與進(jìn)階應(yīng)用在實(shí)際使用自編SST代碼時(shí)你可能會(huì)遇到以下典型問(wèn)題。5.1 時(shí)頻譜出現(xiàn)水平條紋或斷裂現(xiàn)象SST結(jié)果圖中本應(yīng)連續(xù)的信號(hào)軌跡出現(xiàn)明顯的水平斷裂帶或者在整個(gè)時(shí)間軸上出現(xiàn)均勻的水平條紋?？赡茉蚺c排查相位導(dǎo)數(shù)計(jì)算不準(zhǔn)確這是最常見的原因。確保計(jì)算dSTFT_dt時(shí)使用的是中心差分法并且時(shí)間步長(zhǎng)dt計(jì)算正確dt hop / fs。避免使用前向或后向差分它們?cè)谶吔缣幷`差大且整體精度低。邊界效應(yīng)我們的my_sst函數(shù)在計(jì)算時(shí)間導(dǎo)數(shù)時(shí)通過(guò)補(bǔ)零來(lái)近似中心差分但信號(hào)兩端的導(dǎo)數(shù)計(jì)算本身就不準(zhǔn)確。這會(huì)導(dǎo)致時(shí)間第一幀和最后一幀的瞬時(shí)頻率估計(jì)錯(cuò)誤從而產(chǎn)生邊界處的畸變。一種改進(jìn)方法是使用更復(fù)雜的邊界處理或者簡(jiǎn)單地在分析時(shí)舍棄頭尾幾幀。閾值過(guò)高過(guò)高的閾值會(huì)剔除掉構(gòu)成連續(xù)軌跡所必需的、幅度稍低的點(diǎn)導(dǎo)致軌跡斷裂。嘗試降低閾值。5.2 重構(gòu)信號(hào)誤差大現(xiàn)象使用SST的系數(shù)進(jìn)行信號(hào)重構(gòu)逆變換時(shí)重構(gòu)信號(hào)與原始信號(hào)差異顯著。可能原因與排查能量歸一化問(wèn)題SST是一個(gè)線性重分配過(guò)程但簡(jiǎn)單的“投票式”累加會(huì)改變系數(shù)的總能量。嚴(yán)格的可逆SST需要滿足保范數(shù)條件即在再分配過(guò)程中每個(gè)源點(diǎn)貢獻(xiàn)的能量權(quán)重需要精心設(shè)計(jì)使得整個(gè)變換是等距的。我們的基礎(chǔ)實(shí)現(xiàn)未做此處理因此逆變換不完美。若需精確重構(gòu)需查閱文獻(xiàn)實(shí)現(xiàn)“二階”或“可逆”SST。僅使用了SST的模如果只壓縮了能量abs(STFT)^2而丟棄了相位信息則絕對(duì)無(wú)法重構(gòu)。我們的代碼壓縮的是復(fù)數(shù)STFT保留了重構(gòu)的可能性。數(shù)值誤差累積相位導(dǎo)數(shù)的計(jì)算涉及除法對(duì)數(shù)值誤差敏感。確保使用雙精度計(jì)算并對(duì)極小分母進(jìn)行保護(hù)代碼中的eps檢查。5.3 對(duì)多分量信號(hào)中交叉軌跡的處理現(xiàn)象當(dāng)兩個(gè)信號(hào)的時(shí)頻軌跡在某個(gè)時(shí)間點(diǎn)交叉時(shí)SST結(jié)果在交叉點(diǎn)附近可能出現(xiàn)模糊或畸變。原因與對(duì)策這是SST以及大多數(shù)時(shí)頻后處理方法的固有挑戰(zhàn)。在交叉點(diǎn)信號(hào)的局部相位特性變得復(fù)雜瞬時(shí)頻率估計(jì)可能失效。對(duì)于交叉軌跡嘗試更短的窗短窗雖然初始頻率分辨率差但能更好地分離時(shí)間上快速變化的成分可能使交叉點(diǎn)的影響區(qū)域變小。使用方向性SST有研究提出在交叉區(qū)域根據(jù)信號(hào)分量方向進(jìn)行選擇性壓縮的算法但這非常復(fù)雜。接受局限對(duì)于高度非平穩(wěn)、分量交叉的信號(hào)需要認(rèn)識(shí)到時(shí)頻分析工具的局限性結(jié)合其他方法如經(jīng)驗(yàn)?zāi)B(tài)分解EMD先進(jìn)行信號(hào)分離再對(duì)單分量做SST。5.4 在強(qiáng)噪聲環(huán)境下的表現(xiàn)現(xiàn)象信號(hào)信噪比很低時(shí)SST可能無(wú)法清晰提取出信號(hào)軌跡甚至可能因噪聲而產(chǎn)生虛假結(jié)構(gòu)。優(yōu)化策略前置去噪在SST之前先對(duì)信號(hào)進(jìn)行濾波或小波去噪預(yù)處理。閾值調(diào)優(yōu)提高閾值只壓縮能量顯著高于噪聲基底的點(diǎn)。結(jié)合魯棒性估計(jì)使用更魯棒的瞬時(shí)頻率估計(jì)方法例如基于時(shí)頻分布如Wigner-Ville分布的重分配方法但計(jì)算量更大。多次平均如果條件允許對(duì)同一現(xiàn)象進(jìn)行多次測(cè)量對(duì)SST幅值譜進(jìn)行平均可以抑制隨機(jī)噪聲。5.5 計(jì)算效率優(yōu)化我們的雙循環(huán)實(shí)現(xiàn)直觀但較慢。對(duì)于長(zhǎng)信號(hào)或?qū)崟r(shí)處理可以考慮以下優(yōu)化向量化利用Matlab的矩陣運(yùn)算避免雙重循環(huán)?？梢酝瑫r(shí)計(jì)算所有點(diǎn)的瞬時(shí)頻率估計(jì)需處理除以零問(wèn)題并使用accumarray函數(shù)進(jìn)行高效的“投票”累加。這能帶來(lái)數(shù)量級(jí)的速度提升。使用C/MEX編碼將核心循環(huán)用C語(yǔ)言編寫并通過(guò)MEX接口調(diào)用適用于對(duì)性能要求極高的場(chǎng)合。利用GPUSST的并行性很好可以使用Matlab的Parallel Computing Toolbox或直接使用CUDA進(jìn)行GPU加速。最后分享一個(gè)我個(gè)人的深刻體會(huì)同步壓縮變換是一個(gè)極其強(qiáng)大的工具但它不是“銀彈”。它完美解決了單分量調(diào)頻信號(hào)在時(shí)頻譜上頻率擴(kuò)散的問(wèn)題。理解它的原理基于相位導(dǎo)數(shù)的再分配比單純調(diào)用一個(gè)函數(shù)更重要。這能幫助你在面對(duì)復(fù)雜信號(hào)時(shí)正確解讀SST的結(jié)果判斷哪些是真實(shí)的信號(hào)特征哪些可能是方法局限帶來(lái)的偽影。從STFT到SST就像是從一幅模糊的素描到一張清晰的高清照片而掌握拍攝參數(shù)設(shè)置和后期處理SST算法的技巧才能讓你成為真正的“信號(hào)攝影師”。