スペクトル解析
Note
本節のポイント
- FFTの出力結果を物理的な「周波数」と「強度」に対応させる方法を学ぶ。
- スペクトル漏れを防ぐための窓関数の役割と、
ndarrayを用いた適用方法を習得する。plottersを用いて、計算結果を可視化する。
FFTを計算しただけでは、物理的な意味を持つ「周波数」は分かりません。FFTの出力結果をどのように解釈し、実際の信号を分析するか（スペクトル解析）について学びます。
周波数軸との対応
サンプリング周期を （サンプリング周波数 ）、データ数を とします。 FFTの出力 の各インデックス に対応する物理的な周波数 は以下のようになります。
ここで、 が正の周波数成分、 が負の周波数成分に対応します。 ナイキスト周波数（再現可能な最高周波数）は です。
パワースペクトル(Power Spectrum)
各周波数成分の「強さ」を見るために、複素数 の絶対値の2乗をとったものをパワースペクトルと呼びます。
物理学では、エネルギーがどの周波数帯域に分布しているか（パワースペクトル密度, PSD）を調べることで、系の振動特性などを明らかにします。
窓関数(Window Function)
FFTは、入力信号が「周期 で無限に繰り返されている」ことを仮定しています。 実際の有限な信号をそのまま切り取ると、信号の端点で不連続が生じ、スペクトル漏れ (Spectral Leakage) というノイズが発生します。
これを防ぐために、信号の両端を滑らかに
に落とす窓関数を事前に乗算します。ndarrayを使えば、信号ベクトルと窓関数ベクトルの要素ごとの積として簡潔に記述できます。
代表的な窓関数：
- Hann窓(Hann window, Hanning window とも呼ばれる): 最も一般的に使われる窓関数で、周波数分解能とスペクトル漏れの抑制のバランスが良い。
- ハミング窓(Hamming window): Hann窓に似ているが、端点が になりきらない。周波数分解能を少し優先したい場合に用いられる。
- ブラックマン窓(Blackman window): Hann窓よりもスペクトル漏れを強力に抑えるが、メインローブが広くなり周波数分解能は低下する。
実践例：合成信号の解析と可視化
2つの異なる周波数のサイン波を含む信号をFFTし、その周波数を特定してみましょう。数値計算結果の可視化で触れたplottersを用いた可視化も行います。
依存関係
[dependencies]
rustfft = "6.4"
ndarray = "0.17"
num-complex = "0.4"
plotters = "0.3"
実装例
use rustfft::FftPlanner;
use ndarray::Array1;
use num_complex::Complex64;
use std::f64::consts::PI;
use plotters::prelude::{
BitMapBackend, ChartBuilder, IntoDrawingArea, LineSeries, RED, WHITE,
};
fn time_grid(n: usize, dt: f64) -> Array1<f64> {
Array1::from_shape_fn(n, |i| i as f64 * dt)
}
fn synthetic_signal(times: &Array1<f64>) -> Array1<Complex64> {
times.mapv(|ti| {
let val = 1.0 * (2.0 * PI * 50.0 * ti).sin()
+ 0.5 * (2.0 * PI * 120.0 * ti).sin();
Complex64::new(val, 0.0)
})
}
fn hann_window(n: usize) -> Array1<f64> {
Array1::from_shape_fn(n, |i| {
0.5 * (1.0 - (2.0 * PI * i as f64 / (n as f64 - 1.0)).cos())
})
}
fn apply_real_window(signal: &mut Array1<Complex64>, window: &Array1<f64>) {
signal.zip_mut_with(window, |value, &weight| {
*value *= Complex64::new(weight, 0.0);
});
}
fn fft_in_place(signal: &mut Array1<Complex64>) {
let mut planner = FftPlanner::new();
let fft = planner.plan_fft_forward(signal.len());
fft.process(signal.as_slice_mut().expect("Array must be contiguous"));
}
fn positive_frequency_power(signal: &Array1<Complex64>, fs: f64) -> (Vec<f64>, Vec<f64>) {
let n = signal.len();
let freqs: Vec<f64> = (0..n / 2).map(|k| k as f64 * fs / n as f64).collect();
let powers: Vec<f64> = signal.iter().take(n / 2).map(|c| c.norm_sqr()).collect();
(freqs, powers)
}
fn plot_power_spectrum(
freqs: &[f64],
powers: &[f64],
fs: f64,
path: &str,
) -> Result<(), Box<dyn std::error::Error>> {
let y_max = powers.iter().copied().fold(0.0_f64, f64::max);
let root = BitMapBackend::new(path, (800, 600)).into_drawing_area();
root.fill(&WHITE)?;
let mut chart = ChartBuilder::on(&root)
.caption("Power Spectrum", ("sans-serif", 30))
.margin(10)
.x_label_area_size(40)
.y_label_area_size(50)
.build_cartesian_2d(0.0..fs / 2.0, 0.0..y_max)?;
chart.configure_mesh()
.x_desc("Frequency [Hz]")
.y_desc("Power")
.draw()?;
chart.draw_series(LineSeries::new(
freqs.iter().copied().zip(powers.iter().copied()),
&RED,
))?;
root.present()?;
Ok(())
}
fn main() -> Result<(), Box<dyn std::error::Error>> {
let n = 1024;
let fs = 1000.0; // サンプリング周波数 1000Hz
let dt = 1.0 / fs;
// 1. 信号の生成(50Hz + 120Hz)
let t = time_grid(n, dt);
let mut signal = synthetic_signal(&t);
// 2. 窓関数の適用(Hann窓)
let window = hann_window(n);
apply_real_window(&mut signal, &window);
// 3. FFTの実行
fft_in_place(&mut signal);
// 4. 結果の解析(パワースペクトルの計算)
let (freqs, powers) = positive_frequency_power(&signal, fs);
// 5. 可視化
plot_power_spectrum(&freqs, &powers, fs, "spectrum.png")?;
println!("spectrum.png を生成しました。");
Ok(())
}
まとめ
- FFTのインデックス は周波数 に対応する。
- データの端点による誤差（スペクトル漏れ）を抑えるには窓関数が不可欠。
ndarrayのベクトル演算を活用することで、信号処理のコードを簡潔かつ効率的に記述できる。- パワースペクトルを可視化することで、信号に含まれる支配的な周波数成分を一目で特定できる。
参考リンク
本章はこれで終わりです。次は常微分方程式に進みましょう。