【デジタルフィルタの基礎１】まずは、FIRフィルタの世界を一望してみよう

Finite Impulse Response（有限インパルス応答）、略して FIR フィルタとは、入力に対する応答（インパルス応答）が有限の長さで確実にゼロになるデジタルフィルタのことを指します。たとえば、鐘を叩いたときに「カーン」と響く音を思い浮かべてください。これは「叩いた瞬間（入力）」に対して「余韻（応答）」がしばらく続く現象です。この“余韻の長さ”に相当するものが、信号処理ではインパルス応答と呼ばれます。FIR フィルタの特徴は、このインパルス応答が有限時間のうちに完全に消えることです。

　それに対し、応答がずーっとしつこく、時間がたつと小さくはなるけど、無限に続くものは、Infinite Impulse Response（IIR：無限インパルス応答）フィルタと呼ばれます。多くの学生が、深く考えずに「とりあえずバターワース型の IIR フィルタ」を使ってしまうことがあります。IIR は強力で便利な道具ですが、位相歪みが避けられない、数値安定性に注意が必要といった側面もあります。

　デジタルフィルタの理屈を知らず、実際の計算式を見たことがない人は、FIRや、IIRという用語がなんだか難しそうに響くかもしれまっせん。しかし、フィルタリングで行う計算式はとてもシンプルです。

　FIR フィルタは、入力信号を $\displaystyle{ x[n] }$ 、出力信号を $\displaystyle{ y[n] }$ として、

$\displaystyle{ y[n] = \sum_{k} b_k\, x[n-k] }$

という式で表されます。ここで、 $\displaystyle{ \{ b _ k \} }$ は、フィルタを表す係数のセットです。デジタルフィルタの計算式では、高校のときに習った記号と掛け算と足し算しか使ってません。

　また、IIR フィルタは、

$\displaystyle{ y[n] = \sum_{k} b_k\, x[n-k] \;+\; \sum_{k} a_k\, y[n-k] }$

という式で表されます。ここで、 $\displaystyle{ \{ a _ k \} }$ も、フィルタを表す係数のセットです。右辺にも $\displaystyle{ y[n] }$ の過去の値があることがFIRとの違いです。

　どちらも複雑な計算式ではありません。「掛け算して足す」という、ごく基本的な計算の繰り返しで動いています。これらの計算の仕組みを知ると、フィルタは驚くほど身近で扱いやすいものに見えてくるはずです。

　IIRについても理解してほしいですが、今回の主役はFIRです。

　この記事の最初に、インパルス応答を鐘の音に例えて導入しましたが、実際のインパルス応答とは、「単位インパルスを入力したときの出力」のことです。単位インパルスとは、時間 $\displaystyle{ n = 0 }$ で 1、それ以外は 0 の信号

$\displaystyle{ \delta[n] = \begin{cases} 1, & n = 0 \\ 0, & n \neq 0 \end{cases} }$

のことです (下図)。

この信号を FIR フィルタに入力すると、出力はどうなるかというと、

$\displaystyle{ y[n] = \sum_{k} b_k \, \delta[n-k] }$

において、 $\displaystyle{ \delta[n-k] }$ は $\displaystyle{ n = k }$ のときだけ 1 になるので、

$\displaystyle{ y[n] = b_n }$

となります。つまり、係数 $\displaystyle{ \{ b _ n \} }$ をそのまま順に並べたものがインパルス応答です。係数はインパルス応答と同じものですが、あらためて、インパルス応答を $\displaystyle{ h[n] }$ と表すことにすれば、

$\displaystyle{ h[n] = b_n }$

であり、FIR フィルタの式は、インパルス応答を使って

$\displaystyle{ y[n] = \sum_{k} h[k]\, x[n-k] }$

と表されるとも言えます。有限インパルス応答 (FIR) とは、有限の $\displaystyle{ k }$ の範囲でのみ、 $\displaystyle{ h[k] \neq 0 }$ であるということです。

FIRフィルタのインパルス応答 $\displaystyle{ h[k] }$ の例。ここに描かれていない横軸の領域で、 $\displaystyle{ h[k] }$ がすべてゼロになるのがFIRフィルタ。

1. なぜ「有限である」ことが重要なの？

　インパルス応答が有限で終わるという性質は、実はとても重要で、次のような大きなメリットを生みます。ここでは、IIR フィルタではどのような問題が起こりやすいか、という比較をしながら見ていきます。

