ホームに戻る
 FFTによる周波数解析

FFTは日本語で「高速フーリエ変換」のことです。

フーリエ変換の説明をする前に、まず波形の説明から、
通常、音や電気などの信号は波形で表すことができます。
周波数が一定の波は、ずっと一定の音程で音が鳴り続けます。
このとき、高い音は周波数が高いですし、低い音は周波数が低いです。
自然の波形というものは、いろんな周波数の波が混ざりあっています。
混ざり合うと、波は足し算されるので見た目は不規則な波になって、
どの周波数の波が混ざりあっているのかがわからなくなります。
混ざりあった結果の波形を目でみることができたとしても、
どの周波数の波がどれくらい混ざっているか分解して見ることはできません。

フーリエ変換はこれを可能にします。
実際には各周波数のサイン波に分解し、現在の位相までも解析します。
原理は説明しませんが、とにかく数式に当てはめればOKです。
(数式コードの探し方については後述)
例えば、サンプルレート(単位時間当たりのデータ数)が fs のとき、
FFTでは fs / 2 までの周波数解析が可能です。
データ数を N とし実数部 R[N] と虚数部 I[N] を準備します。
R[N] には解析前のデータが入っており、I[N] はすべて 0 です。
これをFFTにかけると R[N] と I[N] ともに左右対称の結果が得られます。
左右対称ですので片方あればもう一方に価値はありません。
解析後の N / 2 のデータから N / 2 個のサイン波が得られます。
サイン波は以下のようです、k は 0 から N / 2 までの整数です。

((R[k] * R[k] + I[k] * I[k]) ^ (1 / 2)) * sin(2 * PAI * (fs / N) * k) + arctan(I[k] / R[k]))

強度:((R[k] * R[k] + I[k] * I[k]) ^ (1 / 2))
位相:sin(2 * PAI * (fs / N) * k) + arctan(I[k] / R[k]))

制限として、解析できる周波数はどうしても間隔が fs / N 単位になります。
つまり、周波数 fk = (fs / N) * k の情報は得ることができます。
では、その間の周波数の結果はというと、最寄の fk に移動してしまいます。
N はできるだけ大きいほうが詳しく解析できることがわかります。

FFTを2回行うことで正確な周波数を特定する方法があります。
1回のFFTでは正確な周波数がたとえ fk + offset であったとしても、
FFTで得られる情報は周波数 fk へ移動してしまいました。 
ここで既知の時間での位相の変化に注目します。
例えば、既知の時間間隔において、
fk の周波数で位相の変化を調べることは可能です。
ただし、この周波数にある結果が実際には fk + offset のものであるとすると、
位相は fk にありながら周波数 fk + offset で位相変化します。
ここで fk の位相変化と fk + offset の位相変化のズレから、
offset をもとめることが可能です。
これでより正確な周波数を確定することができます。

実際にやってみると以下の通り。

サンプルレート 44100 Hz FFTサイズ 512 で行う。
OVERLAP は 8 とする。
OVERLAP とは2回目測定までの間隔を定める。
OVERLAP を 8 とするとき、2回目測定は 512 / 8 = 64 サンプル後とする。

FFT を行うとき結果の間隔は、

44100 / 512 = 86.1328 Hz ごとになる。

ここで 1000 Hz のサイン波をFFTにかけたとき、

86.1328 * 11 = 947.4609 に結果が出たとしよう。

このときFFTの結果は以下の通り

1回目の位相 0.340390
2回目の位相 -3.105986

よって位相のズレは -3.105986 - 0.340390 = -3.446376

さらに 64 サンプル進むのでその時間での位相のズレも計算に入れる。

例えば周波数が 86.1328 Hz であれば 512 サンプルで 2 * 3.14 だけ位相が進む。
よって、周波数が 86.1328 * 11 = 947.4609 Hz で 64 サンプルにおいてはすなわち、

(11 * 2 * 3.14) / (512 / 64) = 8.635 だけ位相が進む。

よって、 最終的な位相のズレは -3.446376 - 8.635 = -12.081376 となる。

これを -12.081376 * 2 * 6.28 = 0.478624 とし -3.14 から 3.14 の間におさめる。

FFTの結果軸が1本違うごとに 512 サンプルの移動で 2 * 3.14 位相がズレる。

よって、 11 + (0.478624 / (2 * 3.14)) のオフセットのズレがあるといえる。

ただし、今の例では OVERLAP が 8 であるので、

オフセットのズレは 11 + (8 * 0.478624 / (2 * 3.14)) となる。

よって、正確な周波数は、

86.1328 * (11 + (8 * 0.478624 / (2 * 3.14))) = 999.9770 Hz

となりほぼ 1000 Hz となる。


