デジタル信号処理5
note Item Type Metadata
note
フーリエ変換とスペクトル解析
この資料について
この資料では、フーリエ変換とスペクトル解析の基礎から応用までを体系的に学習します。理論的な理解と実用的な応用の両方を重視し、段階的に学習できるよう構成されています。
学習目標:
- フーリエ変換の基本原理を理解する
- DTFTとDFTの関係を把握する
- 音響信号解析への応用を学ぶ
- 実装時の考慮事項を理解する
想定学習時間: 約8-10時間
学習手順
アナログ(連続)信号処理
- フーリエ級数 - 周期関数の基本
- 複素フーリエ級数 - オイラーの公式を活用
- フーリエ変換 - 非周期関数への拡張
デジタル(離散)信号処理
- 離散時間フーリエ変換(DTFT) - 離散時間信号の周波数解析
- 離散フーリエ変換(DFT) - コンピュータ実装向け
フーリエ変換とは
フーリエ変換とは、複雑な波や音を「いくつかの単純な波(正弦波)」に分けて考える方法です。
身近な例で理解する
音楽の例:
好きな音楽の音も、実はいろんな高さの音が混ざっています。フーリエ変換を使うと、どんな音がどれくらい混ざっているのかわかりやすくできます。
光の例:
ちょうど、虹が白い光をいろんな色に分けるのと似ています。
基本的な概念
- 用途: 信号の周波数解析に用いられる
- 効果: 信号の波形の見方を変える
- 重要性: デジタル信号処理の基礎技術
スペクトラムアナライザー
スペクトラムアナライザーとは
音や信号に含まれるさまざまな周波数(音の高さ)をグラフで表示して分析できる装置やソフトウェアです。実用のスペクトラムアナライザーは多くの場合 FFT を用い、窓関数、スケール補正(RMS/dBFS/dBV/dBm)、平均化(線形平均/Welch)などを行います。表示は PSD(dB/Hz)またはパワースペクトル(dB)が一般的です。PSD表示は V²/Hz, dBV²/Hz, dBm/Hz など単位を明示(帯域幅依存の誤読防止)。
実装時の重要な注意点
振幅推定の正確な方法
片側スペクトルからの振幅推定:
- 片側表示:$A_{\text{est}} = \frac{|X_{w}^{(1\text{-side})}[k_0]|}{N \cdot \mathrm{CG}}$(内点×2が既に適用済み)
- 二側表示:$A_{\text{est}} = \frac{2|X_w^{(2\text{-side})}[k_0]|}{N \cdot \mathrm{CG}}$(DC/Nyquistは倍化なし)
注意:片側表示では内点を×2して表示するため、そのままCGで割れば振幅が得られます。
PSDスケーリングの正しい方法
Welch/Periodogram標準形:
これをENBW & CGで等価変形すると:
重要:CG²の補正を忘れずに適用してください。
Parsevalの定理とPSDの関係
理解のポイント:窓掛け後のスペクトル $X_w$ を用いると、PSD帯域和は窓エネルギー $U$ を含みます。「元信号のエネルギー」と「窓込みPSDの帯域和」は一般に一致しません。
実装時の注意:測定対象(電圧/電力・帯域)を明記し、CG/ENBW/Uを適切に使い分けてください。
窓補正と ENBW(有効雑音帯域幅)
窓を掛けたスペクトルを物理量に正しく換算するために、以下を明示する:
- コヒーレントゲイン(CG):\(\displaystyle \mathrm{CG} = \frac{1}{N}\sum_{n=0}^{N-1} w[n] \tag{3}\)。
ピーク(正弦1本)の振幅は 振幅補正として \(1/\mathrm{CG}\) を掛けて復元する。
片側スペクトルからピーク振幅を読む場合:\(A_{\text{est}} = \frac{|X_w^{(1\text{-side})}[k_0]|}{N \cdot \mathrm{CG}} \tag{4}\)(片側表示では内点×2が既に適用されているため、どのビンでも同じ式)。
二側スペクトルの場合:\(A_{\text{est}} = \frac{2|X_w^{(2\text{-side})}[k_0]|}{N \cdot \mathrm{CG}} \tag{5}\)(実正弦で対称2本が立つときの各ビン基準)。 - ENBW(有効雑音帯域幅):\(\displaystyle \mathrm{ENBW} = f_s\, \frac{\sum w[n]^2}{\left(\sum w[n]\right)^2} \tag{6}\) [Hz]。
電圧密度(ENBW法):\(\displaystyle S_{vv}[\mathrm{V}^2/\mathrm{Hz}] = \frac{|X[k]|^2}{N^2 \cdot \mathrm{ENBW}} \tag{7}\)
電圧密度(U正規化法):\(\displaystyle S_{vv}[\mathrm{V}^2/\mathrm{Hz}] = \frac{|X_w[k]|^2}{f_s \cdot N \cdot U} \tag{8}\)、ここで \(U = \frac{1}{N}\sum w^2\)
電力密度:\(\displaystyle S_{pp}[\mathrm{W}/\mathrm{Hz}] = \frac{S_{vv}}{R} = \frac{|X[k]|^2}{N^2 \cdot R \cdot \mathrm{ENBW}} \tag{9}\)(R=50 Ω等を明記) - dB表記の変換:電力系(dBm/Hz):\(P[\mathrm{W/Hz}] = S_{pp}[\mathrm{W}/\mathrm{Hz}]\) → **10 log₁₀**。電圧系(dBV/√Hz):\(\sqrt{S_{vv}}[\mathrm{V}/\sqrt{\mathrm{Hz}}]\) → **20 log₁₀**(R=50 Ω等)。
- 凡例に 窓名・CG・ENBW・平均化法・片側/二側・負荷R・参照(FS or 1 V)を併記(再現性確保)。1
1ENBW/CGの代表値は実装差あり(定義、係数、窓設計)。Hann 1.5Δf, Hamming 1.36Δf, Blackman 1.73Δf, Flat-top 3.77Δf等。
実務指針:凡例テンプレート
校正済みスペアナUI:凡例に窓名・CG・ENBW・平均法・片側/二側・負荷R・参照(FS or 1 V)を自動出力→実験ログの再現性/法規適合を確保。
例:「Hann窓, CG=0.5, ENBW=1.5Δf, Welch平均(M=16), 片側表示, R=50Ω, 参照=1V」
窓関数の代表値表
| 矩形窓 | 1.0 | 1.0 | 2.0 | ≈0.89 |
| Hann | 0.5 | 1.5 | 4.0 | ≈1.44 |
| Hamming | 0.54 | 1.36 | 4.0 | ≈1.30 |
| Blackman | 0.42 | 1.73 | 6.0 | ≈1.68 |
| Flat-top | 0.22 | 3.77 | 10.0 | >2.0(設計依存) |
単位変換早見表(R=50Ω基準)
| RMS → dBV | 20 log₁₀(V_rms) | 電圧基準(1V) |
| RMS → dBm | 20 log₁₀(V_rms) + 10 log₁₀(1000/R) | 電力基準(1mW) |
| dBV → RMS | 10^(dBV/20) | 電圧復元 |
| dBm → RMS | √(R/1000) × 10^(dBm/20) | 電圧復元 |
| 片側 → 二側 | 内点×2(DC/Nyquist非倍化) | スペクトル表示 |
| 二側 → 片側 | 内点÷2(DC/Nyquist非倍化) | スペクトル表示 |
正規化規約の比較
| DTFT | 無係数 | 1/(2π) | 理論的 |
| DFT(標準) | 無係数 | 1/N | 実装向け |
| DFT(ユニタリ) | 1/√N | 1/√N | 対称的 |
| FFT | DFTと同じ | DFTと同じ | 高速化のみ |
注意:この資料ではDFT(標準)規約を採用。ユニタリ規約ではParsevalの定理が \(\sum |x|^2 = \sum |X|^2\) となる。
ParsevalとPSDの帯域和の一致(検算)
有限長・採用正規化下での一貫性チェック:
- 時間ドメイン:$\sum_{n=0}^{N-1} |x[n]|^2$
- 周波数ドメイン:$\frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2$
- 窓掛けPSD帯域和:$\sum_{k=0}^{N-1} S_{vv}[k] \cdot \Delta f = \frac{f_s}{N} \sum_{k=0}^{N-1} \frac{|X_w[k]|^2}{N^2 \cdot \mathrm{ENBW} \cdot \mathrm{CG}^2}$
注意:Parsevalの定理により、時間ドメインと周波数ドメインでのエネルギーは等しいが、窓掛け後のPSD帯域和は窓エネルギー $U$ を含むため、一般に一致しない。
数値例(正弦波+白色雑音)
例:$x[n] = A\sin(2\pi f_0 n/N) + w[n]$($A=1$, $f_0=100$Hz, $N=1024$, $f_s=1024$Hz, $w[n]$は白色雑音)
- 時間ドメインエネルギー:$\sum_{n=0}^{N-1} |x[n]|^2 \approx N/2 + \sigma_w^2 N$
- 周波数ドメインエネルギー:$\frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2 \approx \frac{1}{2} + \sigma_w^2$
- 検算:$\sum |x|^2 \stackrel{?}{=} \frac{1}{N}\sum |X|^2$ → 両辺を$N$で割ると一致
注意:この例では $f_0 = f_s/10$ として、適切なサンプリング条件を満たしています。
実用的な応用
- 音響解析: 楽器の音色分析
- 通信技術: 信号品質の評価
- 医療機器: 生体信号の解析
- 産業応用: 機械の振動解析
振幅スペクトル
定義
信号の振幅スペクトルとは、フーリエ変換によって得られる各周波数成分ごとの「大きさ(振幅)」をグラフで表したものです。
特徴
- どの周波数が強く(大きな振幅を持ち)現れているかが一目でわかる
- 分析対象の信号がそれぞれの周波数成分ごとの「大きさ(振幅)」を表したグラフ
グラフの見方
- 横軸: 周波数(Hz)
- 縦軸: 振幅(フーリエ変換の大きさ $|X|$)
応用例
音声信号の振幅スペクトルを見ると、低い音(低周波数)や高い音(高周波数)がどれくらい含まれているかを詳細に分析できます。
振幅スペクトルと時間波形の比較
振幅スペクトル(周波数ドメイン)
- 縦軸: 振幅
- 横軸: 周波数
- ※注意: ここでの「振幅」はフーリエ変換の大きさ $|X|$ を指す。エネルギー/パワーを見たい場合は $|X|^2$(ESD/PSD)を用いる
時間波形(時間ドメイン)
- 縦軸: 信号値 $x(t)$ または $x[n]$
- 横軸: 時刻
- 確認方法: Python(NumPy/SciPy, matplotlib)や Audacity 等の汎用ツールで波形・スペクトルを再現可能
重要な違い
時間波形と振幅スペクトルでは、縦軸の意味が異なります:
- 時間波形(時間ドメイン)の縦軸: 信号値 $x(t)$ または $x[n]$
- 振幅スペクトル(周波数ドメイン)の縦軸: $|X(f)|$ または $|X[k]|$
- エネルギー/パワーを見るとき: $|X|^2$(パワースペクトル/PSD)を用いる
注: 時間波形の縦軸:$x(t)$ または $x[n]$。振幅スペクトルの縦軸:$|X(f)|$ または $|X[k]|$。エネルギー/パワーを見るときは $|X|^2$(パワースペクトル/PSD)を用いる。参考(パーセバル)$\int|x(t)|^2 dt = \frac{1}{2\pi}\int|X(\omega)|^2 d\omega$、$\sum|x[n]|^2 = \frac{1}{N}\sum|X[k]|^2$。
振幅スペクトルとエネルギースペクトルの関係
重要な区別:
- 振幅スペクトル: $|X(\omega)|$ または $|X[k]|$(フーリエ変換の大きさ)
- エネルギースペクトル: $|X(\omega)|^2$ または $|X[k]|^2$(パワースペクトル)
パーセバルの定理により、時間ドメインと周波数ドメインでのエネルギーは等しい:
連続時間:
$$\displaystyle \int_{-\infty}^{\infty} |x(t)|^2 \, dt = \frac{1}{2\pi} \int_{-\infty}^{\infty} |X(\omega)|^2 \, d\omega \tag{10}$$離散時間:
$$\displaystyle \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2 \tag{11}$$フーリエ級数(理論的基礎)
フーリエ級数の定義
周期Tの周期関数f(t)は、以下のように正弦波と余弦波の和で表すことができます:
$$\displaystyle f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} (a_n \cos(n\omega_0 t) + b_n \sin(n\omega_0 t))$$ここで、
- $\omega_0 = \frac{2\pi}{T}$:基本角周波数
- $a_0, a_n, b_n$:フーリエ係数
フーリエ級数の導出
ステップ1: 直交性の確認
まず、三角関数の直交性を確認します:
余弦関数の直交性:
$$\displaystyle \int_0^T \cos(n\omega_0 t) \cos(m\omega_0 t) \, dt = \begin{cases} \frac{T}{2} & (n = m \neq 0) \\ T & (n = m = 0) \\ 0 & (n \neq m) \end{cases}$$正弦関数の直交性:
$$\displaystyle \int_0^T \sin(n\omega_0 t) \sin(m\omega_0 t) \, dt = \begin{cases} \frac{T}{2} & (n = m \neq 0) \\ 0 & (n \neq m) \end{cases}$$余弦と正弦の直交性:
$$\displaystyle \int_0^T \cos(n\omega_0 t) \sin(m\omega_0 t) \, dt = 0 \quad \text{(すべてのn, mに対して)}$$ステップ2: フーリエ係数の導出
$a_0$の導出:
f(t)の両辺を[0, T]で積分:
直交性により、$n \geq 1$の項は0になるので:
$$\displaystyle \int_0^T f(t) \, dt = \int_0^T \frac{a_0}{2} \, dt = \frac{a_0}{2} \cdot T$$したがって:
$$\displaystyle a_0 = \frac{2}{T} \int_0^T f(t) \, dt$$$a_n$の導出:
f(t)に$\cos(m\omega_0 t)$を掛けて[0, T]で積分:
直交性により、$m=n$の場合のみ残る:
$$\displaystyle \int_0^T f(t) \cos(m\omega_0 t) \, dt = a_m \cdot \frac{T}{2}$$したがって:
$$\displaystyle a_n = \frac{2}{T} \int_0^T f(t) \cos(n\omega_0 t) \, dt$$$b_n$の導出:
同様に、f(t)に$\sin(m\omega_0 t)$を掛けて[0, T]で積分:
直交性により:
$$\displaystyle \int_0^T f(t) \sin(m\omega_0 t) \, dt = b_m \cdot \frac{T}{2}$$したがって:
$$\displaystyle b_n = \frac{2}{T} \int_0^T f(t) \sin(n\omega_0 t) \, dt$$複素フーリエ級数
複素フーリエ級数の導出
オイラーの公式を用いて、フーリエ級数を複素形式で表します:
ステップ1: オイラーの公式の適用
$$\displaystyle \cos(n\omega_0 t) = \frac{e^{jn\omega_0 t} + e^{-jn\omega_0 t}}{2}$$ $$\displaystyle \sin(n\omega_0 t) = \frac{e^{jn\omega_0 t} - e^{-jn\omega_0 t}}{2j}$$ステップ2: フーリエ級数の変形
$$\displaystyle f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} (a_n \cos(n\omega_0 t) + b_n \sin(n\omega_0 t))$$オイラーの公式を適用:
$$\displaystyle f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[a_n \frac{e^{jn\omega_0 t} + e^{-jn\omega_0 t}}{2} + b_n \frac{e^{jn\omega_0 t} - e^{-jn\omega_0 t}}{2j}\right]$$項を整理:
$$\displaystyle = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[\frac{a_n}{2} \cdot e^{jn\omega_0 t} + \frac{a_n}{2} \cdot e^{-jn\omega_0 t} + \frac{b_n}{2j} \cdot e^{jn\omega_0 t} - \frac{b_n}{2j} \cdot e^{-jn\omega_0 t}\right]$$ $$\displaystyle = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[\left(\frac{a_n}{2} + \frac{b_n}{2j}\right) \cdot e^{jn\omega_0 t} + \left(\frac{a_n}{2} - \frac{b_n}{2j}\right) \cdot e^{-jn\omega_0 t}\right]$$$\frac{1}{j} = -j$ の関係を使用:
$$\displaystyle = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[\frac{a_n - jb_n}{2} \cdot e^{jn\omega_0 t} + \frac{a_n + jb_n}{2} \cdot e^{-jn\omega_0 t}\right]$$ステップ3: 複素係数の定義
$$\displaystyle c_0 = \frac{a_0}{2}$$ $$\displaystyle c_n = \frac{a_n - jb_n}{2} \quad (n > 0)$$ $$\displaystyle c_{-n} = \frac{a_n + jb_n}{2} \quad (n > 0)$$ステップ4: 最終形
$$\displaystyle f(t) = \sum_{n=-\infty}^{\infty} c_n e^{jn\omega_0 t}$$ここで、
$$\displaystyle c_n = \frac{1}{T} \int_0^T f(t) e^{-jn\omega_0 t} \, dt$$導出の確認:
$$\displaystyle c_n = \frac{1}{T} \int_0^T f(t) e^{-jn\omega_0 t} \, dt$$ $$\displaystyle = \frac{1}{T} \int_0^T \left[\sum_{m=-\infty}^{\infty} c_m e^{jm\omega_0 t}\right] e^{-jn\omega_0 t} \, dt$$ $$\displaystyle = \frac{1}{T} \sum_{m=-\infty}^{\infty} c_m \int_0^T e^{j(m-n)\omega_0 t} \, dt$$$m=n$の場合:
$$\displaystyle \int_0^T e^{j(m-n)\omega_0 t} \, dt = \int_0^T 1 \, dt = T$$$m \neq n$の場合:
$$\displaystyle \int_0^T e^{j(m-n)\omega_0 t} \, dt = \left[\frac{e^{j(m-n)\omega_0 t}}{j(m-n)\omega_0}\right]_0^T = 0$$したがって:
$$\displaystyle c_n = \frac{1}{T} \cdot c_n \cdot T = c_n$$連続時間フーリエ変換
フーリエ変換の定義
連続時間信号 x(t) のフーリエ変換は、次の式で定義されます。
$$\displaystyle X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$$変数の意味
- $X(\omega)$:x(t) のフーリエ変換(周波数ドメインの表現)
- $x(t)$:連続時間の信号
- $t$:時刻(秒)
- $\omega$:角周波数(ラジアン/秒),$\omega = 2\pi f$
- $f$:周波数(Hz)
- $j$:虚数単位($j^2 = -1$)
- $e$:自然対数の底
重要な注意事項
- 積分範囲は $-\infty$ から $\infty$ まで
- 被積分関数は $x(t)$ と $e^{-j\omega t}$ の積
- $\omega$ は連続的な変数
フーリエ変換の導出
ステップ1: フーリエ級数からフーリエ変換への拡張
周期Tの信号を非周期信号に拡張する過程:
周期信号の場合:
$$\displaystyle f(t) = \sum_{n=-\infty}^{\infty} c_n e^{jn\omega_0 t}$$ $$\displaystyle c_n = \frac{1}{T} \int_{-T/2}^{T/2} f(t) e^{-jn\omega_0 t} \, dt$$非周期信号への拡張:
- $T \to \infty$ とすると、$\omega_0 = \frac{2\pi}{T} \to 0$
- 離散的な周波数 $n\omega_0$ が連続的な周波数 $\omega$ に変化
- フーリエ係数 $c_n$ が連続的な関数 $X(\omega)$ に変化
ステップ2: 数学的導出
ステップ2-1: 周期信号の表現
周期Tの信号 f(t) について:
複素フーリエ級数の定義:
$$\displaystyle f(t) = \sum_{n=-\infty}^{\infty} c_n e^{jn\omega_0 t}$$ここで、$c_n$は複素フーリエ係数:
$$\displaystyle c_n = \frac{1}{T} \int_0^T f(t) e^{-jn\omega_0 t} \, dt$$または、対称な積分範囲:
$$\displaystyle c_n = \frac{1}{T} \int_{-T/2}^{T/2} f(t) e^{-jn\omega_0 t} \, dt$$基本角周波数:
$$\displaystyle \omega_0 = \frac{2\pi}{T}$$ステップ2-2: 非周期信号への極限
$T \to \infty$ の極限を取る:
ステップ2-3: $c_n$を代入して整理
$c_n$の定義を代入します。ここで、f(t)の変数tと$c_n$の積分変数が衝突するため、$c_n$の積分変数は$\tau$に置き換えます:
$\frac{1}{T}$を分配:
$$\displaystyle = \lim_{T \to \infty} \sum_{n=-\infty}^{\infty} \left[\int_{-T/2}^{T/2} f(\tau) e^{-j\omega_n \tau} \, d\tau\right] e^{j\omega_n t} / T$$注意:変数$\tau$は積分の内部変数で、最終的な式f(t)の変数tとは区別されます。
ステップ2-4: 周波数変数の変換
$n\omega_0 = \omega_n$ とおき、$\Delta\omega = \omega_0 = \frac{2\pi}{T}$ とすると:
リーマン和の準備:
離散的な周波数間隔$\Delta\omega$を導入することで、総和を積分に変換する準備をします。$\frac{1}{T} = \frac{\Delta\omega}{2\pi}$の関係により:
$\Delta\omega/2\pi$の意味:
- $\Delta\omega = \frac{2\pi}{T}$ は周波数間隔
- $\frac{1}{T} = \frac{\Delta\omega}{2\pi}$ は正規化係数
- この係数により、総和が積分に変換される
ステップ2-6: リーマン和から積分へ
$T \to \infty$ のとき、$\Delta\omega \to 0$ となるので、これはリーマン和となる:
ここで、$X(\omega) = \lim_{T \to \infty} X_T(\omega) = \int_{-\infty}^{\infty} f(t) e^{-j\omega t} \, dt$ となる。
導出の確認:
$T \to \infty$ の極限では:
- $\omega_0 = \frac{2\pi}{T} \to 0$
- $n\omega_0 = n\frac{2\pi}{T} \to$ 連続変数 $\omega$
- $\Delta\omega = \omega_0 \to d\omega$
- $T \cdot c_n = X_T(\omega) \to X(\omega)$
- $\sum_{n=-\infty}^{\infty} \to \int_{-\infty}^{\infty}$
ステップ2-7: 収束条件の確認
フーリエ変換が存在するための条件:
- 信号$f(t)$が絶対可積分:$\int_{-\infty}^{\infty} |f(t)| \, dt < \infty$
- または、エネルギー有限:$\int_{-\infty}^{\infty} |f(t)|^2 \, dt < \infty$
- ディリクレ条件を満たす(有限個の不連続点、有限個の極値)
ステップ3: 最終的な定義
フーリエ変換:
$$\displaystyle X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$$フーリエ逆変換:
$$\displaystyle x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(\omega) e^{j\omega t} \, d\omega$$フーリエ逆変換
フーリエ逆変換は、周波数ドメインの信号 $X(\omega)$ を時間ドメインの信号 $x(t)$ に戻す変換です。
$$\displaystyle x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(\omega) e^{j\omega t} \, d\omega$$この式により、周波数スペクトルから元の時間信号を復元することができます。
逆変換の注意事項
- 正規化係数 $\frac{1}{2\pi}$ が重要
- 積分変数は $\omega$(角周波数)
- 被積分関数は $X(\omega)$ と $e^{j\omega t}$ の積
教科書による定義の違い
※教科書によっては、以下の式でフーリエ変換・逆変換を定義しているものもあります:
対称的な定義(第2の規約)
フーリエ変換:
$$\displaystyle X(\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$$フーリエ逆変換:
$$\displaystyle x(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} X(\omega) e^{j\omega t} \, d\omega$$この定義では、フーリエ変換と逆変換の両方に同じ正規化係数 $\frac{1}{\sqrt{2\pi}}$ が付きます。
重要な注意: どちらの定義を使うかは教科書や分野によって異なります。計算時は一貫して同じ定義を使うことが重要です。
離散時間フーリエ変換(DTFT)
DTFTの定義
離散時間信号 x[n] の離散時間フーリエ変換(DTFT)は、次の式で定義されます。
$$\displaystyle X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$$変数の意味
- $X(e^{j\omega})$:x[n] の離散時間フーリエ変換
- $x[n]$:離散時間の信号(nは整数)
- $n$:離散時間のインデックス(サンプル番号)
- $\omega$:正規化角周波数(ラジアン/サンプル)
- $j$:虚数単位($j^2 = -1$)
- $e$:自然対数の底
重要な注意事項
- 総和範囲は $-\infty$ から $\infty$ まで
- 被総和項は $x[n]$ と $e^{-j\omega n}$ の積
- $\omega$ は連続的な変数($-\pi \leq \omega \leq \pi$)
- $X(e^{j\omega})$ は周期 $2\pi$ の周期関数
DTFTの導出
ステップ1: 連続時間フーリエ変換から離散時間フーリエ変換への変換
連続時間フーリエ変換:
$$\displaystyle X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$$離散時間信号 $x[n]$ を考えます。これは連続時間信号 $x(t)$ をサンプリングしたものです:
$$\displaystyle x[n] = x(nT_s)$$ここで、$T_s$ はサンプリング周期です。
ステップ2: サンプリング定理の適用
サンプリング定理(ナイキスト・シャノンの定理):
- 連続信号を正確に復元するには、サンプリング周波数 $f_s$ が信号の最高周波数 $f_{max}$ の2倍以上である必要がある
- 条件:$f_s \geq 2f_{max}$
- ナイキスト周波数:$f_N = \frac{f_s}{2}$
エイリアシング現象:
- ナイキスト周波数を超える周波数成分は、低い周波数に折り返される
- これにより、元の信号と区別できない偽の周波数成分が生じる
ステップ3: 離散化による変換
連続時間フーリエ変換の積分を離散和に変換:
連続時間の場合:
$$\displaystyle X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$$離散時間の場合:
$$\displaystyle X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$$ここで、$\omega$ は正規化角周波数($\omega = \Omega T_s$)です。
ステップ4: 正規化角周波数の意味
正規化角周波数 $\omega$ の定義:
- $\omega = \Omega T_s = 2\pi f T_s = \frac{2\pi f}{f_s}$
- $\Omega$:連続時間の角周波数(ラジアン/秒)
- $f$:連続時間の周波数(Hz)
- $f_s$:サンプリング周波数(Hz)
- $T_s$:サンプリング周期(秒)
周波数範囲:
- 連続時間:$-\infty < \Omega < \infty$
- 離散時間:$-\pi \leq \omega \leq \pi$(または $0 \leq \omega < 2\pi$)
DTFTの性質
周期性
$X(e^{j\omega})$ は周期 $2\pi$ の周期関数です:
$$\displaystyle X(e^{j(\omega+2\pi)}) = X(e^{j\omega})$$証明:
$$\displaystyle X(e^{j(\omega+2\pi)}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j(\omega+2\pi)n}$$ $$\displaystyle = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} e^{-j2\pi n}$$ $$\displaystyle = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \cdot 1$$ $$\displaystyle = X(e^{j\omega})$$ここで、$e^{-j2\pi n} = 1$(nは整数)であることを利用しました。
対称性
実信号の場合、以下の対称性が成り立ちます:
$$\displaystyle X(e^{j\omega}) = X^*(e^{-j\omega})$$ここで、$*$ は複素共役を表します。
DTFTの逆変換
DTFTの逆変換は、次の式で与えられます:
$$\displaystyle x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega}) e^{j\omega n} \, d\omega$$この式により、周波数スペクトルから元の離散時間信号を復元することができます。
DTFTの注意事項
- 積分範囲は $-\pi$ から $\pi$ まで(または $0$ から $2\pi$ まで)
- 正規化係数 $\frac{1}{2\pi}$ が重要
- 積分変数は $\omega$(正規化角周波数)
- 被積分関数は $X(e^{j\omega})$ と $e^{j\omega n}$ の積
DTFTの応用
音響信号解析
DTFTは音響信号の周波数解析に広く用いられています:
- 楽器の音色分析: 各楽器の特徴的な周波数成分を特定
- 音声解析: 音声の基本周波数やフォルマント周波数の分析
- ノイズ解析: 不要な周波数成分の特定と除去
実装時の考慮事項
計算量の問題:
- DTFTは無限和を扱うため、直接計算は困難
- 実際の実装では、有限長の信号を仮定
- DTFTは連続$\omega$の数値評価:$K$ 点の周波数評価で $O(NK)$
- DFT/FFTは離散$\omega$:$K=N$ の離散周波数点を $O(N \log N)$ で計算
メモリ制約:
- 連続的な周波数変数 $\omega$ を離散化する必要
- 周波数解像度と計算量のトレードオフ
- 実用的にはDFT(離散フーリエ変換)を使用
離散フーリエ変換(DFT)
DFTの定義
N点離散フーリエ変換(DFT)は、次の式で定義されます。
$$\displaystyle X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn}$$変数の意味
- $X[k]$:x[n] の離散フーリエ変換(k番目の周波数成分)
- $x[n]$:離散時間の信号($n = 0, 1, \ldots, N-1$)
- $n$:離散時間のインデックス(サンプル番号)
- $k$:周波数インデックス($k = 0, 1, \ldots, N-1$)
- $N$:DFTの点数
- $j$:虚数単位($j^2 = -1$)
- $e$:自然対数の底
重要な注意事項
- 総和範囲は $0$ から $N-1$ まで
- 被総和項は $x[n]$ と $e^{-j\frac{2\pi}{N}kn}$ の積
- $k$ は離散的な変数($0 \leq k \leq N-1$)
- $X[k]$ は(拡張すれば)周期 $N$ の周期関数($X[k+N]=X[k]$)
DFTの導出
ステップ1: DTFTからDFTへの変換
DTFTの定義:
$$\displaystyle X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$$有限長信号 $x[n]$($n = 0, 1, \ldots, N-1$)を考えます。この場合、DTFTは:
$$\displaystyle X(e^{j\omega}) = \sum_{n=0}^{N-1} x[n] e^{-j\omega n}$$ステップ2: 周波数の離散化
周波数離散化の必要性:
- コンピュータで連続的な周波数変数 $\omega$ を扱うことは不可能
- 有限個の周波数点でのみ計算を行う必要
- 周波数解像度:$\Delta\omega = \frac{2\pi}{N}$
離散化された周波数:
$$\displaystyle \omega_k = k \cdot \frac{2\pi}{N} \quad (k = 0, 1, \ldots, N-1)$$ここで、$k$ は周波数インデックスです。
ステップ3: DFTの定義
離散化された周波数でのDTFT:
$$\displaystyle X[k] = X(e^{j\omega_k}) = \sum_{n=0}^{N-1} x[n] e^{-j\omega_k n}$$$\omega_k = k\frac{2\pi}{N}$ を代入:
$$\displaystyle X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn}$$ステップ4: 回転因子の導入
回転因子 $W_N$ を定義:
$$\displaystyle W_N = e^{-j\frac{2\pi}{N}}$$DFTは回転因子を用いて:
$$\displaystyle X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn}$$回転因子の性質
基本性質
回転因子の定義:
$$\displaystyle W_N = e^{-j\frac{2\pi}{N}} = \cos\left(\frac{2\pi}{N}\right) - j \sin\left(\frac{2\pi}{N}\right)$$オイラーの公式の利用:
$$\displaystyle W_N^{kn} = e^{-j\frac{2\pi}{N}kn} = \cos\left(\frac{2\pi kn}{N}\right) - j \sin\left(\frac{2\pi kn}{N}\right)$$重要な性質
周期性:
$$\displaystyle W_N^{k+N} = W_N^k$$対称性:
$$\displaystyle W_N^{N-k} = W_N^{-k} = (W_N^k)^*$$べき乗の性質:
$$\displaystyle W_N^{kn} = (W_N^k)^n$$DFTの逆変換
DFTの逆変換は、次の式で与えられます:
$$\displaystyle x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j\frac{2\pi}{N}kn}$$または、回転因子を用いて:
$$\displaystyle x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] W_N^{-kn}$$DFTの直交性
直交性の証明
DFTの直交性を証明します:
ステップ1: 直交性の定義
異なる周波数成分は直交します:
ステップ2: 証明
$k = m$ の場合:
$k \neq m$ の場合:
$$\displaystyle \sum_{n=0}^{N-1} e^{j\frac{2\pi}{N}(k-m)n} = \frac{1 - e^{j\frac{2\pi}{N}(k-m)N}}{1 - e^{j\frac{2\pi}{N}(k-m)}}$$ここで、$e^{j\frac{2\pi}{N}(k-m)N} = e^{j2\pi(k-m)} = 1$($k-m$は整数)なので:
$$\displaystyle = \frac{1 - 1}{1 - e^{j\frac{2\pi}{N}(k-m)}} = 0$$ステップ3: 逆変換の確認
逆変換の式を確認します:
$X[k]$ の定義を代入:
$$\displaystyle = \frac{1}{N} \sum_{k=0}^{N-1} \left[\sum_{m=0}^{N-1} x[m] e^{-j\frac{2\pi}{N}km}\right] e^{j\frac{2\pi}{N}kn}$$ $$\displaystyle = \frac{1}{N} \sum_{m=0}^{N-1} x[m] \sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}k(n-m)}$$直交性により:
$$\displaystyle = \frac{1}{N} \sum_{m=0}^{N-1} x[m] \cdot N \cdot \delta[n-m]$$ $$\displaystyle = \sum_{m=0}^{N-1} x[m] \delta[n-m]$$ $$\displaystyle = x[n]$$ここで、$\delta[n-m]$ はクロネッカーのデルタ関数です:
$$\displaystyle \delta[n-m] = \begin{cases} 1 & (n = m) \\ 0 & (n \neq m) \end{cases}$$DFTの性質
周期性
$X[k]$ は周期 $N$ の周期関数です:
$$\displaystyle X[k+N] = X[k]$$対称性
実信号の場合、以下の対称性が成り立ちます:
$$\displaystyle X[k] = X^*[N-k]$$ここで、* は複素共役を表します。実信号 $x[n]$ のとき $X[k]$ は共役対称($X[k]=X^*[N-k]$)。ただし $(k=0)$ と($N$ 偶数時の)$(k=N/2)$ は実数。
パーセバルの定理
時間ドメインと周波数ドメインでのエネルギーは等しい:
$$\displaystyle \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2$$DFTの応用
音響信号解析
DFTは音響信号の周波数解析に広く用いられています:
- 楽器の音色分析: 各楽器の特徴的な周波数成分を特定
- 音声解析: 音声の基本周波数やフォルマント周波数の分析
- ノイズ解析: 不要な周波数成分の特定と除去
実装時の考慮事項
周波数解像度:
- 周波数解像度:$\Delta f = \frac{f_s}{N} \tag{12}$
- $f_s$:サンプリング周波数
- $N$:DFT/FFTの点数
- より細かい周波数解析には大きな$N$が必要
- ゼロパディング:ピーク周波数内挿精度は視覚分解能は上がるが、CRLB観点の推定分散はSNR・窓・観測長で決まる(ゼロパッドは不変)。
統計的根拠:位相連続モデルでの周波数推定のCRLBは \(\operatorname{Var}(\hat f) \ge \frac{6}{(2\pi)^2 \cdot \mathrm{SNR} \cdot T^3} \tag{13}\) で、測定時間 \(T\) とSNRに支配される。ゼロパディングは \(T\) を伸ばさないため分散下限は不変。
誤解を明示:ゼロパディングは推定精度を向上させない。真の周波数推定には相位相補間(parabolic/quinn等)が必要。
高速フーリエ変換(FFT)
FFTとは
高速フーリエ変換(FFT)は、離散フーリエ変換(DFT)を効率的に計算するアルゴリズムです。
FFTの動機
DFTの計算量問題:
- DFTの計算量は $O(N^2)$
- $N=1024$の場合、約100万回の複素乗算が必要
- リアルタイム処理には計算時間が長すぎる
FFTの解決策:
- 典型的な基数2 Cooley–Tukey 法では、複素乗算は概ね $\tfrac{N}{2}\log_2 N$ 回、複素加算は $N\log_2 N$ 回
- 計算量は $O(N \log N)$
FFTの基本原理
分割統治法
FFTは分割統治法(Divide and Conquer)を使用します:
ステップ1: 信号の分割
$N$点の信号を2つの$N/2$点の信号に分割:
ステップ2: 部分DFTの計算
各$N/2$点信号のDFTを計算:
ステップ3: 結合
部分DFTを結合して全体のDFTを構成:
バタフライ演算
FFTの基本演算は「バタフライ演算」と呼ばれます:
$$\displaystyle X[k] = X_{even}[k] + W_N^k X_{odd}[k]$$ $$\displaystyle X[k+N/2] = X_{even}[k] - W_N^k X_{odd}[k]$$バタフライ演算の名前の由来:
- 信号の流れが蝶の形に似ている
- 2つの入力から2つの出力を生成
- FFTの基本構成要素
FFTの計算量解析
計算量の導出
ステップ1: 再帰関係
$N$点FFTの計算量を$T(N)$とすると:
ここで:
- $2T(N/2)$:2つの$N/2$点FFTの計算量
- $O(N)$:結合処理の計算量
ステップ2: 解の導出
再帰関係を解くと:
具体的な計算量比較(複素乗算回数の代表値)
- $N=1024$: $\tfrac{N}{2}\log_2 N = 5120$ 回
- $N=4096$: $\tfrac{N}{2}\log_2 N = 24576$ 回
- $N=16384$: $\tfrac{N}{2}\log_2 N = 114688$ 回
DFTの $O(N^2)$ と比較して、$N$ が大きいほど FFT の優位性は圧倒的になる。
FFTの応用
音響信号解析
FFTは音響信号の周波数解析に広く用いられています:
- 楽器の音色分析: 各楽器の特徴的な周波数成分を特定
- 音声解析: 音声の基本周波数やフォルマント周波数の分析
- ノイズ解析: 不要な周波数成分の特定と除去
実装時の考慮事項
リアルタイム処理:
- FFTによりリアルタイム周波数解析が可能
- オーバーラップ処理による連続解析
- 窓関数の適用によるスペクトル漏れの抑制
周波数解像度:
- 周波数解像度:$\Delta f = \frac{f_s}{N}$
- $f_s$:サンプリング周波数
- $N$:FFTの点数
- より細かい周波数解析には大きな$N$が必要
まとめ
学習の流れ
この資料では、フーリエ変換とスペクトル解析の基礎から応用までを体系的に学習しました:
- フーリエ級数 - 周期関数の基本
- 複素フーリエ級数 - オイラーの公式を活用
- フーリエ変換 - 非周期関数への拡張
- 離散時間フーリエ変換(DTFT) - 離散時間信号の周波数解析
- 離散フーリエ変換(DFT) - コンピュータ実装向け
- 高速フーリエ変換(FFT) - 効率的な計算アルゴリズム
重要な概念
数学的基礎:
- フーリエ変換の定義と導出
- サンプリング定理とエイリアシング
- 回転因子とオイラーの公式
- 直交性とパーセバルの定理
実用的な応用:
- 音響信号の周波数解析
- 楽器の音色分析
- 音声解析とノイズ除去
- リアルタイム信号処理
実装時の考慮事項
計算効率:
- FFTによる計算量の大幅削減(複素乗算 \(\tfrac{N}{2}\log_2 N\) 回程度)
- メモリ効率の最適化
- 数値精度の維持
周波数解析:
- 周波数解像度の調整
- 窓関数の適用と CG/ENBW の明示
- オーバーラップ処理
規約宣言
この資料で採用する規約:
- フーリエ変換:第1規約($X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt$、$x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(\omega) e^{j\omega t} \, d\omega$)
- DFT:正規化係数なし($X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn}$、$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j\frac{2\pi}{N}kn}$)
- 片側スペクトル:内点×2(DC/Nyquist非倍化)
- 単位系:$\sum w^2$ は窓の二乗和、$\sum w$ は窓和、$\mathrm{ENBW} = f_s \cdot \frac{\sum w^2}{(\sum w)^2}$。ParsevalとENBWの対応を確認。
参考文献
一次文献:
- F. J. Harris, "On the use of windows for harmonic analysis with the DFT," Proc. IEEE, vol. 66, no. 1, pp. 51-83, 1978. (窓関数とENBW/CGの基礎)
- P. D. Welch, "The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms," IEEE Trans. Audio Electroacoust., vol. 15, no. 2, pp. 70-73, 1967. (Welch平均法)
- J. W. Cooley & J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series," Math. Comput., vol. 19, no. 90, pp. 297-301, 1965. (FFTアルゴリズム)
- A. V. Oppenheim & R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Prentice Hall, 2010. (学習者向けに最適)
応用・転用例:
- 校正済みスペアナUI:凡例に「窓・CG・ENBW・平均法・片側/二側・負荷R」を自動出力→実験ログの再現性確保(計測工学、品質保証、法規適合)
- ピーク推定器:Hann窓×相位相補間(parabolic/quinn等)+CRLB境界表示→音響F0推定・モータ振動診断・通信キャリアオフセット推定(信号処理、機械保全、通信工学)
- PSD規格変換モジュール:dBFS↔dBV↔dBm/Hz間をENBWと負荷で可逆変換→計測器間の整合(計量管理、計装、テスト自動化)
コメント