① 絶対に安定する

　FIR フィルタはインパルス応答が有限長です。したがって、有限個の点についてしか、畳み込み演算を行いません。そのため、原理的に不安定になることはありません。

　このことは、次のような数学的事実に基づきます。それは、

入力時系列 $\displaystyle{x[n]}$ が有限分散（ $\displaystyle{\mathbb{E}[x[n]^ 2] \lt \infty}$ ）
FIR フィルタのインパルス応答 $\displaystyle{h[n]}$ が有限長（または絶対可和が、 $\displaystyle{\sum _ n |h[n]| \lt \infty}$ ）

であれば、出力 $\displaystyle{y[n] = \sum _ k h[k] x[n-k]}$ の分散も必ず有限になるということです。

　一方、IIR フィルタは、極（pole）の位置が 1 の外側に出ると発散するため、パラメータのわずかなズレでも「暴走」「発散」「ノイズ増幅」が起こりえます。さらに有限精度の計算では、丸め誤差や量子化誤差が極の近傍で増幅され、安定性が大きな問題になります。

　このように、FIR フィルタは

原理的に発散しない（Bounded-Input Bounded-Output Stability (BIBO) 安定）
数値計算でも極が無いため誤差に強い
有限分散の入力に対して、出力も常に有限分散

という意味で、安定性を気にしなくてよい非常に安全なフィルタです。

② 設計しやすく予測しやすい

　FIR のインパルス応答は有限長の係数列なので、「フィルタの形が直感的に理解しやすい」、「設計後の特性が予測しやすい」、「最適化や制約条件の設定が容易」という利点があります。

一方、IIR は、「極と零点の配置」、「安定条件」、「位相歪み」などを同時に考える必要があり、総合的な制御が難しいという側面があります。

③ 完全な線形位相が実現できる

　FIR フィルタは、インパルス応答を対称または反対称に設計することで、理論的に完全な線形位相（群遅延一定）を実現できます。ここでいう線形位相とは、フィルタの位相特性が周波数に対して

$\displaystyle{ \phi(f) = - 2 \pi f D }$

ということです。

　この特性に注目し、フーリエ変換の「時間シフトの性質」を使うと、この意味がはっきりします。

連続時間の場合、 $\displaystyle{x(t)}$ のフーリエ変換を $\displaystyle{X(f)}$ とすれば、

$\displaystyle{ x(t-D) \;\Longleftrightarrow\; e^{-j 2\pi f D}\, X(f) }$

が成り立ちます。

　離散時間の場合も同様に、 $\displaystyle{x[n]}$ のフーリエ変換を $\displaystyle{X(f)}$ とすれば、

$\displaystyle{ x[n-D] \;\Longleftrightarrow\; e^{-j 2\pi f D}\, X(f) }$

が成り立ちます。つまり、スペクトルに、 $\displaystyle{ e^{-j 2\pi f D} }$ を掛ける操作は、時間領域では「波形をそのまま $\displaystyle{ D }$ だけ遅らせる」ことに対応します。

　線形位相 FIR フィルタの周波数特性を

$\displaystyle{ H(f) = |H(f)|\, e^{-j 2\pi f D} }$

と書くと、入力 $\displaystyle{X(f)}$ に対する出力は

$\displaystyle{ Y(f) = H(f)\, X(f) = |H(f)|\, e^{-j 2\pi f D}\, X(f) }$

となります。時間領域に戻すと、

$\displaystyle{|H(f)|}$ ：振幅だけを周波数ごとに調整する（＝平滑化や帯域制限）
$\displaystyle{e^ {-j 2 \pi f D}}$ ：波形全体を $\displaystyle{D}$ サンプルだけシフトする

に対応するので、通過帯域内で $\displaystyle{ |H(f)| \approx 1 }$ なら、出力は「少しなめらかになった元の波形を、そっくりそのまま $\displaystyle{ D }$ だけ遅らせたもの」とみなすことができます。したがって、フィルタを通しても入力波形の「形」はほとんど歪まず、時間方向に平行移動したような出力が得られます。

　「心電図の QRS 複合波の形を崩したくない」「光電脈波の立ち上がりを正確に扱いたい」「音響でドラムのアタックを変えたくない」といったように、波形の形そのものが情報になっている信号では、線形位相 FIR を使うことが決定的なメリットになります。