最後に、窓の話。
以上の方法で計算してみると、意外に結果がバラつきます。
原因はFFTが繰り返しても連続な波形を想定している点にあります。
すなわち解析領域を繰り返したときに波形が切れないほうが良いのです。
ここで「窓関数」というのを使います。
窓関数とは中央が 1 で両端が 0 になっている関数です。
コサイン窓やハミング窓などがありますが、
とにかく両端が収束していれば問題ありません。
計算前にデータに窓関数を掛け算しておくのです。
こうすれば繰り返しても連続する波形として解析が可能です。
この方法である程度バラつきを解決できます。

さて、高速フーリエ変換の「高速」はというと、
解析するデータ数を2の累乗にすると計算が速くなるというものです。
FFTのプログラムはたくさん公開されていますので、
書籍でもWebでもどこでも見つかると思います。

以下のサンプルはFFTのコードをどこかから拝借しており、
ハミング窓を使用し、2回解析による正確な周波数を求めています。
窓にもいろいろあるので探してみてください。

サンプルの結果は、

解析周波数(2回解析による補正周波数): 2回目解析での強度, 2回目解析での位相

となります。

N、OVERLAP、SAMPLERATE などを変えてみる。
WAVEを変えてみる。窓を変えてみる。
などの変更を試してみると良いのではないでしょうか。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.14159265358979323846

/*
/  サイン波を作る
/  ポインタのアドレスを渡すと実体を作ってアドレスを入れる
/  ポインタには NULL を入れておいてください
/  buffer_size はサンプル数 hz は振動数 samplerate はサンプルレート
*/
int make_sinwave(float **pwave, int buffer_size, int hz, int samplerate){
  int i;
  float *p;

  if(*pwave != NULL){
    free(*pwave);
  }
  p = (float *)malloc(buffer_size * sizeof(float));

  if(p == NULL)return 0;

  for(i = 0; i < buffer_size; i++){
    p[i] = (float)sin(2 * PI * i * hz / samplerate);
  }

  *pwave = p;

  return 1;
}

/*
/  FFT の結果よりスペクトルの強度とその周波数を得る
/  n:サンプル数 (0 を入れるとバッファを解放する)
/  x[]:FFT の結果の実数部分
/  y[]:FFT の結果の虚数部分
/  magni[]、freq[]:結果を入れる それぞれ強度と周波数
*/
int getMagnitudeAndFreq(int n, float x[], float y[], float magni[], float phase[]){
  int i;

  for(i = 0; i < n / 2; i++){
    magni[i] = (float)sqrt(x[i] * x[i] + y[i] * y[i]);
    phase[i] = (float)atan2(y[i], x[i]);
  }

  return 1;
}

/*
*  ハミング窓を作成
*/
void make_hamming_wintbl(int n, float wintbl[]){
  int i;

  for(i = 0; i < n; i++){
    wintbl[i] = 0.54f - 0.46f * (float)cos((2.0f * PI * i) / (n - 1));
  }
}

/*
/  ここから「C言語による最新アルゴリズム事典」より抜粋
*/

/*
  関数{\tt fft()}の下請けとして三角関数表を作る.
*/
static void make_sintbl(int n, float sintbl[])
{
  int i, n2, n4, n8;
  double c, s, dc, ds, t;

  n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
  t = sin(PI / n);
  dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
  t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
  for (i = 1; i < n8; i++) {
    c -= dc;  dc += t * c;
    s += ds;  ds -= t * s;
    sintbl[i] = s;  sintbl[n4 - i] = c;
  }
  if (n8 != 0) sintbl[n8] = sqrt(0.5);
  for (i = 0; i < n4; i++)
    sintbl[n2 - i] = sintbl[i];
  for (i = 0; i < n2 + n4; i++)
    sintbl[i + n2] = - sintbl[i];
}

/*
  関数{\tt fft()}の下請けとしてビット反転表を作る.
*/
static void make_bitrev(int n, int bitrev[])
{
  int i, j, k, n2;

  n2 = n / 2;  i = j = 0;
  for ( ; ; ) {
    bitrev[i] = j;
    if (++i >= n) break;
    k = n2;
    while (k <= j) {  j -= k;  k /= 2;  }
    j += k;
  }
}

