#より後ろはコメントなので入力する必要はない
octave:1> t = linspace (0, 5, 256); # 横軸の作成. 0から5の区間を256等分.tは縦ベクトル. octave:2> f = exp (-t); octave:3> g = shift(f, 50); # 横に 50だけずらす操作. octave:4> plot (t, f, ";ext(-t);", t,g, ";ext(-(t-50));") # ";comment;"でグラフに名前を付けられる.
octave:5> F = fft (f); # fft は fast fourier transform(高速フーリエ変換)の略) octave:6> G = fft (g);F,G は複素数値をもつベクトルなので,グラフを描くときには注意が必要.
F = Columns 1 through 3: 51.16136 + 0.00000i 20.23596 - 24.70413i 7.46969 - 17.44847i Columns 4 through 6: 3.85317 - 12.59200i 2.44173 - 9.72234i 1.75952 - 7.88322i ............................... (q を押すと表示終了)256次元のベクトルをフーリエ変換したので,その結果は256次元ベクトルに なっている. 要素が複素数なので,2次元のグラフに描くには, 複素数値を絶対値と偏角に分けて表示するか, 実数部と虚数部に分けて表示する方法がある. 絶対値と偏角に分けて表示する場合:
octave:7> w = linspace (1, 256, 256); # 横軸の作成. 1から256の区間を256等分 octave:8> plot ( w, abs(F), ";F;", w, abs(G), ";G;") # 二つのグラフが同じであると一つしか見えない... octave:9> plot ( w, angle(F), ";F;", w, angle(G), ";G;")横軸の番号は大雑把にいうと周波数に対応する.
octave:10> plot ( w, real(F), ";F;", w, real(G), ";G;" ) octave:11> plot ( w, imag(F), ";F;", w, imag(G), ";G;")
octave: > g = shift(f, 1); # 横に 1だけずらす操作. octave: > G = fft (g); octave: > plot ( w, imag(F), ";F;", w, imag(G), ";G;") octave: > g = shift(f, 2); # 横に 2だけずらす操作. octave: > G = fft (g); octave: > plot ( w, imag(F), ";F;", w, imag(G), ";G;") .....
% octave foo401.octfoo401.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番目の図
[xb, yb] = bar (w, abs(b));
plot(xb, yb)
xlabel("k")
ylabel("|Ck|")
title("Power Spectrum")
subplot(2,1,2)
[xb, yb] = bar (w, angle(b));
plot(xb, yb)
xlabel("k")
ylabel("arg Ck")
title("Phase Spectrum")
input ("Press Return");
a0 = shift(a, 10); # 横に 10だけずらす操作.
b0 = fft (a0);
subplot(2,1,1) # 縦2段, 横1段で図を配置。 1番目の図
[xb, yb] = bar (w, abs(b0));
plot(xb, yb)
xlabel("k")
ylabel("|Ck|")
title("Power Spectrum (shift=10)")
subplot(2,1,2)
[xb, yb] = bar (w, angle(b0));
plot(xb, yb)
xlabel("k")
ylabel("arg Ck")
title("Phase Spectrum (shift=10)")
input ("Press Return");
a0 = shift(a, 20); # 横に 20だけずらす操作.
b0 = fft (a0);
subplot(2,1,1) # 縦2段, 横1段で図を配置。 1番目の図
[xb, yb] = bar (w, abs(b0));
plot(xb, yb)
xlabel("k")
ylabel("|Ck|")
title("Power Spectrum (shift=20)")
subplot(2,1,2)
[xb, yb] = bar (w, angle(b0));
plot(xb, yb)
xlabel("k")
ylabel("arg Ck")
title("Phase Spectrum (shift=20)")
input ("Press Return");
octave:12> t = -3*pi: 0.01: 3*pi;
octave:13> length = size(t,2); # 後々のため,t が何次元のベクトルか計算しておく.この場合は1895次元になっている.
octave:14> T = 2*pi; % 周期
octave:15> A = 0.7; % 振幅
octave:16> t0 = 50; % シフト(横ずらし)する量
octave:17> s=abs(t./T-floor(t./T+0.5));
octave:18> f=0.5.*A.*(sign(0.25-s)+1.0);
octave:19> axis ([-10,10, -0.3, 1]);
octave:20> plot(t,f)
octave:21> g = shift(f, t0);
octave:22> F = fft(f);
octave:23> G = fft(g);
octave:24> w = linspace (1, length, length);
octave:25> axis("auto")
octave:26> plot ( w, abs(F), ";F;", w, abs(G), ";G;")
octave:27> plot ( w, angle(F), ";F;", w, angle(G), ";G;")
octave:24> axis("auto")
としておく.