　一方 IIR フィルタでは、一般に位相 $\displaystyle{ \phi(f) }$ が周波数に対して直線にならず（非線形位相）、群遅延 $\displaystyle{ \tau _ g(f) = -d\phi/df }$ が周波数によって大きく変化します。その結果、

ある周波数成分だけ極端に遅れる
波形のピーク位置が前後にずれる
立ち上がり／立ち下がりの非対称性が生じる

といった、形状の変形が起こり、波形の解釈に致命的な影響を与えることがあります。

④ 安全で扱いやすい

　FIR は極を持たないため、数値誤差が蓄積しても挙動が大きく変わりません。

対して IIR は、「丸め誤差で極が微妙に動く」、「フィルタ特性が劣化する」、「特定のデータでは発振が起こる」など、実装上のデリケートさがあります。

FIRフィルタである移動平均フィルタの例 (Rスクリプト:1)。上段が入力信号、それ以外の段の青実線が出力信号、桃色実線が入力信号に含まれる確定信号。

2. FIR フィルタはどのように働くの？

　FIR フィルタは数学的には「入力信号を一定区間にわたって重み付きで足し合わせる（畳み込む）」操作で定義されます。その基本となるのが畳み込みの式です。

■ 一般の（非因果を含む）FIR フィルタの定義

　FIR フィルタの最も一般的な形は、インパルス応答 $\displaystyle{ h[n] }$ が有限長であることだけを条件として、次のように定義できます。

$\displaystyle{ y[n] = \sum_{k = n_1}^{n_2} h[k]\; x[n - k] }$

ここで、

$\displaystyle{h[k]}$ ：インパルス応答（フィルタ係数）
$\displaystyle{n _ 1, n _ 2}$ ：有限の整数で、 $\displaystyle{h[k]}$ がゼロでない範囲を指定しています。

　この形では、インパルス応答 $\displaystyle{h[k]}$ に正の添字（未来側）が含まれると、和の中に、未来のサンプル $\displaystyle{ x[n+k] }$ （ $\displaystyle{ k\gt0 }$ ）が登場することになります。そのようなフィルタは、非因果（non-causal）と呼ばれます。

　非因果フィルタは、数学的には完全に定義できても、未来の値を必要とするため、リアルタイム処理には使えません。たとえば $\displaystyle{ h[+2] }$ が非ゼロなら、

$\displaystyle{ y[n] \ \text{を求めるのに} \ x[n+2] \ \text{が必要} }$

となり、「まだ観測されていない未来の信号を使って出力を計算する」という状況が生じるからです。

　しかし、非因果性は

オフライン解析（後から処理する場合）
データをまとめて処理する場合
線形位相 FIR の理論を考える場合

には何の問題もありません。

■ 因果型FIR フィルタ

　リアルタイム処理では未来の値を使えないため、そのような場合には因果型のインパルス応答使います。因果型 FIR フィルタは一般に、次の形で表されます。

$\displaystyle{ y[n] = \sum_{k=0}^{N} h[k]\, x[n-k] }$

因果型、非因果型は用途に応じて使い分けてください。

3. FIR フィルタでできること

　FIR フィルタはこの「重み付け平均」を巧みに設計することで、次のような処理が実現できます。

ノイズ除去（平滑化）　高周波成分を抑えて、信号をなめらかにすることができます。
特定の周波数成分の抽出　FIRでも、ローパス、ハイパス、バンドパスができます。
不要成分の抑制（バンドストップ）　電源ノイズ（50/60 Hz）だけを取り除くノッチフィルタにもなります。
波形の形を崩さない処理（線形位相）　生体信号や音響では波形の形を維持することが大切です。

　このように、FIR フィルタは「単なる平均化処理」を超え、係数の選び方次第で、精密な周波数操作が自在に行えるツールになります。

理想ローパスフィルタであるsincを有限長に切り取る方法で作成されたFIRローパスフィルタの例 (Rスクリプト:2)

4. FIR フィルタの世界をどう理解するか

　FIR フィルタは、位相特性・周波数特性・設計法・構造の側面により、その特徴が分類できます。これらを整理して理解するために、この記事では次の4 つの軸で体系化します。

