高速フーリエ変換 FFT

高速フーリエ変換 FFT アルゴリズムはデジタル信号処理に果した役割にははかりしれないものがあります。

Maple にも1次元のFFTプログラムが用意されています。ここではこれを利用した例を示しFFTの価値を改めて見直しましょう。

>    restart;
with(linalg):
with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

 波の合成

いまギターの弦が 1 秒間に 442 回振動しているとします。

Maple で遅くしたアニメーションをみてみましょう。

>    animate(cos(2*Pi*t)*sin(x),x=0..Pi,t=0..1);

[Maple Plot]

弦の中心が時間と共にどのように変化するかをプロットしてみます (0.01 秒間)。

>    plot(sin(2*442*Pi*t),t=0..0.01);

[Maple Plot]

440 ヘルツの波 (1 秒間に 440 回振動) と足したらどうなるでしょうか。(1 秒間)

>    plot(sin(2*442*Pi*t)+sin(2*440*Pi*t),t=0..1,numpoints=200);

[Maple Plot]

高周波なのでプロットはきれいではありませんが”唸り”がよく分かります。このように波を合成(足しあわせる)すると
もとの波は想像できないものになります。しかしフーリエ変換は元の信号を探し出すことを可能にします。
FFTアルゴリズムは計測したデジタルデータからその周波数成分を導きだしてくれます

FFTの使い方

FFT ではまずデータの観測時間を限定して処理されます。観測時間をデータ数で (正確にはデータ数 -1) 割ったものがサンプリング時間です。

通常FFTは 2 のベキ (4, 8, 16, 32, 64, ... , 2^n) のデータ数に対して計算します。
FFTを呼び出すときこの指数nを指定します

>    n:=4;

n := 4

n が 4 であれば 16 個のデータ数です。また呼び出すには実部と虚部のベクトルを使用します。
サンプルのデータを準備してみましょう。

>    d_real:=vector(2^n,[k$k=1..2^n]);

d_real := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16])

虚部には 0 を入れておきます。

>    d_imag:=vector(2^n,0);

