簡(jiǎn)介這是一份用 C 語(yǔ)言實(shí)現(xiàn)的 Parks-McClellan 算法remez-exchange源碼面向數(shù)字信號(hào)處理學(xué)習(xí)者、嵌入式開發(fā)者以及需要設(shè)計(jì)線性相位 FIR 濾波器的工程師。算法采用切比雪夫逼近和交替誤差交換策略通過迭代優(yōu)化濾波器系數(shù)在指定頻段內(nèi)獲得等波紋equiripple頻率響應(yīng)適合教學(xué)實(shí)驗(yàn)、濾波器原型驗(yàn)證、通信系統(tǒng)仿真以及工程移植等場(chǎng)景。壓縮包共 3 個(gè)文件包含 C 語(yǔ)言源文件、Markdown 說(shuō)明文檔和許可證文本整體約 10KB結(jié)構(gòu)精簡(jiǎn)便于快速閱讀和直接編譯。已有 738 人下載學(xué)習(xí)說(shuō)明該源碼在 FIR 濾波器設(shè)計(jì)場(chǎng)景中具有一定參考價(jià)值。使用者可直接運(yùn)行源碼生成濾波器系數(shù)也可結(jié)合說(shuō)明文檔理解迭代優(yōu)化的完整步驟并根據(jù)濾波器階數(shù)、通帶/阻帶邊界和權(quán)重等需求修改參數(shù)從而深入掌握 Parks-McClellan 算法的原理及 C 語(yǔ)言實(shí)現(xiàn)細(xì)節(jié)。1. Parks-McClellan算法為什么“誤差全部相等”才是FIR設(shè)計(jì)的最優(yōu)解理想低通濾波器在物理上不可實(shí)現(xiàn)——頻域是矩形時(shí)域就是無(wú)限長(zhǎng)的sinc。工程上所有FIR設(shè)計(jì)方法本質(zhì)都在回答同一個(gè)問題用一個(gè)有限長(zhǎng)沖激響應(yīng)去逼近那個(gè)不可能得到的理想響應(yīng)逼近得“最好”。窗函數(shù)法靠截?cái)囝l率采樣法靠插值而Parks-McClellan算法給出的答案完全不同它把設(shè)計(jì)問題轉(zhuǎn)化成切比雪夫逼近用Remez交換迭代去解最后得到的濾波器誤差曲線是一條上下均勻起伏的“等波動(dòng)”帶通帶和阻帶里每個(gè)峰值誤差都恰好一樣大。這正是這個(gè)算法最反直覺的地方最優(yōu)解不是讓誤差“平均最小”而是讓“最大誤差”最小并且所有最大誤差點(diǎn)嚴(yán)格相等、符號(hào)交替。要理解為什么這個(gè)條件就是最優(yōu)以及一段用C寫出來(lái)的remez程序憑什么能迭代收斂需要從頭拆開交替定理和Remez交換的機(jī)制。2. 從窗函數(shù)到切比雪夫逼近要“平均誤差小”還是“最大誤差小”2.1 窗函數(shù)法輸在哪旁瓣不是你能控制的窗函數(shù)法的思路是把理想沖激響應(yīng) h_d[n] 直接乘以一個(gè)有限長(zhǎng)窗口 w[n]得到 h[n] h_d[n]·w[n]。時(shí)域相乘等于頻域卷積理想矩形頻譜和窗口頻譜卷積的結(jié)果決定了通帶邊緣的過渡帶寬度和阻帶衰減。問題在于一旦你選了窗口類型過渡帶寬度和阻帶衰減的配對(duì)關(guān)系就被鎖死了。矩形窗阻帶衰減只有 -21dB漢寧窗是 -44dB海明窗是 -53dB布萊克曼窗是 -74dB。你沒法說(shuō)“我要 -60dB 衰減、過渡帶 0.01π”——窗函數(shù)只給你固定配方的折中指標(biāo)是切不出來(lái)的。更麻煩的是通帶和阻帶共享同一個(gè)卷積核通帶波紋和阻帶波紋不可分別調(diào)節(jié)這對(duì)很多實(shí)際系統(tǒng)是不夠用的。換窗能改善旁瓣但改善的機(jī)制是把能量從旁瓣推遠(yuǎn)代價(jià)是主瓣變寬、過渡帶變緩。這是時(shí)域截?cái)噙@個(gè)動(dòng)作本身帶來(lái)的物理限制不是實(shí)現(xiàn)精度問題。所以窗函數(shù)法適合快速出原型不適合做有嚴(yán)格指標(biāo)的濾波器。2.2 頻率采樣法輸在哪只在采樣點(diǎn)上精確頻率采樣法把理想響應(yīng)在N個(gè)等間隔頻率點(diǎn)上采樣然后 IDFT 得到系數(shù)。它的優(yōu)點(diǎn)是設(shè)計(jì)過程零迭代、極快缺點(diǎn)是只在采樣點(diǎn)上嚴(yán)格等于目標(biāo)值兩個(gè)采樣點(diǎn)之間的響應(yīng)完全不受控。當(dāng)理想響應(yīng)在過渡帶邊緣有陡峭跳變時(shí)點(diǎn)間會(huì)產(chǎn)生很大的振蕩誤差峰值甚至比窗函數(shù)法還大?？梢赃@么理解頻率采樣法是一個(gè)在N個(gè)點(diǎn)上精確、在其他頻率放任自流的插值逼近。它沒有把“整個(gè)頻帶上的誤差”作為一個(gè)整體來(lái)優(yōu)化所以談不上最優(yōu)。2.3 極小極大準(zhǔn)則從“平均誤差”換到“最大誤差”既然前兩種方法都缺乏全局誤差控制正確的數(shù)學(xué)化方式是什么逼近論里有一個(gè)現(xiàn)成的框架設(shè)理想響應(yīng)為 D(ω)設(shè)計(jì)出的濾波器頻率響應(yīng)為 A(ω)誤差為 E(ω) W(ω)·(A(ω) ? D(ω))其中 W(ω) 是用戶給的權(quán)重函數(shù)——你對(duì)哪個(gè)頻帶更敏感就給哪個(gè)頻帶更大的權(quán)重。切比雪夫逼近也叫極小極大逼近的目標(biāo)是minimize max |E(ω)|, ω ∈ 通帶 ∪ 阻帶也就是說(shuō)最小化整個(gè)頻帶上的最大加權(quán)誤差。注意這里不含過渡帶——過渡帶沒有目標(biāo)值誤差在那里不參與優(yōu)化這正是它存在的意義過渡帶越窄代價(jià)是通帶和阻帶的波紋越大反之亦然。用戶用過渡帶寬度換取帶內(nèi)誤差這是其他方法給不了的旋鈕。而“最大誤差最小”這件事遠(yuǎn)比聽起來(lái)苛刻。它要求最優(yōu)解在所有局部極值點(diǎn)上的誤差幅度完全一樣——否則某個(gè)大誤差點(diǎn)可以被壓低直到它和其他極值點(diǎn)一樣高。最終誤差曲線必然是等波動(dòng)形狀這正是Parks-McClellan算法產(chǎn)物的視覺簽名。3. 交替定理最佳逼近的“指紋”以及Remez交換為什么能收斂3.1 交替定理究竟說(shuō)了什么對(duì)于抽頭數(shù)為N的線性相位FIRN為奇數(shù)Type I對(duì)稱其幅度響應(yīng) A(ω) 可以寫成 cos(ω) 的 M 次三角多項(xiàng)式其中 M (N?1)/2。于是設(shè)計(jì)問題歸結(jié)為找一個(gè) M 次多項(xiàng)式去加權(quán)逼近給定的目標(biāo)函數(shù) D(ω)使最大加權(quán)誤差最小。這一節(jié)的數(shù)學(xué)結(jié)論是交替定理若加權(quán)誤差 E(ω) 在頻帶上有至少 M2 個(gè)極值點(diǎn)并且這些極值點(diǎn)處誤差的符號(hào)正負(fù)交替、幅度都等于同一個(gè)值 δ那么這個(gè)多項(xiàng)式就是唯一的最佳逼近δ 就是能達(dá)到的最小最大誤差。這個(gè)“M2個(gè)交替符號(hào)的極值點(diǎn)”就是最優(yōu)解的指紋。直觀理解一個(gè)有 M 個(gè)自由參數(shù)的逼近誤差曲線如果只有 M1 個(gè)極值你還能調(diào)整參數(shù)把最高的峰壓下去、導(dǎo)致另一個(gè)峰升高直到峰與峰拉平。只有交替峰數(shù)達(dá)到 M2系統(tǒng)才“卡死”——任何微調(diào)都會(huì)讓某個(gè)已經(jīng)等于理論下限的峰更高最優(yōu)性被鎖定了。注意極值點(diǎn)必須同時(shí)覆蓋通帶和阻帶兩個(gè)頻帶的邊界過渡帶邊緣天然是候選極值點(diǎn)。3.2 Remez交換的四個(gè)步驟網(wǎng)格、插值、誤差、換點(diǎn)交替定理給了判據(jù)但怎么找到這組極值點(diǎn)Remez交換算法是一個(gè)不斷“猜測(cè)→驗(yàn)證→修正”的迭代過程初始化在頻率網(wǎng)格上等間隔取 M2 個(gè)點(diǎn)作為初始極值點(diǎn)集合插值假設(shè)這組極值點(diǎn)處誤差都等于 ±δ符號(hào)交替用拉格朗日插值反解出當(dāng)前的多項(xiàng)式系數(shù)和 δ 值求誤差在非常密集的全網(wǎng)格上計(jì)算加權(quán)誤差 E(ω)找出所有局部極值點(diǎn)交換從這些局部極值點(diǎn)中選出絕對(duì)值最大且符合符號(hào)交替要求的 M2 個(gè)點(diǎn)替換掉舊極值點(diǎn)重復(fù) 24直到連續(xù)兩輪的極值點(diǎn)集合不再變化、δ 的波動(dòng)小于閾值。為什么這個(gè)循環(huán)能收斂關(guān)鍵在于第 2 步解出的 δ 一定小于上一輪的真實(shí)最大誤差——因?yàn)樯弦惠喌臉O值點(diǎn)并非真正的等波動(dòng)點(diǎn)插值得到的“假定δ”不是全局最大被低估了。第 3 步在全網(wǎng)格上掃描必然能找到超過這個(gè) δ 的新極值點(diǎn)把它們換進(jìn)來(lái)下一輪的 δ 就會(huì)增大。這個(gè)單調(diào)過程最終收斂到理論最小最大誤差。實(shí)際程序里通常幾輪迭代 δ 就穩(wěn)定了Remez交換的收斂速度在工程上完全夠用。4. 核心代碼拆解C語(yǔ)言實(shí)現(xiàn)里的網(wǎng)格、插值與極值點(diǎn)搜索一個(gè)典型的P-M算法C實(shí)現(xiàn)主循環(huán)看起來(lái)并不長(zhǎng)難在幾個(gè)細(xì)節(jié)函數(shù)。下面這段是主循環(huán)的核心骨架注釋標(biāo)出了每步的機(jī)制。/* remez 主迭代 * L (N-1)/2, N為抽頭數(shù) * ext[] 存極值點(diǎn)在網(wǎng)格上的下標(biāo)大小 L2 */ static int remez_iter(double grid[], double des[], double wt[], int gridsize, int L, int ext[], double *delta, int max_iter) { int iter, i, n_ext, converged 0; for (iter 0; iter max_iter; iter) { /* 1. 由當(dāng)前極值點(diǎn)解切比雪夫方程組算出假定誤差上界 delta * 經(jīng)典實(shí)現(xiàn)不顯式解矩陣改用重心拉格朗日插值一次求出 delta */ *delta calc_delta(ext, L, grid, des, wt); /* 2. 在全網(wǎng)格上計(jì)算誤差函數(shù) err[i] wt[i] * (H(i) - des[i]) */ calc_error_curve(grid, des, wt, err, gridsize, ext, L); /* 3. 掃描全網(wǎng)格找符號(hào)交替、幅度超過 |delta| 的候選極值點(diǎn) */ n_ext find_extrema(ext, grid, err, gridsize, L, *delta); /* 4. 極值點(diǎn)集合不再變化 - 收斂 */ if (n_ext L 2 same_extrema(ext, prev_ext)) { converged 1; break; } /* 5. 交換用新極值點(diǎn)集合覆蓋舊集合進(jìn)入下一輪 */ memcpy(prev_ext, ext, (L 2) * sizeof(int)); } return converged ? iter : -1; }calc_delta承擔(dān)的是交替定理里“假定所有極值點(diǎn)誤差相等”的角色它把當(dāng)前 M2 個(gè)極值點(diǎn)看作已經(jīng)處于等波動(dòng)狀態(tài)用插值倒推多項(xiàng)式在某個(gè)參考點(diǎn)處的值和 δ。這里用重心形式的拉格朗日插值barycentric interpolation而不是直接解矩陣方程因?yàn)楹笳咴跇O值點(diǎn)位置靠近時(shí)矩陣條件數(shù)會(huì)爆炸而重心插值數(shù)值上穩(wěn)得多。find_extrema則是交替定理的“驗(yàn)證器”它只接受符號(hào)正負(fù)交替的點(diǎn)連續(xù)兩個(gè)同號(hào)極值點(diǎn)會(huì)被丟棄因?yàn)樗鼈儾豢赡軜?gòu)成最優(yōu)解。網(wǎng)格密度和分段方式也是機(jī)制的一部分實(shí)現(xiàn)把 0 到 0.5歸一化頻率單位cycles/sample分成各頻帶每個(gè)頻帶內(nèi)部再細(xì)分網(wǎng)格點(diǎn)默認(rèn)每頻帶取約 16 倍的點(diǎn)數(shù)過渡帶邊緣點(diǎn)會(huì)強(qiáng)制加入極值點(diǎn)候選。網(wǎng)格太疏會(huì)漏掉真正的局部極值點(diǎn)導(dǎo)致迭代停滯在錯(cuò)誤解太密則純?cè)鲇?jì)算量。這是這個(gè)算法里“工程直覺”最集中的地方——數(shù)學(xué)上要求掃描所有頻率工程上只能在足夠密的離散網(wǎng)格上近似。5. 避坑Remez迭代翻車的幾個(gè)典型現(xiàn)場(chǎng)1. 迭代不收斂delta 上下抖動(dòng)現(xiàn)象程序跑滿 max_iter 次誤差曲線仍然不是等波動(dòng)delta 在相鄰輪次之間反復(fù)跳變。原因最常見的是網(wǎng)格太稀疏局部極值點(diǎn)被漏掉導(dǎo)致 find_extrema 返回的候選點(diǎn)集合不穩(wěn)定也可能是初始極值點(diǎn)個(gè)數(shù)不對(duì)比如抽頭數(shù)給錯(cuò)導(dǎo)致 L 與網(wǎng)格大小不匹配。解決把每頻帶的網(wǎng)格密度從 16 倍翻到 32 倍再跑同時(shí)檢查調(diào)用端傳入的抽頭數(shù)和頻帶數(shù)組是否一致尤其是奇數(shù)抽頭數(shù)Type I不要誤傳成偶數(shù)。2. 阻帶衰減始終達(dá)不到指標(biāo)誤差曲線“一邊高一邊低”現(xiàn)象設(shè)計(jì)結(jié)果通帶波紋很小但阻帶衰減差得遠(yuǎn)或者反過來(lái)阻帶很好通帶超差。原因加權(quán)函數(shù) W(ω) 沒有正確表達(dá)你的偏好。P-M算法最小化的是加權(quán)誤差不是絕對(duì)誤差——你不給阻帶加權(quán)重算法就默認(rèn)阻帶誤差和通帶誤差“同罪”阻帶目標(biāo)值接近0同樣的絕對(duì)誤差在阻帶看來(lái)占比反而小算法自然懶得照顧阻帶。解決按“阻帶衰減需要比通帶波紋高多少dB”來(lái)設(shè)定權(quán)重比例。例如要通帶波紋 0.1dB、阻帶衰減 60dB兩者誤差比約 1:1000權(quán)重就要反過來(lái)設(shè)約 1000:1。3. 系數(shù)出現(xiàn)不合理的異常大值濾波器輸出發(fā)飄現(xiàn)象h[n] 序列里某些抽頭值有幾十甚至上百的數(shù)量級(jí)明明想要的通帶增益只是 1。原因極值點(diǎn)之間靠得太近尤其是把過渡帶邊緣兩點(diǎn)同時(shí)列為強(qiáng)制極值點(diǎn)時(shí)插值矩陣接近奇異數(shù)值解被放大。解決檢查過渡帶邊緣點(diǎn)是否被重復(fù)加入候選使用重心插值版本必要時(shí)提高網(wǎng)格密度讓相鄰極值點(diǎn)間距不小于一個(gè)網(wǎng)格步長(zhǎng)。4. 在 fs/2 處響應(yīng)不為零設(shè)計(jì)結(jié)果和預(yù)期偏差很大現(xiàn)象用偶數(shù)抽頭數(shù)Type II設(shè)計(jì)高通或帶阻濾波器Nyquist頻率處頻率響應(yīng)不為零甚至很差。原因Type II 對(duì)稱FIR在 ωπ 處幅度響應(yīng)被結(jié)構(gòu)強(qiáng)制為零數(shù)學(xué)上不可能在那一點(diǎn)給出非零目標(biāo)——再迭代也突破不了這個(gè)結(jié)構(gòu)性約束。解決設(shè)計(jì)高通、帶阻這類在 Nyquist 處要求非零響應(yīng)的濾波器時(shí)強(qiáng)制使用奇數(shù)抽頭數(shù)Type I或者在指標(biāo)里避免在π處定義非零目標(biāo)。5. 過渡帶邊緣位置總差幾個(gè)采樣點(diǎn)現(xiàn)象測(cè)出來(lái)的 -3dB 點(diǎn)、阻帶起始頻率和設(shè)計(jì)指標(biāo)對(duì)不上位置整體偏移。原因頻帶數(shù)組的定義方式問題。實(shí)現(xiàn)里通帶和阻帶的邊界必須精確落在網(wǎng)格點(diǎn)上如果換算歸一化頻率時(shí)用了錯(cuò)誤的采樣率或者過渡帶邊緣沒有對(duì)齊網(wǎng)格算法只能近似移動(dòng)邊界。解決確認(rèn)頻帶數(shù)值是按“cycles/sample”還是“rad/sample”計(jì)算的設(shè)計(jì)前先打印網(wǎng)格落點(diǎn)確認(rèn)每個(gè)頻帶邊緣對(duì)應(yīng)的網(wǎng)格下標(biāo)是整數(shù)。6. 設(shè)計(jì)取舍階數(shù)、波紋、過渡帶以及“最優(yōu)”并不等于“好用”P-M算法給出了數(shù)學(xué)意義上的最優(yōu)但“最優(yōu)”是相對(duì)你給的權(quán)重和階數(shù)而言的。權(quán)重如何分配、階數(shù)取多少這些選擇里藏著真正的工程判斷。首先階數(shù)是最大的取舍變量。階數(shù)越高逼近多項(xiàng)式的次數(shù)越高理論能達(dá)到的最小最大誤差 δ 越小但代價(jià)是指標(biāo)收益遞減且數(shù)值穩(wěn)定性變差——高次多項(xiàng)式的極值對(duì)系數(shù)微小擾動(dòng)極度敏感C語(yǔ)言里 double 的舍入誤差會(huì)開始主導(dǎo)設(shè)計(jì)結(jié)果。實(shí)際經(jīng)驗(yàn)是階數(shù)一到幾百盡量用浮點(diǎn)設(shè)計(jì)加量化分析兩步走。其次等波動(dòng)結(jié)構(gòu)本身是一種偏好。它讓最壞誤差最小但代價(jià)是誤差在整個(gè)頻帶上均勻分布。如果你需要的是“通帶極度平坦阻帶只要夠用”等波動(dòng)設(shè)計(jì)反而浪費(fèi)了阻帶的裕量——這時(shí)候用最小二乘設(shè)計(jì)誤差能量最小得到的通帶更平只是沒有嚴(yán)格的最壞誤差保證。P-M算法的“最優(yōu)”只在極小極大準(zhǔn)則下成立換一個(gè)目標(biāo)函數(shù)最優(yōu)解就換人了。第三加權(quán)函數(shù)是設(shè)計(jì)者唯一表達(dá)偏好的地方但權(quán)重比和實(shí)際波紋比不是線性對(duì)應(yīng)。你設(shè)阻帶權(quán)重是通帶的 100 倍實(shí)際阻帶衰減和通帶波紋的差并不會(huì)正好是 40dB因?yàn)檫^渡帶寬度、極值點(diǎn)分布都會(huì)影響最終結(jié)果。正確做法是先跑一輪不帶權(quán)重的設(shè)計(jì)量出各頻帶實(shí)際誤差再反推權(quán)重比迭代兩三次才能對(duì)齊指標(biāo)。我自己現(xiàn)在的習(xí)慣是設(shè)計(jì)一個(gè)新濾波器時(shí)先用一個(gè)寬松階數(shù)跑通P-M程序打印每輪 delta 的收斂軌跡和最終極值點(diǎn)符號(hào)序列確認(rèn)交替定理的條件確實(shí)滿足了再按指標(biāo)去收緊階數(shù)。這樣踩過幾次上面說(shuō)的坑之后我再也沒讓一個(gè)看起來(lái)“不收斂”的設(shè)計(jì)白白浪費(fèi)過時(shí)間。這個(gè)驗(yàn)證習(xí)慣幫我少走了很多彎路希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取