位相特性：インパルス応答の対称性（Type I〜IV）
周波数特性：LPF、HPF、BPF、Notch、微分器、Hilbert など
設計法：窓関数法、等リプル法、最小位相設計、周波数サンプリング法
構造：直接形、多相フィルタ、線形位相構造、FFT畳み込み

　これによって、FIR フィルタがどのような存在で、どのように作られ、どのように使われるのかが一気に見渡せるようになります。

■ 位相特性から見る FIR フィルタの本質

　FIR フィルタの最大の特長は、完全な線形位相（群遅延が一定）を実現できることです。この線形位相特性は、インパルス応答の対称性によって保証され、FIR フィルタは次の 4 種類に分類されます。

　以下の説明で登場する FIR フィルタの「次数」とは、インパルス応答の長さ（タップ数）から 1 を引いたものを意味します。たとえば、インパルス応答の長さが 5 なら次数は 4、長さが 4 なら次数は 3 になります。

　なぜ、次数はインパルス応答の長さより 1 だけ小さいのかというと、インパルス応答を z 変換して得られる

$\displaystyle{ H(z) = h[0] + h[1]z^{-1} + \cdots + h[N]z^{-N} }$

が $\displaystyle{ z^ {-1} }$ に関する N 次多項式になるためです。多項式の次数が係数の個数より 1 小さくなるのと同じ理由で、FIR フィルタの次数もインパルス応答の長さより 1 小さくなるのです。

Type I：偶数次数・対称　たとえば、 $\displaystyle{ h = [1, 3, 5, 3, 1] }$ のように、奇数長で左右対称のパターン。最も汎用的で、ローパスからハイパス、バンドパスやノッチまで幅広く使える。
Type II：奇数次数・対称　たとえば、 $\displaystyle{ h=[1,2,2,1] }$ のように、偶数長で左右対称のパターン。 $\displaystyle{\omega = \pi}$ で必ず 0 になり、高域通過フィルタには向かない。
Type III：偶数次数・反対称　たとえば、 $\displaystyle{ h = [1,2,0,−2,−1] }$ のように、奇数長で左右反対称のパターン。直流成分が必ず 0 になる。微分器や高速ハイパスフィルタに適する。
Type IV：奇数次数・反対称　たとえば、 $\displaystyle{ h = [1,1,−1,−1] }$ のように、偶数長で左右反対称のパターン。Hilbert 変換器や分数遅延フィルタに使われる特殊タイプ。

　この分類は、「何が設計できるか」を決める最も本質的な基盤であり、FIR フィルタの設計は、まず Type を決めることから始まります。

■ 周波数特性から見るフィルタの機能

　FIR フィルタは、周波数領域においてさまざまな機能を実現できる柔軟な構造を持っています。ここでは、代表的な 7 種類の FIR フィルタとその用途を紹介します。

低域通過フィルタ（Low-Pass Filter） 高周波成分を抑えることで、ノイズ除去や信号の平滑化に用います。
高域通過フィルタ（High-Pass Filter） 低周波成分を取り除くことで、トレンドの除去やベースラインドリフト補正に役立ちます。
帯域通過フィルタ（Band-Pass Filter） 特定の周波数帯域のみを通過させるフィルタで、特定周波数の成分抽出に有効です。
帯域阻止フィルタ（Band-Stop Filter / Notch Filter） 特定の狭い周波数帯のみを抑制するフィルタで、50 Hz / 60 Hz の電源ノイズ除去に広く使用されます。
微分フィルタ（Differentiator） 信号の変化を強調し、急峻な立ち上がりやピークの検出に用います。
ヒルベルト変換フィルタ（Hilbert Transformer） 信号の振幅を変えず、全周波数成分の位相を一律に 90° シフトさせるフィルタで、解析信号生成や瞬時位相解析に不可欠です。
分数遅延フィルタ（Fractional Delay Filter） 信号を 0.1 サンプル、0.5 サンプルなど、非整数サンプル数だけ遅延させるためのフィルタで、音響、通信、波形整形などで重要な役割を果たします。

　以下に、これらの FIR フィルタのインパルス応答例と、その周波数特性の図を示します。これまで、深く考えずに IIR フィルタで処理していた方も、FIR フィルタでここまで多様な特性を実現できることに驚くはずです。

低域通過フィルタ（Low-Pass Filter, Hamming-windowed）。左： sinc 形のインパルス応答に Hamming 窓を掛けたもので、中心付近が最も大きく、外側に行くほど減衰していく典型的な非因果 LPF の形。右：低周波が通過し、高周波成分が滑らかに抑制される低域通過特性が見られる。通過帯域は平坦で、カットオフ付近で緩やかに減衰している。

高域通過フィルタ（High-Pass Filter, Spectral Inversion）。左：低域通過フィルタのインパルス応答をスペクトル反転したもので、中心成分が負になり、両端に正の振幅が現れる特徴的な形。右：低周波側が強く減衰し、高周波側が通過する高域通過特性。カットオフ周波数より上の成分が平坦に残る。

帯域通過フィルタ（Band-Pass Filter, LPF の差分）。左：低カットオフ LPF と高カットオフ LPF の差分によるインパルス応答。0 付近で符号変化があり、外側に振動を持つ帯域通過特有の形状。右：指定した帯域のみが通過し、その外側（低周波・高周波）が抑制される帯域通過特性。通過帯域が中央に山として現れる。

帯域阻止フィルタ（Band-Stop / Notch Filter, 電源ノイズ除去用）。左：特定帯域だけを打ち消すため、中央に鋭いゼロクロッシングが現れる複雑なインパルス応答。右：中央の周波数付近で鋭く深い減衰が現れ、それ以外の周波数はほぼ通過する典型的なノッチ特性。

微分フィルタ（Differentiator FIR, 反対称）。左：中心が 0、左右で大きさが反対符号になる反対称インパルス応答。Hamming 窓によって外側が滑らかに抑制されている。右：振幅が周波数に比例して増加する微分特性が見られる。低周波は抑制され、高周波でゲインが伸びる。

Hilbert 変換フィルタ（Hilbert Transformer, Type III）。左：インパルス応答は完全な反対称で、奇数インデックスにのみ値を持つ Hilbert 変換器特有の形状。中心は必ず 0。右：振幅は周波数全体でほぼ一定で、位相特性がすべての周波数で +90°（または −90°）になる Hilbert 変換に対応する振幅応答。

分数遅延フィルタ（Fractional Delay Filter, Lagrange 型）。左： Lagrange 補間による滑らかな係数列で、整数サンプルでは表現できない「0.5 サンプル遅延」を再現するための形状になっている。右：振幅はほぼ平坦で、位相が周波数に対して直線的に変化する特性（＝分数サンプル遅延の実現）を反映した周波数応答。

■ 設計手法：どのように FIR を作るか

　FIR フィルタは設計自由度が高く、目的に応じて複数の方法が選べます。

（1）窓関数法（Window Method）— 最も簡単で実用的

　理想フィルタに窓を掛けて有限長化する方法。Hamming、Blackman、Kaiserなどがよく使われる。

長所：安定・設計容易
短所：遷移帯域の制御が弱い

（2）等リプル設計（Parks–McClellan 法）— 最適なフィルタ長

　Chebyshev 近似を用いて、周波数領域の誤差を等リプル化する方法。

長所：最短長で仕様を満たす「最適フィルタ」
短所：設計が複雑

（3）周波数サンプリング法 — 周波数軸から直接設計

　周波数応答をサンプリングして IDFT。直感的で柔軟。

（4）最小位相フィルタ — 遅延を最小にしたいとき

　線形位相を捨て、群遅延を最小化。音響やリアルタイムデバイスで有効。

（5）最適化ベース設計

　L1/L2 最小化や制約付き最適化でフィルタ特性を最適化。レーダー・医療など特殊用途で使用。

■ 構造：どのように FIR を実装するか

　FIR フィルタは構造の工夫によって計算量や性能が大きく変わります。

（1）直接形（Direct Form）

　最もオーソドックス。実装が簡単でハードウェアにも適する。

（2）線形位相構造（対称性を利用）

　対称性（ $\displaystyle{h[n] = h[N–n]}$ を利用して掛け算を半分に削減。DSP・FPGA では標準手法。

（3）多相フィルタ（Polyphase Structure）

　サンプリングレート変換（アップ/ダウン）で必須。オーディオ圧縮・画像処理・通信で重要な技術。

（4）FFT畳み込み（Overlap-Add/Save）

　長い FIR（数千点以上）を高速処理するための技術。音響畳み込みリバーブや画像フィルタで必須。

5. まとめ：深く理解する人だけが見える景色がある

