複素フーリエ係数 C_k とフーリエ変換 F(ω) は,性質が類似している. 今の段階では混同してもよいので,まず第一に,複素フーリエ級数, フーリエ変換の雰囲気をつかむ. (手計算では連続関数を扱うが, コンピュータで扱うのは離散データであるので, 以下に示す octave のコマンドでは離散フーリエ変換をおこなっている.)
#より後ろはコメントなので入力する必要はない
octave:1> t = linspace (0, 3, 256); # 横軸の作成. 0から3の区間を256等分 octave:2> a = sin(2*pi*t); # 基本周期1 の正弦波 を作る octave:3> plot (t, a, ";sin2πt;")
octave:4> b = fft (a); # fft は fast fourier transform(高速フーリエ変換)の略)
octave:5> size(b)
ans =
1 256
これは b ベクトルが 1 x 256 行列,
つまり 256次元の縦ベクトルであることを示している.
octave:6> b
b =
Columns 1 through 3:
-0.00000 + 0.00000i 0.00458 - 0.37298i 0.02914 - 1.18712i
Columns 4 through 6:
4.70103 - 127.63366i -0.08530 + 1.73639i -0.05798 + 0.94374i
................................ (q を押すと表示終了)
複素数なので,グラフを描くときには注意が必要.
octave:7> w = linspace (1, 256, 256); # 横軸の作成. 1から256の区間を256等分 octave:8> plot ( w, abs(b) )横軸の番号はフーリエ係数 Ck の kとは違うが,今のところは 気にしなくてよい.
octave:9> a = sin(2*pi*t) + sin(6*pi*t) ; octave:10> plot (t, a, ";sin2πt+sin6πt;") octave:11> b = fft (a); octave:12> plot ( w, abs(b) ) # 振幅スペクトル octave:13> plot ( w, angle(b) ) # 位相スペクトル振幅スペクトルは左右対称で, 二つのピークがでることを確認する. 位相スペクトルは,左右で符合が反対になっていることを確認する.
octave:14> plot ( w, real(b) ) # (フーリエ係数) 実部 octave:15> plot ( w, imag(b) ) # 虚部
以下おまけ.
データの個数は256個. したがってそれをフーリエ変換した結果も 256個のデータからなる. 上の図では横軸を k としているが, 縦軸の値は C_k であるが,C_(-k) はC_{256-k} と表示されている と考えればよい.
以下のようなファイルをエディターで作っておけば再実験が楽にできる.
% octave foo301.octfoo301.oct
clear
t = linspace (0, 3, 256);
a = sin(2*pi*t) + sin(6*pi*t) ; # 基本周期1 の正弦波 を作る
grid "on"
title("2つの正弦波の重ね合わせ波形")
xlabel("sin2πt+sin6π")
plot(t,a)
input ("Press Return");
b = fft (a);
w = linspace (0, 255, 256);
subplot(2,1,1) # 縦2段, 横1段で図を配置。 1番目の図
plot(w, abs(b) )
xlabel("k")
ylabel("|Ck|")
title("Power Spectrum")
subplot(2,1,2)
plot(w, angle(b) )
title("Phase Spectrum")
xlabel("k")
ylabel("arg Ck")
input ("Press Return");
octave:16> a = sin(2*pi*t) + sin(6*pi*t) + 0.2*randn(size(t));先の結果と比べて,いろいろな周波数成分が含まれていることを確認できる.
octave:17> a = sin(2*pi*t) + sin(6*pi*t) ; octave:18> plot (t, a, ";sin2πt+sin6πt;") octave:19> b = fft (a); octave:20> plot ( abs(b) ) octave:21> b = b + 2*randn(size(b)); octave:22> plot ( abs(b) ) octave:23> c = ifft(b); octave:24> plot (t,c ) octave:25> b = b + 50*randn(size(b)); octave:26> plot ( abs(b) ) octave:27> c = ifft(b); octave:28> plot (t,c )