/*
  高速Fourier変換 (Cooley--Tukeyのアルゴリズム).
  標本点の数 {\tt n} は2の整数乗に限る.
  {\tt x[$k$]} が実部, {\tt y[$k$]} が虚部 ($k = 0$, $1$, $2$,
  \ldots, $|{\tt n}| - 1$).
  結果は {\tt x[]}, {\tt y[]} に上書きされる.
  ${\tt n} = 0$ なら表のメモリを解放する.
  ${\tt n} < 0$ なら逆変換を行う.
  前回と異なる $|{\tt n}|$ の値で呼び出すと,
  三角関数とビット反転の表を作るために多少余分に時間がかかる.
  この表のための記憶領域獲得に失敗すると1を返す (正常終了時
  の戻り値は0).
  これらの表の記憶領域を解放するには ${\tt n} = 0$ として
  呼び出す (このときは {\tt x[]}, {\tt y[]} の値は変わらない).
*/
int fft(int n, float x[], float y[])
{
  static int    last_n = 0;    /* 前回呼出し時の {\tt n} */
  static int   *bitrev = NULL; /* ビット反転表 */
  static float *sintbl = NULL; /* 三角関数表 */
  int i, j, k, ik, h, d, k2, n4, inverse;
  float t, s, c, dx, dy;

  /* 準備 */
  if(n < 0){
    n = -n;
    inverse = 1;  /* 逆変換 */
  }
  else{
    inverse = 0;
  }

  n4 = n / 4;

  if (n != last_n || n == 0) {
    last_n = n;
    if (sintbl != NULL) free(sintbl);
    if (bitrev != NULL) free(bitrev);
    if (n == 0) return 0;  /* 記憶領域を解放した */
    sintbl = malloc((n + n4) * sizeof(float));
    bitrev = malloc(n * sizeof(int));
    if (sintbl == NULL || bitrev == NULL) {
      fprintf(stderr, "記憶領域不足\n");  return 1;
    }
    make_sintbl(n, sintbl);
    make_bitrev(n, bitrev);
  }
  for (i = 0; i < n; i++) {    /* ビット反転 */
    j = bitrev[i];
    if (i < j) {
      t = x[i];  x[i] = x[j];  x[j] = t;
      t = y[i];  y[i] = y[j];  y[j] = t;
    }
  }
  for (k = 1; k < n; k = k2) {    /* 変換 */
    h = 0;  k2 = k + k;  d = n / k2;
    for (j = 0; j < k; j++) {
      c = sintbl[h + n4];
      if (inverse) s = - sintbl[h];
      else         s =   sintbl[h];
      for (i = j; i < n; i += k2) {
        ik = i + k;
        dx = s * y[ik] + c * x[ik];
        dy = c * y[ik] - s * x[ik];
        x[ik] = x[i] - dx;  x[i] += dx;
        y[ik] = y[i] - dy;  y[i] += dy;
      }
      h += d;
    }
  }
  if (! inverse)    /* 逆変換でないならnで割る */
    for (i = 0; i < n; i++) {  x[i] /= n;  y[i] /= n;  }
  return 0;  /* 正常終了 */
}

/*
/  ここまで「C言語による最新アルゴリズム事典」より抜粋
*/

#define N 512            // N 個のデータで解析 2の累乗であること
#define OVERLAP 8        // 2回目の解析はサンプルを N / OVERLAP だけ進めてから解析
#define SAMPLERATE 44100 // サンプルレート

int main()
{
  int i, s;
  static float win[N], x[N], y[N], magni[N / 2], phase[N / 2], old_phase[N / 2];
  float *pwave = NULL;

  // 4649Hz の波形を N + N / OVERLAP 個作成
  make_sinwave(&pwave, N + N / OVERLAP, 4649, SAMPLERATE);

  // ハミング窓を作成
  make_hamming_wintbl(N, win);

  // データの準備
  for (i = 0; i < N; i++){
    x[i] = pwave[i];
    x[i] *= win[i]; // ハミング窓をかける
    y[i] = 0;
  }

  // FFT
  if(fft(N, x, y))return EXIT_FAILURE;

  // FFT の結果を強度と周波数に分解
  getMagnitudeAndFreq(N, x, y, magni, phase);

  // 前回の位相を記録しておく
  for(i = 0; i < N / 2; i++){
    old_phase[i] = phase[i];
  }

  // 波形を N / OVERLAP だけ進めます
  for(i = 0; i < N; i++){
    x[i] = pwave[i + N / OVERLAP];
    x[i] *= win[i]; // ハミング窓をかける
    y[i] = 0;
  }

  // FFT
  if(fft(N, x, y)) return EXIT_FAILURE;

  // FFT の結果を強度と周波数に分解
  getMagnitudeAndFreq(N, x, y, magni, phase);

  for(i = 0; i < N / 2; i++){
    float OffsetPhase = (phase[i] - old_phase[i]) - (i * 2.0f * PI / OVERLAP); // (FFT計算結果の位相変化) - (SAMPLERATE / N 間隔での位相変化)

    // 以降 -PI < OffsetPhase < PI にするための処理
    s = (int)(OffsetPhase / PI); 
    if(s >= 0){
      s += s & 1;
    }
    else{
      s -= s & 1;
    }
    OffsetPhase -= PI * (float)s;
    // ここまで -PI < OffsetPhase < PI にするための処理

    // 位相のズレから正確な周波数を算出
    float Offset = OVERLAP * OffsetPhase / (2.0 * PI);
    float freq = (float)((i + Offset) * SAMPLERATE / N);

    // 解析周波数(2回解析による補正周波数): 2回目解析での強度, 2回目解析での位相
    printf("%012f (%012f): %012f, %012f\n", (float)SAMPLERATE / N * i, freq, magni[i], phase[i]);
  }

  getchar();

  // メモリの解放
  free(pwave);

  return EXIT_SUCCESS;
}

inserted by FC2 system