　最近は、「目的さえ達成できればよい」、「道具が使えれば中身は知らなくていい」という風潮が強まっている気がします。高校生までの勉強で言えば、使える知識、役に立つ知識を身に着けるのではなく、大学入試で得点をとれるテクニックを優先する感じです。フィルタ設計でも、ブラックボックスの関数を呼び出せば結果が出るため、大学入学後は、数学を避け、理屈を避け、手間のかかる理解そのものを省いてしまう人も少なくありません。

　もちろん、やりたくないことに無理に深入りする必要はありません。誰にでも限られた時間があり、優先すべきことは人それぞれです。「とりあえず動けばよい」という選択も、場面によっては合理的でしょう。

　しかし一方で、仕組みを理解した人にだけ見える世界があることも、ぜひ覚えておいてください。FIR フィルタをはじめとする信号処理の理論は、単なる計算の手順ではなく、「なぜそうなるのか」「どのように動いているのか」という深い構造と美しい整合性の上に成り立っています。

　もしあなたが興味を持ち、「もっと知りたい」「自分でも理解してみたい」と思ったなら、その気持ちを大切にしてください。そして、学びの歩みを進めてください。深く学ぶことは、遠回りに見えて、実は最も「近道」かもしれません。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

[付録] Rでやってみる

Rスクリプト1：移動平均フィルタ

########################################
# Comparison of Moving Average Filters
# Input signal: Gaussian peak train + noise
########################################

### 1. Generate Input Signal ---------------------------------------

n.peaks <- 10
sig     <- exp(seq(log(200), log(2), length = n.peaks))
mu      <- max(sig) * 10 * (1:n.peaks)
dt      <- 1

t.range <- seq(0, max(mu) + max(sig) * 5, by = dt)
n.t     <- length(t.range)

# Clean signal
x.clean <- numeric(n.t)
for (i in 1:n.peaks) {
  x.clean <- x.clean +
    dnorm(t.range, mean = mu[i], sd = sig[i]) /
    dnorm(mu[i],    mean = mu[i], sd = sig[i])
}

# Noisy signal
x    <- x.clean + rnorm(n.t, sd = 0.05)
time <- seq(0, (n.t - 1) * dt, by = dt)


### 2. Define Moving Average Filters -------------------------------

filter_lengths <- c(11, 51, 201)
h <- lapply(filter_lengths, function(M) rep(1 / M, M))

### 3. Apply Filters -----------------------------------------------

y <- lapply(h, function(hf) stats::filter(x, hf, sides = 2))


### 4. Plot Results -------------------------------------------------

par(mfrow = c(4, 1), mar = c(5, 5, 4, 2))

## (0) Input signal
plot(time, x, type = "l", col = gray(0.7),
     main = "Input Signal: Gaussian Peaks + Noise",
     xlab = "Time", ylab = "Amplitude",
     xaxs = "i", ylim = c(-0.2, 1.5))
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Observed", "Clean"),
       col = c(gray(0.7), "red"), lty = 1, bty = "n")

## (1)-(3) Filtered signals
for (i in seq_along(filter_lengths)) {

  plot(time, y[[i]], type = "l", col = "blue",
       main = paste0("Moving Average Filter (length: ", filter_lengths[i], ")"),
       xlab = "Time", ylab = "Amplitude",
       xaxs = "i", ylim = c(-0.2, 1.5))

  # Overlay true clean signal (light red)
  lines(time, x.clean, col = rgb(1, 0.2, 0.3, alpha = 0.3), lwd = 2)

  legend("topright",
         legend = c("Filtered", "Clean"),
         col = c("blue", rgb(1, 0.2, 0.3, alpha = 0.3)),
         lty = 1, bty = "n")
}

Rスクリプト2：Windowed-sinc FIR ローパスフィルタ

########################################
# Comparison of Hamming-window FIR LPFs
# Varying cutoff frequency (fc)
# fc = 0.01, 0.03, 0.10
# Filter length is fixed (M = 101)
########################################

### 0. Utility: sinc function ---------------------------------------

sinc <- function(x) {
  y <- numeric(length(x))
  y[x == 0] <- 1
  y[x != 0] <- sin(pi * x[x != 0]) / (pi * x[x != 0])
  return(y)
}

### 1. Generate Input Signal ---------------------------------------

n.peaks <- 10
sig     <- exp(seq(log(200), log(2), length = n.peaks))
mu      <- max(sig) * 10 * (1:n.peaks)
dt      <- 1

