簡(jiǎn)介《數(shù)字信號(hào)處理理論、算法與實(shí)現(xiàn)》胡廣書(shū)著的配套Matlab代碼與參考文獻(xiàn)合集面向數(shù)字信號(hào)處理課程學(xué)習(xí)者、考研復(fù)習(xí)者及需要?jiǎng)邮址抡骝?yàn)證的工程技術(shù)人員。壓縮包含127個(gè)文件以114個(gè).m腳本為主覆蓋離散信號(hào)與系統(tǒng)、z變換、傅里葉變換、統(tǒng)計(jì)信號(hào)處理等經(jīng)典章節(jié)的算法實(shí)現(xiàn)另含7個(gè)wav音頻樣本和1個(gè)mat數(shù)據(jù)文件便于直接運(yùn)行和觀察處理效果還附帶3個(gè)pdf參考文獻(xiàn)及1個(gè)doc光盤(pán)使用說(shuō)明幫助讀者快速上手并對(duì)照教材學(xué)習(xí)。資源包整體僅4.5MB小巧但內(nèi)容集中已有1663人學(xué)習(xí)下載。通過(guò)運(yùn)行這些代碼可直觀理解濾波器設(shè)計(jì)、頻譜分析、LOT變換等核心方法的編程實(shí)現(xiàn)思路節(jié)省自行編寫(xiě)和調(diào)試基礎(chǔ)算法的時(shí)間是將教材理論轉(zhuǎn)化為可執(zhí)行程序的實(shí)用參考資料。1. 數(shù)字信號(hào)處理為什么照著教材寫(xiě)代碼頻譜圖還是比書(shū)上難看“數(shù)字信號(hào)處理”是那種看著推導(dǎo)容易、動(dòng)手復(fù)現(xiàn)就露怯的領(lǐng)域。有一回我按教材里的周期圖法估一個(gè) 10 Hz 正弦的功率譜書(shū)上的圖干凈得像一根針我跑出來(lái)的譜線卻毛刺橫生連續(xù)調(diào)了幾天參數(shù)都?jí)翰黄健Ｗ詈蟀l(fā)現(xiàn)錯(cuò)不在算法本身而在我沒(méi)搞懂 DFT 點(diǎn)數(shù)、采樣頻率和隱式周期延拓三者之間的耦合關(guān)系。教材配套代碼的價(jià)值就在這里它把每一段數(shù)學(xué)推導(dǎo)變成可執(zhí)行程序逼著你把公式里的符號(hào)和程序里的變量一一對(duì)上號(hào)。搞懂這套代碼等于給 FFT、濾波器設(shè)計(jì)、功率譜估計(jì)這些基礎(chǔ)機(jī)制都裝上一個(gè)可調(diào)試的心理模型。之后再看任何頻域問(wèn)題你都能直接定位到是“哪兒算錯(cuò)了”而不是“玄學(xué)問(wèn)題”。2. 離散傅里葉變換點(diǎn)數(shù)決定分辨率補(bǔ)零不是后悔藥2.1 從連續(xù)到離散周期延拓如何決定頻譜外觀先看 DFT 的定義X[k] sum_{n0}^{N-1} x[n] * exp(-j*2*pi*k*n/N)。這個(gè)式子里藏著三個(gè)隱含假設(shè)序列有限長(zhǎng)、頻域離散、時(shí)域周期延拓。特別是最后一條——N 點(diǎn) DFT 默認(rèn)x[n]是以 N 為周期的信號(hào)被取了一個(gè)主值區(qū)間哪怕你喂進(jìn)去的只是一段錄音、一段振動(dòng)波形算法也強(qiáng)制認(rèn)為它首尾相接、循環(huán)播放。大多數(shù)真實(shí)信號(hào)截?cái)嗟?N 點(diǎn)之后首尾的值對(duì)不上延拓出來(lái)的周期信號(hào)在邊界處產(chǎn)生跳變頻域上就表現(xiàn)為能量從真實(shí)譜線位置向整個(gè)頻帶蔓延這就是頻譜泄漏。教材里的示例信號(hào)幾乎都是精心挑選的頻率恰好落在整數(shù)譜線 bin 上、采樣點(diǎn)數(shù)剛好覆蓋整數(shù)個(gè)周期畫(huà)出來(lái)當(dāng)然漂亮。你自己隨便采一段數(shù)據(jù)頻率落在兩個(gè) bin 之間譜峰就會(huì)變矮、變寬旁邊還多出一串旁瓣。所以書(shū)上的圖“好看”不是算法的問(wèn)題是信號(hào)選取的問(wèn)題。這個(gè)認(rèn)知直接決定了后來(lái)的一切判斷。你看到頻譜圖上出現(xiàn)不該有的分量第一反應(yīng)不應(yīng)該是“換一個(gè)更高級(jí)的窗”而是先算一下信號(hào)周期和采樣點(diǎn)數(shù)是否匹配。周期延拓這個(gè)前提一旦沒(méi)建立起來(lái)后面所有頻域操作都會(huì)在錯(cuò)誤的坐標(biāo)系里打轉(zhuǎn)。2.2 分辨率公式 dffs/N補(bǔ)零插值改變不了物理極限要分辨兩個(gè)頻率挨得很近的正弦關(guān)鍵指標(biāo)是頻率分辨率df fs / N。fs 固定時(shí)只有增加采樣點(diǎn)數(shù) N——也就是延長(zhǎng)采集時(shí)長(zhǎng)——才能把譜線間距壓小。舉個(gè)例子fs 1000N 512 時(shí)df ≈ 1.95 Hz兩個(gè)相隔 2 Hz 的譜峰理論上勉強(qiáng)分開(kāi)但矩形窗主瓣寬度約兩倍 df實(shí)際上糊在一起。如果你把信號(hào)補(bǔ)零到 4096 再 FFT譜線確實(shí)變得更平滑但兩個(gè)峰還是疊著因?yàn)檠a(bǔ)零沒(méi)有增加任何新信息。fs 1000; N 512; % 實(shí)際采樣點(diǎn)數(shù) t (0:N-1) / fs; x sin(2*pi*100*t) sin(2*pi*102*t); % 間隔 2 Hz 的兩個(gè)正弦 X1 fft(x); % 512 點(diǎn) FFT X2 fft(x, 4096); % 補(bǔ)零到 4096 點(diǎn)再做 FFT f1 (0:N/2) * fs / N; f2 (0:2048) * fs / 4096; subplot(2,1,1); plot(f1, abs(X1(1:N/21))); title(512 點(diǎn)直接 FFT); subplot(2,1,2); plot(f2, abs(X2(1:2049))); title(補(bǔ)零到 4096 點(diǎn));補(bǔ)零后的頻譜曲線看起來(lái)更“圓潤(rùn)”每個(gè)峰之間多了很多采樣點(diǎn)但這兩個(gè)峰能不能被區(qū)分取決于原始 512 點(diǎn)數(shù)據(jù)的時(shí)長(zhǎng)N/fs 0.512 s。1/0.512 ≈ 1.95 Hz正好對(duì)應(yīng) df。你要真的分開(kāi) 2 Hz 的間隔就得把 N 加到 1000 以上讓采集時(shí)長(zhǎng)超過(guò) 1 秒而不是在補(bǔ)零上做文章。換個(gè)說(shuō)法補(bǔ)零改的是顯示分辨率物理分辨率在你決定采多少個(gè)點(diǎn)的瞬間就已經(jīng)固定了。注意補(bǔ)零改變的是顯示分辨率不是物理分辨率。物理分辨率在數(shù)據(jù)采集完成那一刻就固定了后續(xù)任何插值算法都補(bǔ)不回沒(méi)有采到的時(shí)間。2.3 一版能直接跑的頻譜分析流程理解了前兩節(jié)再看一套完整的頻譜分析流程每一步都有了落點(diǎn)。以下是我常用的單邊幅度譜模板fs 1000; % 采樣率 1 kHz N 1024; % 采樣點(diǎn)數(shù) t (0:N-1) / fs; x sin(2*pi*50*t) 0.5*sin(2*pi*123.4*t); x x - mean(x); % 去直流防止 0 Hz 大分量壓低譜峰 w hann(N, periodic).; % 周期 Hann 窗注意轉(zhuǎn)置成行向量 X fft(x .* w); % 加窗后做 N 點(diǎn) FFT f (0:N/2) * fs / N; % 單邊譜頻率軸 A 2 * abs(X(1:N/21)) / sum(w); % 幅度歸一化因子用窗和不是 N plot(f, 20*log10(A)); xlabel(頻率 (Hz)); ylabel(幅度 (dB)); grid on;這里每個(gè)細(xì)節(jié)都有講究。去直流是因?yàn)?0 Hz 處的大分量會(huì)通過(guò)窗函數(shù)的旁瓣污染附近頻點(diǎn)選周期 Hann 窗而不是對(duì)稱(chēng) Hann 窗是因?yàn)?DFT 隱含周期延拓周期窗能保證窗本身首尾連續(xù)降低泄漏幅度歸一化除以sum(w)而不是 N是因?yàn)榧哟昂笮盘?hào)總能量變了用 N 去歸一化會(huì)把幅度低估將近一半單邊譜只取前N/21個(gè)點(diǎn)并把除直流和奈奎斯特點(diǎn)外的幅度乘 2把雙邊譜能量折回單邊。教材配套代碼里很多程序直接畫(huà)abs(fft(x))那是為了展示算法結(jié)構(gòu)不是工程標(biāo)準(zhǔn)對(duì)照著看能明白作者簡(jiǎn)化了什么又為什么敢這么簡(jiǎn)化。3. FIR 濾波器設(shè)計(jì)窗函數(shù)法為什么總在邊界處翻車(chē)3.1 理想濾波器的沖擊響應(yīng)是無(wú)限長(zhǎng)的 sinc窗函數(shù)法設(shè)計(jì) FIR 濾波器的起點(diǎn)是一個(gè)理想化的頻域門(mén)通帶內(nèi)幅度為 1阻帶內(nèi)幅度為 0。對(duì)這個(gè)理想頻率響應(yīng)做逆傅里葉變換得到的時(shí)域沖擊響應(yīng)是一個(gè)無(wú)限長(zhǎng)的 sinc 函數(shù)。實(shí)際能實(shí)現(xiàn)的 FIR 濾波器必然有限長(zhǎng)于是只能把這個(gè)無(wú)限長(zhǎng)序列截?cái)?。fs 8000; % 采樣率 8 kHz fc 1000; % 截止頻率 1 kHz M 31; % 濾波器長(zhǎng)度奇數(shù) n 0:M-1; hd 2*fc/fs * sinc(2*fc*(n-(M-1)/2)/fs); % 理想低通沖擊響應(yīng)注意sinc在 MATLAB 里是歸一化形式sin(pi*x)/(pi*x)所以這里系數(shù)直接寫(xiě)2*fc/fs單位是“周期/樣本”。(M-1)/2是為了把沖擊響應(yīng)移到濾波器中間保證因果性。為什么 M 取奇數(shù)因?yàn)榕紨?shù)長(zhǎng)度時(shí)群延遲(M-1)/2是半樣本相位響應(yīng)會(huì)引入非線性教材里絕大多數(shù)窗函數(shù)法程序都刻意避免這種情況。這段代碼展開(kāi)的 hd 在中心點(diǎn)附近值最大向兩邊衰減震蕩——這就是頻域“門(mén)”在時(shí)域里的真實(shí)形狀。3.2 窗函數(shù)的取舍旁瓣與過(guò)渡帶之間的交易截?cái)啾旧砭褪浅艘跃匦未熬匦未暗牡谝慌园曛挥屑s -13 dB阻帶衰減遠(yuǎn)達(dá)不到工程要求。所以標(biāo)準(zhǔn)做法是截?cái)嘀笤俪艘粋€(gè)低旁瓣窗讓時(shí)域序列兩端平滑衰減到零代價(jià)是主瓣變寬過(guò)渡帶變寬。窗類(lèi)型主瓣寬度第一旁瓣衰減矩形窗4π/M-13 dBHann 窗8π/M-31 dBHamming 窗8π/M-43 dB這里有個(gè)初學(xué)者很難繞出來(lái)的彎加大濾波器長(zhǎng)度 M 會(huì)讓過(guò)渡帶變窄但通帶波紋不會(huì)跟著變小——通帶波紋的極限由窗類(lèi)型決定M 只影響波紋的“抖動(dòng)密度”。教科書(shū)里那段“增加 M 不能改善阻帶衰減”的結(jié)論代碼里驗(yàn)證方式很簡(jiǎn)單把 M 從 31 改成 127用 freqz 看阻帶衰減仍然是同樣的 -53 dB 左右只是過(guò)渡帶明顯變陡。w hamming(M).; % 轉(zhuǎn)置成行向量保證和 hd 維度一致 h hd .* w; % 加窗后的實(shí)際濾波器系數(shù)3.3 用 freqz 把設(shè)計(jì)指標(biāo)翻譯回頻域設(shè)計(jì)完系數(shù)第一件事不是直接濾波而是用 freqz 看頻率響應(yīng)[h_resp, f] freqz(h, 1, 1024, fs); % 1024 點(diǎn)橫軸單位為 Hz plot(f, 20*log10(abs(h_resp)));這一行把頻域坐標(biāo)全部拉直橫軸是真實(shí) Hz縱軸是 dB?？?-3 dB 點(diǎn)是否落在 1000 Hz、阻帶是否達(dá)到設(shè)計(jì)指標(biāo)、通帶波紋是否在可接受范圍內(nèi)。教材配套代碼里這一節(jié)的通常做法是打印20*log10(abs(h))直接觀察但你自己畫(huà)一遍 freqz 會(huì)看到教科書(shū)不會(huì)強(qiáng)調(diào)的細(xì)節(jié)過(guò)渡帶中心處的衰減速度、通帶邊緣的吉布斯余暉以及不同窗函數(shù)作用在同一個(gè) hd 上的差異。工程上真正跑濾波之前這套檢查是必須的不然信號(hào)進(jìn)去什么樣、出來(lái)什么樣全靠黑匣子猜。4. 功率譜估計(jì)周期圖方差大的根源與 Welch 平均4.1 周期圖法的方差不隨 N 增大而收斂周期圖法是最直觀的譜估計(jì)對(duì)一段數(shù)據(jù)做 FFT取幅度平方再除以 N就得到功率譜。它的問(wèn)題非常反直覺(jué)——數(shù)據(jù)越長(zhǎng)頻率分辨率越高譜線的方差卻不下降。一段白噪聲不管采 1024 點(diǎn)還是 8192 點(diǎn)周期圖在任意一個(gè)頻點(diǎn)上的相對(duì)標(biāo)準(zhǔn)差都在 100% 附近晃蕩因?yàn)槊總€(gè)頻點(diǎn)本質(zhì)上只有兩個(gè)自由度實(shí)部、虛部增加點(diǎn)數(shù)只是讓頻點(diǎn)變得更密每個(gè)頻點(diǎn)擁有的數(shù)據(jù)量并不增加。N 1024; Pall zeros(50, N/21); for trial 1:50 x randn(1, N); % 高斯白噪聲 X fft(x); P abs(X(1:N/21)).^2 / N; % 周期圖 Pall(trial,:) P; end rel_std std(Pall, 0, 1) ./ mean(Pall, 1); fprintf(周期圖頻點(diǎn)相對(duì)標(biāo)準(zhǔn)差均值: %.2f\n, mean(rel_std));跑完這段代碼你算出來(lái)的相對(duì)標(biāo)準(zhǔn)差均值通常接近 1.0也就是說(shuō)譜估計(jì)的隨機(jī)起伏和信號(hào)本身差不多大。這就是為什么直接把a(bǔ)bs(fft(x)).^2/N當(dāng)功率譜用畫(huà)出來(lái)的圖難看得沒(méi)法看。明白這一點(diǎn)才能理解 Welch 平均法到底在解決什么問(wèn)題——它犧牲了分辨率換取了方差的收斂。4.2 手寫(xiě) Welch 平均法分幀、加窗、平均Welch 法的核心是把 N 點(diǎn)數(shù)據(jù)切成若干段每段加窗后做周期圖再把所有段的周期圖平均。段數(shù)越多方差越小每段越短頻率分辨率越差。這兩者是直接對(duì)立的。一段常規(guī)配置是 50% 重疊、Hann 窗相鄰幀共享一半數(shù)據(jù)在方差改善和計(jì)算量之間取平衡。N 1024; fs 1000; M 256; % 每幀長(zhǎng)度 overlap 128; % 50% 重疊 nseg fix((N-M)/(M-overlap)) 1; % 7 幀 w hann(M, periodic); Pxx zeros(M/21, 1); for k 1:nseg idx (k-1)*(M-overlap) (1:M); xk x(idx) .* w; % 逐幀加窗 Xk fft(xk); P abs(Xk(1:M/21)).^2 / (fs * sum(w.^2)); % PSD 量綱 Pxx Pxx P; end Pxx Pxx / nseg; % 平均 f (0:M/2) * fs / M;這段代碼里最容易忽略的是最后那個(gè)歸一化因子。fs * sum(w.^2)是等效噪聲帶寬用來(lái)把 FFT 輸出轉(zhuǎn)換為“每赫茲功率”——功率譜密度的量綱。如果只是要比較相對(duì)形狀用sum(w)^2做歸一化也行但如果你想把譜密度和信號(hào)實(shí)際功率對(duì)應(yīng)上就必須用這個(gè)因子。另外注意分辨率已經(jīng)退化成fs/M 1000/256 ≈ 3.9 Hz相比整段 1024 點(diǎn)的 0.98 Hz每個(gè)峰都會(huì)被抹得更寬。4.3 用重復(fù)實(shí)驗(yàn)驗(yàn)證噪聲抑制效果我自己驗(yàn)證 Welch 效果的方式是對(duì)同一個(gè)隨機(jī)過(guò)程做 30 次獨(dú)立實(shí)驗(yàn)每次分別用周期圖法和 Welch 法估計(jì)同一頻點(diǎn)的功率然后比較兩者的標(biāo)準(zhǔn)差。std_periodogram zeros(1, 30); std_welch zeros(1, 30); for trial 1:30 x randn(1, N); X fft(x); P_period abs(X(1:N/21)).^2 / N; std_periodogram(trial) std(P_period(10:100)); % Pxx 來(lái)自上一節(jié)的 Welch 循環(huán) std_welch(trial) std(Pxx(10:50)); end fprintf(周期圖平均標(biāo)準(zhǔn)差: %.2f\n, mean(std_periodogram)); fprintf(Welch 平均標(biāo)準(zhǔn)差: %.2f\n, mean(std_welch));Welch 的平均標(biāo)準(zhǔn)差大致是周期圖的1/sqrt(nseg)量級(jí)和理論預(yù)期吻合。這套驗(yàn)證思路可以直接復(fù)用到任何譜估計(jì)算法上先造一個(gè)已知統(tǒng)計(jì)特性的信號(hào)用程序跑出估計(jì)結(jié)果再對(duì)比理論值而不是拿真實(shí)測(cè)量數(shù)據(jù)“感覺(jué)差不多就行”。后面在資源里的參考文獻(xiàn)部分關(guān)于譜估計(jì)的那幾篇老文章專(zhuān)門(mén)討論過(guò)這種統(tǒng)計(jì)驗(yàn)證的必要性建議對(duì)照著看。5. 常見(jiàn)問(wèn)題避坑指南五條最容易踩的坑5.1 把歸一化頻率當(dāng)成 Hz 用濾波器完全對(duì)不上現(xiàn)象設(shè)計(jì)一個(gè)截止 1 kHz 的低通寫(xiě)成fir1(30, fc/fs)實(shí)測(cè)截止點(diǎn)在 500 Hz 附近整個(gè)濾波器只有預(yù)期的一半帶寬。原因fir1的歸一化截止頻率是相對(duì)奈奎斯特頻率fs/2的所以 Wn 應(yīng)該寫(xiě)成fc/(fs/2)等于手寫(xiě) sinc 公式里的2*fc/fs。誤寫(xiě)成fc/fs相當(dāng)于把 Wn 縮小了一半。解決統(tǒng)一約定。寫(xiě)fir1參數(shù)時(shí)用fc_norm fc/(fs/2)手寫(xiě)理想低通沖擊響應(yīng)時(shí)用2*fc/fs。畫(huà)頻率響應(yīng)圖時(shí)永遠(yuǎn)帶上fs參數(shù)讓橫軸顯示 Hz別讓程序猜單位。5.2 指望補(bǔ)零提高分辨率兩個(gè)譜峰還是糊在一起現(xiàn)象數(shù)據(jù)只采了 512 點(diǎn)兩個(gè)頻率差 2 Hz 的信號(hào)分不開(kāi)補(bǔ)零到 8192 點(diǎn)后頻譜曲線變得很細(xì)但兩個(gè)峰依然重疊成一個(gè)寬峰。原因補(bǔ)零是對(duì)已有頻譜做 sinc 插值讓曲線看起來(lái)更平滑它不產(chǎn)生任何新信息。能否分開(kāi)兩個(gè)峰取決于原始數(shù)據(jù)的時(shí)長(zhǎng)N/fs是否大于1/Δf。解決先算最小需求N_min fs / Δf然后回頭檢查采集時(shí)長(zhǎng)。數(shù)據(jù)長(zhǎng)度不夠就去補(bǔ)采軟件手段解決不了信息缺失。類(lèi)似地不要用“補(bǔ)零后峰值變高了”來(lái)評(píng)估幅度那是插值函數(shù)的貢獻(xiàn)。5.3 浮點(diǎn)下穩(wěn)穩(wěn)的 IIR 濾波器一換定點(diǎn)就爆炸現(xiàn)象雙線性變換設(shè)計(jì)的巴特沃斯低通在浮點(diǎn)仿真里表現(xiàn)完美移植到定點(diǎn) DSP 上輸出逐漸發(fā)散最終飽和成一條直線。原因IIR 濾波器極點(diǎn)離單位圓太近系數(shù)用 16 bit 量化后極點(diǎn)的微小位移可能直接越過(guò)單位圓濾波器從穩(wěn)定變成不穩(wěn)定。階數(shù)越高、極點(diǎn)越密集對(duì)量化誤差越敏感。這是我在某個(gè)傳感器項(xiàng)目里踩得最深的一次浮點(diǎn)圖和定點(diǎn)圖完全是兩個(gè)世界。解決不要直接量化傳遞函數(shù)的多項(xiàng)式系數(shù)先轉(zhuǎn)換成二階節(jié)SOS級(jí)聯(lián)每節(jié)單獨(dú)量化量化后畫(huà)出零極點(diǎn)圖確認(rèn)所有極點(diǎn)仍在單位圓內(nèi)最后用定點(diǎn)化系數(shù)重新跑一遍仿真對(duì)比阻帶衰減是否還在指標(biāo)內(nèi)。這個(gè)檢查流程要固化成腳本每次設(shè)計(jì)完 IIR 都強(qiáng)制走一遍。5.4 filter 和 conv 的長(zhǎng)度與瞬態(tài)不一致現(xiàn)象同一個(gè)濾波器用filter(h,1,x)和conv(x,h)得到的結(jié)果長(zhǎng)度不一樣逐點(diǎn)對(duì)比前幾個(gè)點(diǎn)也完全對(duì)不上于是懷疑程序有問(wèn)題。原因conv輸出完整線性卷積長(zhǎng)度是len(x)len(h)-1filter是差分方程仿真輸出長(zhǎng)度和輸入相同。對(duì)單次有限長(zhǎng)輸入filter的前l(fā)en(x)個(gè)點(diǎn)其實(shí)和conv對(duì)應(yīng)位置相等但兩者尾部不同——conv多出的部分是脈沖響應(yīng)越過(guò)輸入末尾后的“拖尾”。解決比較響應(yīng)時(shí)用conv且注意截取對(duì)齊做流式處理時(shí)一定用filter的第四個(gè)參數(shù)傳遞初始狀態(tài)[y, zf] filter(h, 1, x, zi);把上一幀的zf傳給下一幀的zi否則每幀開(kāi)頭的瞬態(tài)會(huì)污染輸出常見(jiàn)的“分幀濾波咔嗒聲”就是這么來(lái)的。5.5 功率譜估計(jì)漏掉 fs 的量綱縱軸對(duì)不上現(xiàn)象用pwelch算出來(lái)的譜密度和手寫(xiě)abs(fft(x)).^2/N結(jié)果不在一個(gè)數(shù)量級(jí)差出好幾倍甚至幾十倍。原因功率譜密度PSD的單位是“每 Hz 的功率”需要除以等效噪聲帶寬fs * sum(w.^2)而abs(fft(x)).^2/N是周期圖法下的總功率譜兩者量綱本就不同。教材里周期圖程序常直接畫(huà)后者因?yàn)樗魂P(guān)心譜峰形狀和相對(duì)位置不關(guān)心絕對(duì)功率。解決在代碼注釋里寫(xiě)清“這里算的是 PSD 還是功率譜”要和其他工具對(duì)齊就用標(biāo)準(zhǔn)函數(shù)手寫(xiě)版本用我第 4 章給過(guò)的歸一化方式。別在同一個(gè)圖里混著對(duì)比兩種量綱的數(shù)據(jù)縱軸對(duì)不上是第一眼就能發(fā)現(xiàn)的信號(hào)。6. 浮點(diǎn)與定點(diǎn)算法驗(yàn)證之后的設(shè)計(jì)取舍6.1 為什么教材代碼清一色浮點(diǎn)教材里的 MATLAB 代碼全部跑在雙精度浮點(diǎn)上這不是偷懶而是刻意為之。算法驗(yàn)證階段的目標(biāo)是確認(rèn)“算法本身能不能 work”此時(shí)把量化效應(yīng)混進(jìn)來(lái)會(huì)分不清誤差來(lái)自算法還來(lái)自數(shù)值實(shí)現(xiàn)。FIR 濾波器系數(shù)量化相對(duì)溫和IIR 則完全不同——極點(diǎn)的位置對(duì)系數(shù)誤差極其敏感尤其高階設(shè)計(jì)里極點(diǎn)密集在一起時(shí)微小擾動(dòng)就可以讓一對(duì)共軛極點(diǎn)移出單位圓。浮點(diǎn)仿真里永遠(yuǎn)看不到這個(gè)問(wèn)題換成定點(diǎn)實(shí)現(xiàn)立刻爆炸。6.2 一個(gè)簡(jiǎn)單的定點(diǎn)化預(yù)演我現(xiàn)在的習(xí)慣是任何濾波器設(shè)計(jì)在交付前都強(qiáng)制做一遍定點(diǎn)化預(yù)演哪怕目標(biāo)平臺(tái)是浮點(diǎn)處理器也照做h fir1(60, 0.2); % 浮點(diǎn) FIR 原型 h_fix round(h / max(abs(h)) * 32767) / 32767; % 16 bit 量化 [H1, f] freqz(h, 1, 1024, 1); [H2, ~] freqz(h_fix, 1, 1024, 1); plot(f, 20*log10(abs(H1)), b, f, 20*log10(abs(H2)), r); legend(浮點(diǎn), 16-bit 定點(diǎn));FIR 量化后通常只是阻帶衰減略微變差還能接受IIR 要更嚴(yán)苛地檢查零極點(diǎn)位置用 zplane 看量化后極點(diǎn)是否還在單位圓內(nèi)。有一次我在某溫控系統(tǒng)里用了一個(gè)四階巴特沃斯帶通浮點(diǎn)仿真輸出平滑穩(wěn)定定點(diǎn)化后系統(tǒng)開(kāi)始自激振蕩查了一整天才定位到是一對(duì)共軛極點(diǎn)被量化誤差推出了單位圓。從那以后我每次從這套教材配套代碼里拿到濾波器設(shè)計(jì)都會(huì)先跑一遍定點(diǎn)化仿真再談部署這個(gè)習(xí)慣救過(guò)我很多次。類(lèi)似的驗(yàn)證做法也可以反向用如果你拿到一個(gè)嵌入式上的濾波器實(shí)現(xiàn)先把它轉(zhuǎn)成雙精度浮點(diǎn)在 MATLAB 里重建一遍看看和原設(shè)計(jì)差在哪能定位出不少隱藏問(wèn)題。需要把代碼在本地逐行跑一遍的讀者這套教材的配套 MATLAB 代碼和參考文獻(xiàn)都在對(duì)應(yīng)資源包里建議從第 2 章的 FFT 程序開(kāi)始對(duì)著本文提到的幾個(gè)驗(yàn)證思路逐章核對(duì)希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取