FIR（Finite Impulse Response）フィルタは、有限長のインパルス応答を持つディジタルフィルタです。IIRフィルタと異なり常に安定であり、設計次第で正確な線形位相特性が得られるため、音声処理・通信・計測など幅広い分野で利用されています。
FIRフィルタの設計法にはいくつかのアプローチがありますが、最も直感的で広く使われるのが 窓関数法 です。窓関数法では、理想フィルタのインパルス応答を有限長に打ち切り、その際に生じるスペクトル漏れ（ギブス現象）を窓関数によって制御します。
本記事の内容
- 理想ローパスフィルタのインパルス応答（sinc関数）の導出
- 窓関数法の原理（時間領域での打ち切りと周波数領域での畳み込み）
- ギブス現象の数学的説明
- 各種窓関数（矩形・ハニング・ハミング・ブラックマン・カイザー）の定義と特性
- カイザー窓のパラメータ設計公式
- 線形位相条件と対称性
- Pythonによる設計・比較・フィルタリング実験
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
理想ローパスフィルタのインパルス応答
理想ローパスフィルタの周波数応答
理想ローパスフィルタ（LPF）は、カットオフ周波数 $\omega_c$ 以下の周波数成分をすべて通過させ、それ以上の成分を完全に遮断するフィルタです。その周波数応答 $H_d(e^{j\omega})$ は次のように定義されます。
$$ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega| \le \omega_c \\ 0, & \omega_c < |\omega| \le \pi \end{cases} $$
ここで $\omega$ は正規化角周波数（$-\pi \le \omega \le \pi$）、$\omega_c$ はカットオフ角周波数です。
sinc関数の導出
この理想フィルタのインパルス応答 $h_d[n]$ を逆離散時間フーリエ変換（IDTFT）で求めます。
$$ h_d[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} H_d(e^{j\omega}) e^{j\omega n} d\omega $$
$H_d(e^{j\omega}) = 1$ となるのは $|\omega| \le \omega_c$ の範囲だけなので、
$$ h_d[n] = \frac{1}{2\pi} \int_{-\omega_c}^{\omega_c} 1 \cdot e^{j\omega n} d\omega $$
この積分を実行します。
$$ \begin{align} h_d[n] &= \frac{1}{2\pi} \left[ \frac{e^{j\omega n}}{jn} \right]_{-\omega_c}^{\omega_c} \\ &= \frac{1}{2\pi} \cdot \frac{e^{j\omega_c n} – e^{-j\omega_c n}}{jn} \\ &= \frac{1}{2\pi} \cdot \frac{2j\sin(\omega_c n)}{jn} \quad (\because \text{オイラーの公式}) \\ &= \frac{\sin(\omega_c n)}{\pi n} \end{align} $$
$n = 0$ の場合はロピタルの定理より、
$$ h_d[0] = \lim_{n \to 0} \frac{\sin(\omega_c n)}{\pi n} = \frac{\omega_c}{\pi} $$
まとめると、理想LPFのインパルス応答は次のようになります。
$$ \begin{equation} h_d[n] = \frac{\sin(\omega_c n)}{\pi n} = \frac{\omega_c}{\pi} \mathrm{sinc}\left(\frac{\omega_c n}{\pi}\right) \end{equation} $$
ここで $\mathrm{sinc}(x) = \sin(\pi x)/(\pi x)$ です。
この結果は重要です。理想LPFのインパルス応答は sinc関数 であり、$n \to \pm\infty$ でゆっくりと減衰しながら永遠に続きます。つまり、理想LPFは 因果的（causal）でなく、無限長のインパルス応答 を持つため、そのままでは実現できません。
窓関数法の原理
時間領域での打ち切り
理想フィルタのインパルス応答 $h_d[n]$ は無限長なので、これを有限長 $M$ に打ち切ります。最も単純な方法は、矩形窓 $w[n]$ を掛けることです。
$$ h[n] = h_d[n] \cdot w[n] $$
矩形窓は次のように定義されます。
$$ w_{\text{rect}}[n] = \begin{cases} 1, & 0 \le n \le M \\ 0, & \text{otherwise} \end{cases} $$
因果的なフィルタにするため、$h_d[n]$ を $M/2$ サンプルだけ遅延させます。すなわち、
$$ h[n] = h_d[n – M/2] \cdot w[n] = \frac{\sin(\omega_c (n – M/2))}{\pi (n – M/2)} \cdot w[n] $$
フィルタ長は $N = M + 1$ タップです。
周波数領域での畳み込み
時間領域での乗算は周波数領域での畳み込みに対応します。窓関数のフーリエ変換を $W(e^{j\omega})$ とすると、
$$ \begin{equation} H(e^{j\omega}) = \frac{1}{2\pi} H_d(e^{j\omega}) * W(e^{j\omega}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} H_d(e^{j\theta}) W(e^{j(\omega – \theta)}) d\theta \end{equation} $$
この畳み込みにより、理想フィルタの急峻な遷移帯域が 滑らかに なります。窓関数のスペクトル $W(e^{j\omega})$ の形状が、設計されたフィルタの周波数応答の品質を決定します。
ギブス現象
矩形窓で単純に打ち切ると、理想フィルタの不連続点（$\omega = \omega_c$）の近傍で振動（リンギング）が発生します。これが ギブス現象 です。
矩形窓のフーリエ変換は、長さ $N$ のディリクレ核で、
$$ W_{\text{rect}}(e^{j\omega}) = e^{-j\omega M/2} \cdot \frac{\sin(\omega N/2)}{\sin(\omega/2)} $$
この関数は大きなサイドローブを持ちます。畳み込みの際、このサイドローブが理想フィルタの不連続点にぶつかり、約 8.95%（約 -21 dB）のオーバーシュートを生じます。
重要なのは、フィルタ長 $N$ を増やしてもオーバーシュートの 振幅比率 は変わらないという点です。$N$ を大きくするとリンギングの幅は狭くなりますが、ピーク値は 8.95% のまま残ります。これがギブス現象の本質です。
各種窓関数の定義と特性
ギブス現象を抑制するために、矩形窓以外の様々な窓関数が考案されています。窓関数の設計では、主ローブ幅（遷移帯域幅に影響）と サイドローブ減衰量（阻止域の減衰に影響）のトレードオフがあります。
以下、すべて長さ $N = M + 1$ の窓関数 $w[n]$（$0 \le n \le M$）として定義します。
矩形窓（Rectangular Window）
$$ w[n] = 1 $$
- 主ローブ幅: $4\pi / N$
- 第一サイドローブ減衰: $-13$ dB
- 最小阻止域減衰: $-21$ dB
最も狭い主ローブ幅を持ちますが、サイドローブが大きく、実用的なフィルタ設計にはほとんど使いません。
ハニング窓（Hanning / von Hann Window）
$$ w[n] = 0.5 – 0.5\cos\left(\frac{2\pi n}{M}\right) $$
- 主ローブ幅: $8\pi / N$
- 第一サイドローブ減衰: $-31$ dB
- 最小阻止域減衰: $-44$ dB
余弦の持ち上げ（raised cosine）で両端が滑らかにゼロに近づきます。サイドローブ減衰率は $1/\omega^3$ で、矩形窓の $1/\omega$ より速く減衰します。
ハミング窓（Hamming Window）
$$ w[n] = 0.54 – 0.46\cos\left(\frac{2\pi n}{M}\right) $$
- 主ローブ幅: $8\pi / N$
- 第一サイドローブ減衰: $-43$ dB
- 最小阻止域減衰: $-53$ dB
ハニング窓の係数を調整して、第一サイドローブを最小化したものです。係数 $0.54$ と $0.46$ は、ディリクレ核のサイドローブのピークと余弦関数のサイドローブのピークがちょうど打ち消し合うように選ばれています。
具体的には、ハミング窓のスペクトルは
$$ W_{\text{Ham}}(e^{j\omega}) = 0.54 \cdot W_{\text{rect}}(e^{j\omega}) + 0.23 \cdot W_{\text{rect}}(e^{j(\omega – 2\pi/M)}) + 0.23 \cdot W_{\text{rect}}(e^{j(\omega + 2\pi/M)}) $$
と分解でき、矩形窓のサイドローブをシフトした余弦成分で打ち消しています。
ブラックマン窓（Blackman Window）
$$ w[n] = 0.42 – 0.5\cos\left(\frac{2\pi n}{M}\right) + 0.08\cos\left(\frac{4\pi n}{M}\right) $$
- 主ローブ幅: $12\pi / N$
- 第一サイドローブ減衰: $-58$ dB
- 最小阻止域減衰: $-74$ dB
3項の余弦和で構成され、非常に小さなサイドローブを実現します。ただし主ローブ幅が広がるため、遷移帯域幅は大きくなります。
カイザー窓（Kaiser Window）
$$ \begin{equation} w[n] = \frac{I_0\left(\beta \sqrt{1 – \left(\frac{2n}{M} – 1\right)^2}\right)}{I_0(\beta)} \end{equation} $$
ここで $I_0(\cdot)$ は第0種変形ベッセル関数、$\beta$ は形状パラメータです。
カイザー窓は $\beta$ を連続的に変化させることで、主ローブ幅とサイドローブ減衰量のトレードオフを自在に制御できる、非常に柔軟な窓関数です。
- $\beta = 0$: 矩形窓に一致
- $\beta \approx 5.44$: ハミング窓に近い特性
- $\beta \approx 8.96$: ブラックマン窓に近い特性
カイザー窓のパラメータ設計公式
カイザーは、所望の阻止域減衰量 $A_s$ [dB] と遷移帯域幅 $\Delta\omega$ [rad/sample] から、$\beta$ とフィルタ次数 $M$ を求める経験的公式を与えています。
$\beta$ の決定
阻止域減衰量 $A_s = -20\log_{10}(\delta_s)$ dB（$\delta_s$ はリプル振幅）に対して、
$$ \begin{equation} \beta = \begin{cases} 0.1102(A_s – 8.7), & A_s > 50 \\ 0.5842(A_s – 21)^{0.4} + 0.07886(A_s – 21), & 21 \le A_s \le 50 \\ 0, & A_s < 21 \end{cases} \end{equation} $$
フィルタ次数 $M$ の決定
$$ \begin{equation} M = \frac{A_s – 7.95}{2.285 \cdot \Delta\omega} \end{equation} $$
ここで $\Delta\omega$ はラジアン/サンプルの遷移帯域幅です。
この公式により、設計仕様（通過域リプル、阻止域減衰、遷移帯域幅）から直接フィルタパラメータを決定できます。
線形位相条件と対称性
FIRフィルタの重要な利点は 正確な線形位相 を実現できることです。線形位相とは、位相応答が周波数に対して直線的になることを意味し、すべての周波数成分が同じ時間遅延を受けます。
線形位相の条件
FIRフィルタ $h[n]$（$n = 0, 1, \ldots, M$）が線形位相を持つための十分条件は、インパルス応答が 対称 または 反対称 であることです。
Type I（偶数対称、奇数長）:
$$ h[n] = h[M – n], \quad M: \text{偶数} $$
Type II（偶数対称、偶数長）:
$$ h[n] = h[M – n], \quad M: \text{奇数} $$
Type III（奇数対称、奇数長）:
$$ h[n] = -h[M – n], \quad M: \text{偶数} $$
Type IV（奇数対称、偶数長）:
$$ h[n] = -h[M – n], \quad M: \text{奇数} $$
群遅延
線形位相FIRフィルタの位相応答は $\angle H(e^{j\omega}) = -\omega M/2 + \phi_0$ の形をとり、群遅延は
$$ \tau_g = -\frac{d\angle H(e^{j\omega})}{d\omega} = \frac{M}{2} \quad \text{(一定)} $$
となります。群遅延が一定ということは、すべての周波数成分が $M/2$ サンプルだけ遅延して出力されることを意味し、波形の歪みが発生しません。
窓関数法での対称性の保証
窓関数法で設計する場合、sinc関数（$h_d[n – M/2]$）は $n = M/2$ に関して自然に対称です。一般的な窓関数 $w[n]$ もまた $n = M/2$ に関して対称です。対称な関数同士の積は対称なので、窓関数法で得られるFIRフィルタは自動的に Type I（$M$ 偶数）または Type II（$M$ 奇数）の線形位相条件を満たします。
Pythonでの実装
各種窓関数の比較
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# パラメータ
N = 51 # フィルタ長（タップ数）
M = N - 1 # フィルタ次数
# 各種窓関数の生成
n = np.arange(N)
windows = {
'Rectangular': np.ones(N),
'Hanning': 0.5 - 0.5 * np.cos(2 * np.pi * n / M),
'Hamming': 0.54 - 0.46 * np.cos(2 * np.pi * n / M),
'Blackman': 0.42 - 0.5 * np.cos(2 * np.pi * n / M) + 0.08 * np.cos(4 * np.pi * n / M),
'Kaiser (β=5.0)': np.kaiser(N, beta=5.0),
'Kaiser (β=8.6)': np.kaiser(N, beta=8.6),
}
# --- 窓関数の波形 ---
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
for name, w in windows.items():
axes[0].plot(n, w, label=name)
axes[0].set_xlabel('Sample n')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Window Functions (N=51)')
axes[0].legend()
axes[0].grid(True)
# --- 窓関数の周波数応答（対数スケール） ---
for name, w in windows.items():
# ゼロパディングしてFFT
nfft = 4096
W = np.fft.fft(w, nfft)
W_shift = np.fft.fftshift(W)
freq = np.linspace(-np.pi, np.pi, nfft)
W_dB = 20 * np.log10(np.abs(W_shift) / np.max(np.abs(W_shift)) + 1e-12)
axes[1].plot(freq / np.pi, W_dB, label=name)
axes[1].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[1].set_ylabel('Magnitude (dB)')
axes[1].set_title('Window Spectra')
axes[1].set_ylim([-120, 5])
axes[1].legend(fontsize=8)
axes[1].grid(True)
plt.tight_layout()
plt.show()
このコードでは、6種類の窓関数の時間波形と周波数スペクトルを比較します。矩形窓は主ローブが最も狭い一方でサイドローブが大きく、ブラックマン窓やカイザー窓（$\beta$ 大）はサイドローブが小さい代わりに主ローブが広いことが確認できます。
FIRローパスフィルタの設計
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
def design_fir_lpf_window(N, wc, window_type='hamming', beta=5.0):
"""窓関数法によるFIR LPF設計
Parameters
----------
N : int
フィルタ長（タップ数、奇数推奨）
wc : float
カットオフ角周波数 [rad/sample]
window_type : str
窓関数の種類
beta : float
カイザー窓のβパラメータ
Returns
-------
h : ndarray
フィルタ係数
"""
M = N - 1
n = np.arange(N)
# 理想LPFのインパルス応答（M/2だけ遅延）
alpha = M / 2.0
h_d = np.where(
n == alpha,
wc / np.pi, # n = alpha のとき（ロピタルの定理）
np.sin(wc * (n - alpha)) / (np.pi * (n - alpha))
)
# 窓関数の選択
if window_type == 'rectangular':
w = np.ones(N)
elif window_type == 'hanning':
w = 0.5 - 0.5 * np.cos(2 * np.pi * n / M)
elif window_type == 'hamming':
w = 0.54 - 0.46 * np.cos(2 * np.pi * n / M)
elif window_type == 'blackman':
w = 0.42 - 0.5 * np.cos(2 * np.pi * n / M) + 0.08 * np.cos(4 * np.pi * n / M)
elif window_type == 'kaiser':
w = np.kaiser(N, beta)
else:
raise ValueError(f'Unknown window type: {window_type}')
# 窓関数を適用
h = h_d * w
return h
# 設計パラメータ
N = 51 # フィルタ長
wc = 0.4 * np.pi # カットオフ周波数 (0.4π rad/sample)
# 各種窓関数で設計
window_types = ['rectangular', 'hanning', 'hamming', 'blackman', 'kaiser']
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
for wtype in window_types:
h = design_fir_lpf_window(N, wc, window_type=wtype, beta=5.44)
# 周波数応答
w_freq, H = signal.freqz(h, worN=4096)
H_dB = 20 * np.log10(np.abs(H) + 1e-12)
axes[0].plot(w_freq / np.pi, H_dB, label=wtype)
axes[1].plot(w_freq / np.pi, np.abs(H), label=wtype)
# 理想特性を表示
axes[0].axvline(x=0.4, color='k', linestyle='--', alpha=0.5, label='Cutoff')
axes[0].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title(f'FIR LPF Frequency Response (N={N}, ωc=0.4π)')
axes[0].set_ylim([-100, 5])
axes[0].legend(fontsize=8)
axes[0].grid(True)
axes[1].axvline(x=0.4, color='k', linestyle='--', alpha=0.5, label='Cutoff')
axes[1].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[1].set_ylabel('Magnitude')
axes[1].set_title('FIR LPF Frequency Response (Linear Scale)')
axes[1].legend(fontsize=8)
axes[1].grid(True)
plt.tight_layout()
plt.show()
FIRバンドパスフィルタの設計
バンドパスフィルタ（BPF）は、理想BPFのインパルス応答から設計します。理想BPFの周波数応答は
$$ H_{\text{BPF}}(e^{j\omega}) = \begin{cases} 1, & \omega_{c1} \le |\omega| \le \omega_{c2} \\ 0, & \text{otherwise} \end{cases} $$
であり、対応するインパルス応答は
$$ h_{\text{BPF}}[n] = \frac{\sin(\omega_{c2}(n – M/2))}{\pi(n – M/2)} – \frac{\sin(\omega_{c1}(n – M/2))}{\pi(n – M/2)} $$
です。これは「上側カットオフのLPF」と「下側カットオフのLPF」の差として得られます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
def design_fir_bpf_window(N, wc1, wc2, window_type='hamming', beta=5.0):
"""窓関数法によるFIR BPF設計"""
M = N - 1
n = np.arange(N)
alpha = M / 2.0
# 理想BPFのインパルス応答 = LPF(wc2) - LPF(wc1)
h_d = np.where(
n == alpha,
(wc2 - wc1) / np.pi,
np.sin(wc2 * (n - alpha)) / (np.pi * (n - alpha))
- np.sin(wc1 * (n - alpha)) / (np.pi * (n - alpha))
)
# 窓関数を適用
if window_type == 'hamming':
w = 0.54 - 0.46 * np.cos(2 * np.pi * n / M)
elif window_type == 'kaiser':
w = np.kaiser(N, beta)
else:
w = np.ones(N)
return h_d * w
# BPF設計: 通過域 0.3π ~ 0.6π
N = 71
wc1 = 0.3 * np.pi
wc2 = 0.6 * np.pi
h_bpf = design_fir_bpf_window(N, wc1, wc2, window_type='hamming')
# 周波数応答の計算
w_freq, H_bpf = signal.freqz(h_bpf, worN=4096)
H_bpf_dB = 20 * np.log10(np.abs(H_bpf) + 1e-12)
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# インパルス応答
axes[0].stem(np.arange(N), h_bpf, markerfmt='C0o', basefmt='k-')
axes[0].set_xlabel('Sample n')
axes[0].set_ylabel('h[n]')
axes[0].set_title(f'FIR BPF Impulse Response (N={N})')
axes[0].grid(True)
# 周波数応答
axes[1].plot(w_freq / np.pi, H_bpf_dB)
axes[1].axvline(x=0.3, color='r', linestyle='--', alpha=0.7, label='ωc1=0.3π')
axes[1].axvline(x=0.6, color='r', linestyle='--', alpha=0.7, label='ωc2=0.6π')
axes[1].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[1].set_ylabel('Magnitude (dB)')
axes[1].set_title('FIR BPF Frequency Response')
axes[1].set_ylim([-80, 5])
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
フィルタリング実験
実際の信号にFIRフィルタを適用してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# テスト信号: 3つの周波数成分（0.1π, 0.4π, 0.8π）
np.random.seed(42)
n_samples = 300
n = np.arange(n_samples)
fs = 1.0 # 正規化サンプリング周波数
# 信号成分
x1 = np.sin(0.1 * np.pi * n) # 低周波 (0.1π)
x2 = 0.8 * np.sin(0.4 * np.pi * n) # 中間周波 (0.4π)
x3 = 0.5 * np.sin(0.8 * np.pi * n) # 高周波 (0.8π)
noise = 0.3 * np.random.randn(n_samples) # ガウスノイズ
x = x1 + x2 + x3 + noise # 合成信号
# FIR LPFの設計 (ωc = 0.25π, ハミング窓, N=61)
N_filt = 61
wc = 0.25 * np.pi
M = N_filt - 1
nn = np.arange(N_filt)
alpha = M / 2.0
h_d = np.where(
nn == alpha,
wc / np.pi,
np.sin(wc * (nn - alpha)) / (np.pi * (nn - alpha))
)
w_ham = 0.54 - 0.46 * np.cos(2 * np.pi * nn / M)
h_lpf = h_d * w_ham
# フィルタリング
y = signal.lfilter(h_lpf, 1.0, x)
# 群遅延分の補正
delay = M // 2
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
# 入力信号
axes[0].plot(n, x, 'b', alpha=0.7, label='Input x[n]')
axes[0].set_xlabel('Sample')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Input Signal (0.1π + 0.4π + 0.8π + Noise)')
axes[0].legend()
axes[0].grid(True)
axes[0].set_xlim([0, 200])
# フィルタリング結果（群遅延補正済み）
axes[1].plot(n, x, 'b', alpha=0.3, label='Input')
axes[1].plot(n[delay:], y[delay:n_samples], 'r', linewidth=2, label='Filtered (delay compensated)')
axes[1].plot(n, x1, 'g--', alpha=0.7, label='Original 0.1π component')
axes[1].set_xlabel('Sample')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('FIR LPF Output (ωc=0.25π, Hamming, N=61)')
axes[1].legend()
axes[1].grid(True)
axes[1].set_xlim([0, 200])
# フィルタの周波数応答
w_freq, H_filt = signal.freqz(h_lpf, worN=4096)
axes[2].plot(w_freq / np.pi, 20 * np.log10(np.abs(H_filt) + 1e-12))
axes[2].axvline(x=0.1, color='g', linestyle='--', alpha=0.7, label='0.1π')
axes[2].axvline(x=0.4, color='orange', linestyle='--', alpha=0.7, label='0.4π')
axes[2].axvline(x=0.8, color='r', linestyle='--', alpha=0.7, label='0.8π')
axes[2].axvline(x=0.25, color='k', linestyle='-', alpha=0.5, label='Cutoff 0.25π')
axes[2].set_xlabel('Normalized Frequency (×π rad/sample)')
axes[2].set_ylabel('Magnitude (dB)')
axes[2].set_title('FIR LPF Frequency Response')
axes[2].set_ylim([-80, 5])
axes[2].legend()
axes[2].grid(True)
plt.tight_layout()
plt.show()
フィルタリングの結果、0.1π の低周波成分のみが通過し、0.4π と 0.8π の成分およびノイズが除去されていることが確認できます。また、FIRフィルタの線形位相特性により、出力波形は群遅延 $M/2 = 30$ サンプル分だけ遅延するものの、波形の形状は保存されます。
カイザー窓による仕様駆動設計
最後に、カイザー窓の設計公式を使って、仕様（阻止域減衰量・遷移帯域幅）から自動的にフィルタを設計する例を示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.special import i0 as bessel_i0
def kaiser_fir_design(wp, ws, As, fs=1.0):
"""カイザー窓によるFIR LPF設計
Parameters
----------
wp : float
通過域端周波数 [Hz]
ws : float
阻止域端周波数 [Hz]
As : float
阻止域減衰量 [dB]
fs : float
サンプリング周波数 [Hz]
Returns
-------
h : ndarray
フィルタ係数
N : int
フィルタ長
beta : float
カイザー窓パラメータ
"""
# 正規化角周波数に変換
wp_rad = 2 * np.pi * wp / fs
ws_rad = 2 * np.pi * ws / fs
# 遷移帯域幅
delta_w = ws_rad - wp_rad
# カットオフ周波数（通過域と阻止域の中間）
wc = (wp_rad + ws_rad) / 2
# β の決定
if As > 50:
beta = 0.1102 * (As - 8.7)
elif As >= 21:
beta = 0.5842 * (As - 21)**0.4 + 0.07886 * (As - 21)
else:
beta = 0.0
# フィルタ次数 M の決定
M = int(np.ceil((As - 7.95) / (2.285 * delta_w)))
if M % 2 == 1: # 奇数タップにするため偶数次数に
M += 1
N = M + 1
# フィルタ設計
n = np.arange(N)
alpha = M / 2.0
h_d = np.where(
n == alpha,
wc / np.pi,
np.sin(wc * (n - alpha)) / (np.pi * (n - alpha))
)
w = np.kaiser(N, beta)
h = h_d * w
return h, N, beta
# 設計仕様
# 通過域: 0 ~ 1000 Hz, 阻止域: 1500 Hz ~, 減衰: 60 dB
# サンプリング周波数: 8000 Hz
fs = 8000
wp = 1000
ws = 1500
As = 60
h, N_filt, beta = kaiser_fir_design(wp, ws, As, fs)
print(f'フィルタ長: {N_filt}, β: {beta:.4f}')
# 周波数応答
w_freq, H = signal.freqz(h, worN=4096, fs=fs)
H_dB = 20 * np.log10(np.abs(H) + 1e-12)
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
axes[0].plot(w_freq, H_dB)
axes[0].axvline(x=wp, color='g', linestyle='--', label=f'Pass edge: {wp} Hz')
axes[0].axvline(x=ws, color='r', linestyle='--', label=f'Stop edge: {ws} Hz')
axes[0].axhline(y=-As, color='orange', linestyle=':', label=f'-{As} dB')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title(f'Kaiser Window FIR LPF (N={N_filt}, β={beta:.2f})')
axes[0].set_ylim([-100, 5])
axes[0].legend()
axes[0].grid(True)
# 位相応答
phase = np.unwrap(np.angle(H))
axes[1].plot(w_freq, phase)
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Phase (rad)')
axes[1].set_title('Phase Response (Linear Phase)')
axes[1].grid(True)
plt.tight_layout()
plt.show()
カイザー窓の設計公式を使うことで、仕様を満たすフィルタが自動的に設計されます。位相応答が直線的であること（線形位相）も確認できます。
まとめ
本記事では、FIRフィルタの窓関数法による設計について解説しました。
- 理想LPFのインパルス応答は sinc 関数であり、逆DTFTから導出される
- 窓関数法は理想フィルタのインパルス応答を有限長に打ち切る手法であり、時間領域の乗算は周波数領域の畳み込みに対応する
- 矩形窓での打ち切りではギブス現象が発生し、サイドローブ（約 -21 dB）が問題になる
- ハニング窓・ハミング窓・ブラックマン窓は余弦和の係数調整でサイドローブを抑制する
- カイザー窓はパラメータ $\beta$ により主ローブ幅とサイドローブ減衰のトレードオフを連続的に制御でき、設計公式により仕様から直接パラメータを決定できる
- FIRフィルタはインパルス応答の対称性により正確な線形位相を保証し、群遅延が一定になる
次のステップとして、以下の記事も参考にしてください。