t.range <- seq(0, max(mu) + max(sig) * 5, by = dt)
n.t     <- length(t.range)

# Clean signal: sum of Gaussian peaks
x.clean <- numeric(n.t)
for (i in 1:n.peaks) {
  x.clean <- x.clean +
    dnorm(t.range, mu[i], sig[i]) /
    dnorm(mu[i],    mu[i], sig[i])
}

# Observed signal with noise
x    <- x.clean + rnorm(n.t, sd = 0.05)
time <- seq(0, (n.t - 1) * dt, by = dt)

### 2. FIR design (Hamming-window LPF) ------------------------------

design_lpf_hamming <- function(M, fc) {
  # M : filter length (odd)
  # fc: cutoff frequency (cycles/sample)
  n   <- 0:(M - 1)
  mid <- (M - 1) / 2

  # Ideal LPF (sinc)
  h_ideal <- 2 * fc * sinc(2 * fc * (n - mid))

  # Hamming window
  w <- 0.54 - 0.46 * cos(2 * pi * n / (M - 1))

  # Windowed FIR
  h <- h_ideal * w

  # Normalize DC gain = 1
  h <- h / sum(h)
  return(h)
}

# Fixed filter length
M <- 101

# Three cutoff frequencies
fc.values <- c(0.10, 0.03, 0.01)

# Create 3 filters
h <- lapply(fc.values, function(fc) design_lpf_hamming(M, fc))
names(h) <- paste0("fc=", fc.values)

### 3. Filtering ----------------------------------------------------

y <- lapply(h, function(hf) stats::filter(x, hf, sides = 2))

### 4. Plot results -------------------------------------------------

par(mfrow = c(4, 1), mar = c(5, 4, 4, 2))

## (0) Input signal
plot(time, x, type = "l", col = gray(0.7),
     main = "Input Signal: Gaussian Peaks + Noise",
     xlab = "Time", ylab = "Amplitude",
     xaxs = "i", ylim = c(-0.2, 1.5))
lines(time, x.clean, col = "red", lwd = 2)
legend("topright",
       legend = c("Observed", "Clean"),
       col    = c(gray(0.7), "red"),
       lty    = 1, bty = "n")

## (1)-(3) Filters with different fc
for (i in seq_along(fc.values)) {

  plot(time, y[[i]], type = "l", col = "blue",
       main = paste0("Hamming-window FIR LPF (fc: ",
                     fc.values[i], ", length: ", M, ")"),
       xlab = "Time", ylab = "Amplitude",
       xaxs = "i", ylim = c(-0.2, 1.5))

  lines(time, x.clean, col = rgb(1, 0.2, 0.3, alpha = 0.3), lwd = 2)

  legend("topright",
         legend = c("Filtered", "Clean"),
         col    = c("blue", rgb(1, 0.2, 0.3, alpha = 0.3)),
         lty    = 1, bty = "n")
}

ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【デジタルフィルタの基礎１】まずは、FIRフィルタの世界を一望してみよう

1. なぜ「有限である」ことが重要なの？

① 絶対に安定する

② 設計しやすく予測しやすい

③ 完全な線形位相が実現できる

④ 安全で扱いやすい

2. FIR フィルタはどのように働くの？

■ 一般の（非因果を含む）FIR フィルタの定義

■ 因果型FIR フィルタ

3. FIR フィルタでできること

4. FIR フィルタの世界をどう理解するか

■ 位相特性から見る FIR フィルタの本質

■ 周波数特性から見るフィルタの機能

■ 設計手法：どのように FIR を作るか

（1）窓関数法（Window Method）— 最も簡単で実用的

（2）等リプル設計（Parks–McClellan 法）— 最適なフィルタ長

（3）周波数サンプリング法 — 周波数軸から直接設計

（4）最小位相フィルタ — 遅延を最小にしたいとき

（5）最適化ベース設計

■ 構造：どのように FIR を実装するか

（1）直接形（Direct Form）

（2）線形位相構造（対称性を利用）

（3）多相フィルタ（Polyphase Structure）

（4）FFT畳み込み（Overlap-Add/Save）

5. まとめ：深く理解する人だけが見える景色がある

[付録] Rでやってみる

Rスクリプト1：移動平均フィルタ

Rスクリプト2：Windowed-sinc FIR ローパスフィルタ