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