d_imag := vector([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

listplot を使ってデータをプロットしてみましょう。

>    listplot(d_real);

[Maple Plot]

FFTを呼び出して計算してみましょう。

>    FFT(n,d_real,d_imag);

16

データ数が出力されます。結果のデータをプリントしてみましょう。

>    print(d_real,d_imag);

vector([136, -7.999999980, -8.000000003, -7.999999999, -8., -8.000000001, -7.999999997, -7.999999999, -8, -8.000000002, -7.999999997, -7.999999999, -8., -8.000000001, -8.000000003, -8.000000019]), vect...
vector([136, -7.999999980, -8.000000003, -7.999999999, -8., -8.000000001, -7.999999997, -7.999999999, -8, -8.000000002, -7.999999997, -7.999999999, -8., -8.000000001, -8.000000003, -8.000000019]), vect...
vector([136, -7.999999980, -8.000000003, -7.999999999, -8., -8.000000001, -7.999999997, -7.999999999, -8, -8.000000002, -7.999999997, -7.999999999, -8., -8.000000001, -8.000000003, -8.000000019]), vect...

プロットしてみましょう。

>    p1:=listplot(d_real):
p2:=listplot(d_imag):
display(array(1..2,[p1,p2]));

[Maple Plot]

処理の例

データ数を増やして (n=10) もっと実際的な計算を行ってみましょう。

観測時間は 1 秒間としましょう。

>    n:=12:
2^n;

4096

>    d_r:=vector(2^n,0):d_i:=vector(2^n,0):

>    k:=0:
for t from 0 by 1/(2^n-1) to 1 do
  k:=k+1:
  d_r[k]:=evalf(sin(2*Pi*442*t)):
end do:

>    listplot([seq(d_r[k],k=1..50)]);

[Maple Plot]

これは 442 ヘルツの波をプロットしたものです。サンプリング間隔によって山の部分が滑らかではなくなっています。
高周波になればなるほどデータ数を多くする必要があることが想像できます。FFT を計算しましょう。

>    FFT(n,d_r,d_i);

4096

結果の実部をlistplotでプロットしてみましょう。

>    listplot(d_r);

[Maple Plot]

ピークが 2 個所でています。FFT のスケールはデータの次元に正規化されていないことに注意してください。
ピーク付近を拡大してプロットしてみましょう。

>    p1:=listplot([seq(d_r[k],k=430..450)]):
p2:=listplot([seq(d_r[k],k=3630..3700)]):
display(array(1..2,[p1,p2]));

[Maple Plot]

同様に虚部をプロットしてみましょう。

>    p1:=listplot([seq(d_i[k],k=430..450)]):
p2:=listplot([seq(d_i[k],k=3630..3700)]):
display(array(1..2,[p1,p2]));

[Maple Plot]

実部と虚部の平方和の平方根がパワーになります。スケールは正規化されていないことに注意してください。

>    pow:=vector(2^n,0):

>    for k from 1 to 2^n do
  pow[k]:=evalf(sqrt(d_r[k]^2+d_i[k]^2));
od:

>    listplot(pow);

[Maple Plot]

ピーク付近を拡大してみましょう。

>    listplot([seq(pow[k],k=3630..3700)]);

[Maple Plot]

例2(ノイズ)

Maple では乱数を発生させる関数 rand() が用意されています。
これは 12 桁の乱数を発生しますので次のように -0.5 から 0.5 の間の乱数を発生することができます。

>    n:=12:
d_r:=vector(2^n,0):
d_i:=vector(2^n,0):
k:=0:for t from 0 by 1/(2^n-1) to 1 do
  k:=k+1:
  d_r[k]:=evalf(rand()/10^12-0.5):
end do:
listplot([seq(d_r[k],k=1..100)]);

[Maple Plot]

これをノイズとして先ほどのデータに足し込んでみましょう。

>    k:=0:for t from 0 by 1/(2^n-1) to 1 do
  k:=k+1:
  d_r[k]:=evalf(sin(2*Pi*442*t)+(rand()/10^12)-0.5):
end do:

>    listplot([seq(d_r[k],k=1..100)]);

[Maple Plot]

これを FFT 計算してみましょう。

>    FFT(n,d_r,d_i):
listplot(d_r);

[Maple Plot]

全帯域にノイズが広がっているのが解ります。

今度は周波数成分 (442 ヘルツ) にノイズを加えてみましょう。周波数が 5 ヘルツの帯域分広がります。

>    d_r:=vector(2^n,0):
d_i:=vector(2^n,0):

>    k:=0:
for t from 0 by 1/(2^n-1) to 1 do
  k:=k+1:
  d_r[k]:=evalf(sin(2*Pi*(442+5*((rand()/10^12)-0.5))*t)):
end do:

>    listplot([seq(d_r[k],k=1..100)]);

[Maple Plot]

これをFFT計算してみます

>    FFT(n,d_r,d_i):
listplot(d_r);
pow:=vector(2^n,0):

[Maple Plot]

同様にパワーをプロットしてみます。

>    for k from 1 to 2^n do
  pow[k]:=evalf(sqrt(d_r[k]^2+d_i[k]^2));
end do:

>    listplot(pow);

[Maple Plot]

全体にノイズが大きくなっているように感じますが、ピークも落ちていることに注意してください。

>    p1:=listplot([seq(d_i[k],k=435..450)]):
p2:=listplot([seq(d_i[k],k=3640..3675)]):
display(array(1..2,[p1,p2]));

[Maple Plot]

帯域

計算時間を少なくするために n を 11 にします。周波数 430 ヘルツから 450 ヘルツまで 1 ヘルツ間隔で足し込んでみましょう。

>    n:=11;
d_r:=vector(2^n,0):d_i:=vector(2^n,0):
k:=0:
for t from 0 by 1/(2^n-1) to 1 do
  k:=k+1:
  d_r[k]:=evalf(sum(sin(2*Pi*w*t),w=430..450)):
end do:

n := 11

>    listplot([seq(d_r[k],k=1..200)]);

[Maple Plot]

>    FFT(n,d_r,d_i);

2048

パワーを計算してみましょう。

>    pow:=vector(2^n):

>    for k from 1 to 2^n do
  pow[k]:=evalf(sqrt(d_r[k]^2+d_i[k]^2));
end do:

>    listplot(pow);

[Maple Plot]

>    listplot([seq(pow[k],k=300..700)]);

[Maple Plot]

バンド幅(帯域)が確認できます。

ソースプログラム

Maple では FFT のプログラムのソースを参照できます。以下のようにすればリストされます。

>    interface(verboseproc=2):

>    print(FFT);

proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...
proc (m, x, y) local n, l, le, le1, ur, ui, wr, wi, j, i, ip, tr, ti, sr, si, nv2, nm1, k; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; n := 2^m; le1 := n; for l to ...