デジタル信号処理9
note Item Type Metadata
note
## 1. ディジタル信号処理の背景
- 音声などのアナログ信号は一定周期 $T$ でサンプリング・量子化され、列 $x(n)$ になる。サンプル周波数は
$$
F_s = \frac{1}{T}
$$
で定義される。
- 例: 従来電話では
$$
T = \frac{1}{8000}\,\text{s}
$$
音楽 CD では
$$
T = \frac{1}{44100}\,\text{s}
$$
- ナイキスト標本化定理:
$$
F_s \ge 2F_{\max}
$$
を満たせば帯域制限信号は損失なく復元可能。$T$ を小さくして $F_s$ を上げると可聴帯域全体（約 20 kHz）をカバーできる。
- 量子化は振幅を有限ビットで表現する過程で、量子化雑音（白色雑音近似）が付加される。ビット数 $B$ を増やすと雑音電力は
$$
\sigma_q^2 \approx \frac{1}{12}\Delta^2
$$
（$\Delta$: 量子化ステップ）に低減する。
- DSP ではサンプル列に離散畳み込みや周波数変換を施し、改善後に D/A 変換・ポストフィルタを経てアナログ領域へ戻す。計算は固定小数点/浮動小数点プロセッサや FPGA で行う。
---
## 2. システム構成におけるフィルタ
1. アナログ入力 → プリフィルタ（アンチエイリアシング / ブリッジングネットワーク）
2. サンプリング & A/D 変換 → デジタル信号 `x(n)`（必要に応じてデジタルダウンコンバータで IF→ベースバンド変換）
3. DSP 上のディジタルフィルタ → 出力 `y(n)`（FIR/IIR、複素処理を含む）
4. D/A 変換 → ポストフィルタ（スペクトルイメージ除去） → アナログ出力
- 処理フローを図示すると以下の通り。
```mermaid
flowchart LR
A[アナログ入力] --> B[プリフィルタ<br>アンチエイリアシング]
B --> C[A/D変換<br>サンプリング・量子化]
C --> D[DSP処理<br>ディジタルフィルタ]
D --> E[D/A変換]
E --> F[ポストフィルタ<br>イメージ抑制]
F --> G[アナログ出力]
```
- プリフィルタは
$$
\frac{F_s}{2}
$$
以上の高域成分を除去し、折り返し雑音を防ぐ。多くはアクティブローパス+アナログアンチエイリアスで構成。
- ポストフィルタは D/A 後に生じるイメージスペクトル（`Fs` の高調波）を抑え、再び連続信号として滑らかにする。
- DSP 内部では加算・乗算・遅延素子で任意のインパルス応答を合成する。FIR は有限インパルス応答で厳密線形位相を実現しやすく、IIR はより鋭い遷移帯域を少次数で実現できる。
- **理解度チェック:** プリフィルタとポストフィルタの役割・設計制約（遷移帯域、群遅延、実装コスト）の違いを説明できるか?
---
- 定義: 直近 $M$ サンプルの平均を出力する等重み FIR ローパス。畳み込み核は
$$
h(n) = \frac{1}{M}, \quad 0 \le n \le M-1
$$
で、それ以外は $0$。
- 数式:
$$
y(n) = \frac{1}{M} \sum_{k=0}^{M-1} x(n-k)
$$
未定義の過去サンプルは $0$ と置く（ゼロ初期条件）。$z$ 領域伝達関数は
$$
H(z) = \frac{1}{M}\left(1 + z^{-1} + \dots + z^{-(M-1)}\right)
= \frac{1}{M}\frac{1 - z^{-M}}{1 - z^{-1}}
$$
であり、単位円上には $M-1$ 個の等間隔ゼロを持つ。
- 特徴:
- 実装容易（加算とスカラー倍のみ）。固定小数点実装ではシフト演算で $1/M$ を近似できる場合も多い（$M = 2^k$ など）。
- 位相は線形で群遅延
$$
\tau_g = \frac{M-1}{2}
$$
サンプルを生むため、リアルタイム制御では遅延許容度の検討が必須。
- 周波数応答はディリクレ核に相当し、主ローブ幅
$$
\approx \frac{4\pi}{M}
$$
サイドローブ減衰 $\approx 13 \text{ dB}$ 程度。滑らかなローパスだがストップバンド減衰は限定的。
- **理解度チェック:** `M` を 2 倍にした際、群遅延と主ローブ幅がどう変化するか説明できるか?
---
## 4. 3点平均フィルタの逐次計算例
入力 $x(n)$:
| n | 0 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| x(n) | 1 | -1 | 3 | 4 | 2 |
出力（$M = 3$）:
| n | 計算式 | y(n) |
|---|---|---|
| 0 | $$\frac{1}{3}(x_0 + x_{-1} + x_{-2})$$ | $$\frac{1}{3}$$ |
| 1 | $$\frac{1}{3}(x_1 + x_0 + x_{-1})$$ | $$0$$ |
| 2 | $$\frac{1}{3}(x_2 + x_1 + x_0)$$ | $$1$$ |
| 3 | $$\frac{1}{3}(x_3 + x_2 + x_1)$$ | $$2$$ |
| 4 | $$\frac{1}{3}(x_4 + x_3 + x_2)$$ | $$3$$ |
逐次実装では、「新サンプルを加算」「最古サンプルを減算」「$1/M$ を乗算」の 3 ステップで $O(1)$ 更新できる。累積和
$$
S(n) = x(n) + \dots + x(n-M+1)
$$
を保持すれば
$$
S(n) = S(n-1) + x(n) - x(n-M)
$$
で更新でき、乗算器が制約となる場合は $M = 2^k$ でシフト実装する。FPGA ではシフトレジスタと加算器ツリーで 1 サイクル遅延のパイプラインを構成できる。
- **実験課題:** Python/NumPy や MATLAB でホワイトノイズに 3 点移動平均を適用し、時間波形・周波数応答を確認せよ。
---
## 5. 雑音抑制の理由
- **時間領域の直感**: 雑音は急峻・ランダムな変動を含むため平均化で相殺され、滑らかになる。独立同分布雑音 $w(n)$ に対して出力雑音分散は
$$
\frac{\sigma_w^2}{M}
$$
まで低減し、
$$
10 \log_{10} M\ \text{[dB]}
$$
の SNR 向上が得られる。
- **$M$ と性能のトレードオフ**: $M$ を大きくすると雑音抑制は強力だが、信号の立ち上がりや変調成分も平滑化される。エッジ再現が重要な画像では $M$ を小さく保ち、スペクトル的に限定されたセンサ雑音除去では大きな $M$ を許容するなど、用途依存の最適化が必要。
- **周波数領域の視点**: フィルタ係数 $[1/M, \dots, 1/M]$ のフーリエ変換は
$$
H(e^{j\omega}) = \frac{1}{M} \frac{\sin(M\omega/2)}{\sin(\omega/2)} e^{-j\omega(M-1)/2}
$$
となり、主ローブ内で通過、サイドローブで減衰するローパス特性を示す。入力スペクトル $X(e^{j\omega})$ とフィルタ $H(e^{j\omega})$ の積が出力スペクトル $Y(e^{j\omega})$ であり、高周波雑音が抑圧される。
- **PSD 解析**: 入力雑音のパワースペクトル密度 $S_w(\omega)$ に対して出力 PSD は
$$
\left|H(e^{j\omega})\right|^2 S_w(\omega)
$$
。白色雑音なら $|H|^2$ に比例して減衰する。低周波寄りのピンク雑音では完全には落ちないため、追加の重み付けや適応フィルタ（LMS、RLS）との組み合わせが検討される。
---
## 6. Python 実装例
以下は `M=7` の移動平均フィルタでホワイトノイズ混入信号を平滑化し、SNR 改善と周波数応答を確認する最小実験コード。
```python
import numpy as np
from numpy.fft import rfft, rfftfreq
import matplotlib.pyplot as plt
def moving_average(x: np.ndarray, M: int) -> np.ndarray:
kernel = np.ones(M) / M
return np.convolve(x, kernel, mode="same")
Fs = 48_000
duration = 0.05
t = np.arange(0, duration, 1 / Fs)
signal = np.sin(2 * np.pi * 1_000 * t)
noise = 0.4 * np.random.randn(t.size)
x = signal + noise
M = 7
y = moving_average(x, M)
snr_in = 10 * np.log10(np.mean(signal**2) / np.mean(noise**2))
snr_out = 10 * np.log10(np.mean(signal**2) / np.mean((y - signal)**2))
freq = rfftfreq(t.size, 1 / Fs)
X = np.abs(rfft(x))
Y = np.abs(rfft(y))
print(f"SNR in = {snr_in:.2f} dB")
print(f"SNR out = {snr_out:.2f} dB")
plt.figure(figsize=(10, 4))
plt.plot(freq, X, label="Input |X|")
plt.plot(freq, Y, label="Output |Y|")
plt.xlim(0, 10_000)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.legend()
plt.tight_layout()
plt.show()
```
$$
\text{SNR}_{\text{out}} - \text{SNR}_{\text{in}}
$$
が理論値
$$
10 \log_{10} M \approx 8.45\ \text{dB}
$$
に近いこと、スペクトルで高周波が抑圧されることを確認できる。$M$ や信号周波数、雑音統計を変えれば課題 1〜4 の検証を自動化できる。
同じ処理フローを図示すると以下の通り。プログラムに馴染みがなくても、データ生成から可視化までの段階が追える。
```mermaid
flowchart TD
A[パラメータ設定<br>Fs, M, duration] --> B[時間軸 t 生成]
B --> C[信号生成<br>sin波]
B --> D[雑音生成<br>正規乱数]
C --> E[信号+雑音で x を作成]
D --> E
E --> F[移動平均フィルタ<br>np.convolve]
F --> G[SNR計算<br>入力または出力]
E --> H[FFTで |X| を計算]
F --> I[FFTで |Y| を計算]
G --> J[数値表示]
H --> K[周波数軸 freq]
I --> K
K --> L[プロット表示]
```
flowchart TD
A[パラメータ設定<br>Fs, M, duration] --> B[時間軸 t 生成]
B --> C[信号生成<br>sin波]
B --> D[雑音生成<br>正規乱数]
C --> E[信号+雑音で x を作成]
D --> E
E --> F[移動平均フィルタ<br>np.convolve]
F --> G[SNR計算<br>入力または出力]
E --> H[FFTで |X| を計算]
F --> I[FFTで |Y| を計算]
G --> J[数値表示]
H --> K[周波数軸 freq]
I --> K
K --> L[プロット表示]
```
Collection
Citation
unjuno, “デジタル信号処理9,” unjuno'sResearchLibrary, accessed October 8, 2026, https://archive.unjuno.org/items/show/131.
コメント