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