2020/06/07

FFTの周波数ビンからずれた正弦波のパワースペクトルとサイドローブとか

FFT（とDFT）を使用した信号処理のコードを書いていて、半端な周波数をFFTしたらどういった結果になるのか試してみたので纏めました。

よくあるサンプリングレート44.1kHzでサンプリングされた音声波形を、窓長4096サンプルでFFTした結果は、44100/4096 = 約10.76[Hz] の周波数ビン単位で離散化された各周波数(0 , 10.76, 21.52,... ,22050[Hz])に分解されて表現されます。

ここで入力信号として正弦波をFFTする場合で、正弦波の周波数がちょうど周波数ビンと一致するケース、
・1001.29...[Hz]（=10.76*93 ちょうど）
・1012.06...‬[Hz]（=10.76*94 ちょうど）
では良いですが、周波数ビンに一致しない中途半端な周波数
・1006.67...[Hz]（=10.76*93.5）
の時はどうなるんだったか、隣接する周波数ビンにパワーが分散するんだか、昔デジタル信号処理とかで学んだような記憶もありますが既にうろ覚えです。

FFTの周波数ビンからずれた正弦波のパワースペクトルは誤差が大きい

というわけで、100ビンから101ビンまで0.25周波数ビン刻みの周波数の正弦波（振幅1.0 = 0[dB]）をFFTしてみたのが以下の結果です。（44.1kHz,4096sample,Hanning窓）
周波数ビンからずれた正弦波のFFT結果
ちょうど100bin,101binの正弦波をFFTした場合、メインローブは±1binに収まり、サイドローブもグラフ範囲外のかなり低いレベルとなっています。一方、少しでも周波数ビンからずれた周波数の正弦波をFFTするとスペクトルが広がっています。（スペクトル・リーク、Spectral leakage）
ずれ量が最大となる0.5ビン差の場合が最もスペクトルが広がります。


周波数ビンからずれた正弦波のFFT結果
ピークパワーについても、周波数ビンちょうどの場合は0dBですが、0.5ビンずれた周波数をFFTした結果は-1.4dBまで低下しています。これはスカラップ・ロス（Scalloping loss）によるものです。

FFTに使用する窓関数の周波数特性を計算する

FFTの周波数分解能やサイドローブの大きさは、窓関数に依存します。（周波数ビンからのずれの影響も窓関数によって異なります）
今回使用した窓関数はハニング窓（ハン窓、Hanning window)なので、ハニング窓の周波数特性を計算します。
ハニング窓の周波数特性
以下のコード（java)でハニング窓関数の周波数特性を計算しました。

/**
 * Hanning窓の周波数特性を計算する
 * @param sampleRate_Hz
 * @param fftSampleLength
 * @param frqOffset_Hz 周波数ビンと一致するなら0
 * @return 振幅（0.0-1.0）
 */
public static double calcMainLobeByHanningWindow(final float sampleRate_Hz,final float fftSampleLength,float frqOffset_Hz) {
  final float len = fftSampleLength/sampleRate_Hz;
  if(frqOffset_Hz == 0f) {
    return 1.0;
  }
  final double w = 2.0*Math.PI*frqOffset_Hz;
  double amp;
  if(frqOffset_Hz == sampleRate_Hz/fftSampleLength || frqOffset_Hz == -sampleRate_Hz/fftSampleLength) {
    amp = len/4.0;
  }else {
    amp = ((1.0/w) - ((len*len*w)/((len*len*w*w) - (4.0*Math.PI*Math.PI))))*Math.sin((len*w)/2.0);
  }
  if(amp<0) {
    return -amp/(len/2);
  }
  return amp/(len/2);
}

// 上記関数で得られるのは振幅なので、10*Math.log10(amp*amp) でパワー(dB)に換算しています。

計算式参考元：Allisone：窓関数の周波数特性

ハニング窓の周波数特性とスペクトルリークの関係
計算したハニング窓のスペクトルを見ると、メインローブの左右にサイドローブが広がっていますが、各ローブ間のディップがちょうど周波数ビンと一致しています。
このため、周波数ビンの整数倍と完全に一致する周波数の正弦波をFFTした場合、FFT結果のメインローブは三点（上図中央の緑の点）となり、更にサイドローブはディップと一致するので非常に低いレベル（グラフ範囲外）となります。

一方、周波数ビンの中間に位置する半端な周波数の場合は、メインローブは中央のピークを外れた4点（上図中央の橙の点）となり、スカラップ・ロスが発生してピーク値が若干下がって、またメインローブも広がります。サイドローブも0.5ビンずれた位置で観測されるのでちょうどサイドローブのピークの部分の値（上図左右の橙の点）となり、スペクトル・リークが大きくなります。
周波数ビンからずれた正弦波のFFT結果
FFTにおける窓関数の必要性

そもそもFFTは連続波形の離散時間信号の一部を有限長のサンプル数で切り出し、無限に繰り返す波形と見做してフーリエ変換するわけですが、「切り出すサンプル数」と「サンプル中に含まれる全ての周波数成分の周期」が完全に一致していないと、サンプルの両端を並べて繋ぎ合わせた際に必ず不連続な部分が発生してしまいます。

普通の音声信号は複数の異なる周波数成分が混ざり合っているので、サンプル中に含まれるすべての周波数成分の周期単位で完全一致するようにサンプルを切り出すことは不可能です。（更に、DFTではなく高速なFFTを使用するにはサンプル数は2^Nに限定されますし、大抵は信号に含まれる未知の周波数成分を知りたいのでFFTする訳です）

結果として任意サンプル長で切り出した波形を繰り返し並べた際に生じる「本来の波形に存在しない不連続な部分」もフーリエ変換時に各周波数成分に分解され、スペクトル上に信号以外の周波数成分がメインローブ、サイドローブの形で現れてしまいます。
その結果周波数分解能やダイナミックレンジが低下することになるので、これを抑えるため（繋ぎ目の不連続さを低減するため）切り出したサンプルの両端で0になるような窓関数を掛けて、サンプルを繰り返し波形と見做して並べた際の繋ぎ目で不連続さが目立たないようにしています。（しかし窓関数を使用しても周波数分解能とダイナミックレンジはトレードオフの関係になってしまうので最強の窓関数は存在しない）

「FFTの周波数ビンの整数倍の周波数の正弦波」は特別

結局、なぜ周波数ビンの整数倍の周波数の正弦波をFFTした場合だけメインローブが細くなりスペクトルリークも無視できるほど抑えられるかというと、
「周波数ビン(サンプリングレート/FFTサンプル数)[Hz]のN倍の周波数の正弦波」の１周期は(FFTサンプル数/(N*サンプリングレート))[Sec]、これにサンプリングレートを掛けて1周期の長さをサンプル数に直すと、(FFTサンプル数/N)[Sample]となります。

つまり「FFTの周波数ビンのN倍の周波数の正弦波」は、FFTのサンプル長に綺麗にN個収まります。

結果的に「FFTの周波数ビンの整数倍の周波数の正弦波」は、FFTサンプル数で切り出しても両端で全く不連続さが発生しない特別な場合となり、メインローブは細くなり、サイドローブも発生しなくなる訳です。よくよく考えてみれば当たり前の事でした。
時刻: