プチコン 非公式コミュニティ
トピック
こういち
◆ou0jbJnEJ0Kb
2025/9/2 7:39
2025/9/2 7:39
質問
最小位相FIRフィルタの設計方法
なんとなく自己流で最小位相FIRフィルタっぽいものを作ってみたんですけど、実際の特性が元のFIRフィルタとは違うものになってしまいます。(伝達関数の根のうち、絶対値が1を超えるものはその根を削除して、逆数を新たな根として変換してました)
ヒルベルト変換というものを使うと近似的に最小位相フィルタが求まるっぽいのですが、ヒルベルト変換がよく分かりません。
そして、厳密な最小位相フィルタを設計したいんですけど、それは可能でしょうか。
あとついでに変形元フィルタの設計に使われるParks-McClellanのアルゴリズムについても知っておきたい。
そういえば、任意長FFTって持ってたっけな…(確かプチコンだと速度出せなそうで作ってなかった気がする)
ヒルベルト変換というものを使うと近似的に最小位相フィルタが求まるっぽいのですが、ヒルベルト変換がよく分かりません。
そして、厳密な最小位相フィルタを設計したいんですけど、それは可能でしょうか。
あとついでに変形元フィルタの設計に使われるParks-McClellanのアルゴリズムについても知っておきたい。
そういえば、任意長FFTって持ってたっけな…(確かプチコンだと速度出せなそうで作ってなかった気がする)
コメント
こういち
2025/9/2 7:47 ネタバレ
◆ou0jbJnEJ0Kb
こういち
2025/9/2 7:48 ネタバレ
◆ou0jbJnEJ0Kb
こういち
2025/9/2 7:49 ネタバレ
◆ou0jbJnEJ0Kb
こういち
2025/9/2 7:50 ネタバレ
◆ou0jbJnEJ0Kb
こういち
2025/9/9 7:43
◆ou0jbJnEJ0Kb
まず、Remezのアルゴリズムについて理解しようと思います。
https://ja.wikipedia.org/wiki/Remez%E3%81%AE%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0
本来はある関数(平方根や三角関数など)を多項式で近似するアルゴリズムらしい。
難しいのは、連立方程式を解く部分と極値を求める部分ぐらい。
https://ja.wikipedia.org/wiki/Remez%E3%81%AE%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0
本来はある関数(平方根や三角関数など)を多項式で近似するアルゴリズムらしい。
難しいのは、連立方程式を解く部分と極値を求める部分ぐらい。
こういち
2025/9/9 7:43
◆ou0jbJnEJ0Kb
連立方程式はガウスの消去法を使うとO(N^3)で解ける(ちなみに、ガウスの消去法よりハウスホルダー変換の方が計算量同じで精度良いらしい)けど、ニュートン補間という方法を使うことで、O(N^2)で精度良く計算できるっぽい。
ニュートン補間
https://risalc.info/src/newtons-interpolation.html#ex
N+1点を通るN次多項式を求めるアルゴリズムらしい。
言い換えると
係数行列がヴァンデルモンド行列という形の連立方程式を高速に解くアルゴリズム(Levinson-Durbin法で解けるのはテプリッツ行列で別の形っぽい)
こんなのあったんだ…
ニュートン補間
https://risalc.info/src/newtons-interpolation.html#ex
N+1点を通るN次多項式を求めるアルゴリズムらしい。
言い換えると
係数行列がヴァンデルモンド行列という形の連立方程式を高速に解くアルゴリズム(Levinson-Durbin法で解けるのはテプリッツ行列で別の形っぽい)
こんなのあったんだ…
こういち
2025/9/9 7:44
◆ou0jbJnEJ0Kb
Remez法で解きたい連立方程式はヴァンデルモンド行列ではないけど、2つのヴァンデルモンド行列に分解することで高速に解けるらしい。
あとは、極値を求める。ちょい面倒ぐらい。
最急降下法や三分探索などが使える。微分が求まるなら最急降下法が楽。
符号が交互に変わるので、極値はXi-1〜Xi+1の間に必ず存在する。
ただ、極値が1つとは限らないので、三分探索は厳しくて、最急降下法を使う必要がある。
n次関数は容易に微分可能なので、元関数の微分を渡せばOK
n次関数はn-1個の極値を持つので、X0,Xn+1は固定で極値を探すのはx1〜Xnのn-1個探せば良いはず。
この極値が求まらないといけないので、おそらくRemez法はノイズに弱い。
あとは、極値を求める。ちょい面倒ぐらい。
最急降下法や三分探索などが使える。微分が求まるなら最急降下法が楽。
符号が交互に変わるので、極値はXi-1〜Xi+1の間に必ず存在する。
ただ、極値が1つとは限らないので、三分探索は厳しくて、最急降下法を使う必要がある。
n次関数は容易に微分可能なので、元関数の微分を渡せばOK
n次関数はn-1個の極値を持つので、X0,Xn+1は固定で極値を探すのはx1〜Xnのn-1個探せば良いはず。
この極値が求まらないといけないので、おそらくRemez法はノイズに弱い。
こういち
2025/9/9 7:45
◆ou0jbJnEJ0Kb
ニュートン補間の実装
まずは商差。
とりあえず、再帰で実装
DEF DIVIDED(X[],Y[],S,N)
IF N>=1 THEN
RETURN (DIVIDED(X,Y,S+1,N-1)-DIVIDED(X,Y,S,N-1)))/(X[S+N]-X[S])
ELSE
RETURN Y[S]
ENDIF
END
これをまともに計算すると計算量がO(2^N)になるので、動的計画法に直す。
そのままS,Nを添え字とするDPテーブルを作ると良い。
漸化式もそのまま
DP[S,N]=(DP[S+1,N-1]-DP[S,N-1])/(X[S+N]-X[S])
でやるだけ。
DIM DP#[100,100]
DEF DIVIDED(X[],Y[],S,N)
VAR I%,J%
FOR I%=0 TO N+S
DP#[I%,0]=Y[I%]
NEXT I%
FOR I%=1 TO N
FOR J%=0 TO N+S-I%
DP#[J%,I%]=(DP#[J%+1,I%-1]-DP#[J%,I%-1])/(X[J%+I%]-X[J%])
NEXT J%
NEXT I%
RETURN DP#[S,N]
END
実はNの次元は要らなかったりする。
まずは商差。
とりあえず、再帰で実装
DEF DIVIDED(X[],Y[],S,N)
IF N>=1 THEN
RETURN (DIVIDED(X,Y,S+1,N-1)-DIVIDED(X,Y,S,N-1)))/(X[S+N]-X[S])
ELSE
RETURN Y[S]
ENDIF
END
これをまともに計算すると計算量がO(2^N)になるので、動的計画法に直す。
そのままS,Nを添え字とするDPテーブルを作ると良い。
漸化式もそのまま
DP[S,N]=(DP[S+1,N-1]-DP[S,N-1])/(X[S+N]-X[S])
でやるだけ。
DIM DP#[100,100]
DEF DIVIDED(X[],Y[],S,N)
VAR I%,J%
FOR I%=0 TO N+S
DP#[I%,0]=Y[I%]
NEXT I%
FOR I%=1 TO N
FOR J%=0 TO N+S-I%
DP#[J%,I%]=(DP#[J%+1,I%-1]-DP#[J%,I%-1])/(X[J%+I%]-X[J%])
NEXT J%
NEXT I%
RETURN DP#[S,N]
END
実はNの次元は要らなかったりする。
コメントを書く
- こちらは「プチコン3号」「プチコンBIG」など、プチコンシリーズに関する話題を扱ったコミュニティです
- プチコンシリーズにまったく関係ない書き込みはご遠慮下さい。削除の対象となります
- こちらにはその他のゲームや雑談のコミュニティはなく、作る予定もありません (ひとりで管理できないため)。ごめんなさい
- ユーザー登録なしで書き込みができます
- 秘密の合い言葉は成りすましの防止 (トリップ機能)、書き込みの編集時の本人認証に使用します
- 秘密の合い言葉に他人に推測されやすい言葉、他サービスと同じパスワードは入力しないでください。
- 書き込むと、投稿時に入力したお名前と秘密の暗号が記憶され、ログイン状態になります