周期2πの関数 f(t) = |t| は次のようにすればできます.
octave:1> t = -3*pi: 0.01: 3*pi; octave:2> x = 2*pi*abs(t./(2*pi) - round(t./(2*pi) ) ); octave:3> plot(t,x);
しかけは以下のとおり:
octave:4> plot(t, t, t, round(t) );round(t) は四捨五入して整数値を求める関数. t と round(t) の差は常に -0.5 と 0.5 の間にある. t が0 から走ると, 0.5 までは差が t だけあり,単調に増加する. 0.5 から 1 までは差(絶対値)が 1-t で, 単調に減少する. 1 を越えると,また初めから同じことの繰り返しになる. つまり,この場合,周期は 1.
octave:5> plot(t, (t - round(t)) );これがノコギリ波. これに絶対値をとると,三角波になる.
octave:6> plot(t, abs(t - round(t)) );
上の波の周期Tをみてみよう. ノコギリ波も三角波も周期 T=1 となっている. 例えば周期を 2π にしたい場合は,どうしたらいいだろうか. 2倍にしたい分,互いをゆっくり変化させればよいので,
octave:7> plot(t, abs(t./(2*pi) - round(t./(2*pi)) ) );とする. t./(2*pi) というのは,ベクトル t の各要素を 2π で割る操作のこと.
octave:8> plot(t, abs( 1/2/pi*t - round( 1/2/pi * t )) );これと同じこと.t はベクトルで,それ以外は単なる数である ことに注意したい.
ノコギリ波は±0.5の間を,三角波は0と0.5の間に波がある. つまりノコギリ波の振幅は 0.5, 三角波の振幅は 0.25. この振幅の大きさを変えるのは,全体を何倍かすればよい. 例えば三角波の振れ幅を1にするには
octave:9> plot(t, 2*abs(t - round(t)) );とすればよい.
波を作るには,縦にどれだけずれるか調整しておく必要がある. これは
octave:10> plot(t, 2*abs(t - round(t)) .+0.5 );で,縦に0.5シフトする.
横軸も調整しておく必要があれば
octave:11> plot(t, 2*abs( (t-0.1) - round( (t-0.1) ) .+0.5 ) );で,横に0.1シフトする.
例の[-π: t : π] の区間で f(t) = |t| と定義されている周期2πの関数を描くには, 縦軸は 0 からπの間で変化し, f(0)=0 であることに注意しておく.
octave:12> t = -3*pi: 0.01: 3*pi; octave:13> x = 2*pi*abs(t./(2*pi) - round(t./(2*pi) ) ); octave:14> plot(t,x);これでできた. フーリエ級数と比較してみよう.
octave:15> y0 = pi/2*cos(0*t); # 定数ベクトルをつくった octave:16> y2 = pi/2 - 4/pi*( cos(t) + 1/9*cos(3*t) ); octave:17> plot (t,x, ";f(t)=|t|;", t, y2, ";fourier2;");図に複数の線を書くと,どれがどれだかわからなくなる. こうやって(;ではさむ;)曲線の意味をグラフに書くことができる.
f(t)=|t| と y2 で,各要素の差が どのくらいあるかを見るには
octave:18> plot(t,x-y2,';error;')とすればよい.
octave:19> y3 = pi/2 - 4/pi*( cos(t) + 1/9*cos(3*t) + 1/25*cos(5*t) ); octave:20> plot(t,x-y2,';error2;', t, x-y3,';error3;' )
octave:21> y4 = pi/2 - 4/pi*( cos(t) + 1/9*cos(3*t) + 1/25*cos(5*t) + 1/49*cos(7*t) ); octave:22> plot(t,x-y2,';error2;', t, x-y3,';error3;', t, x-y4,';error4;' )
octave:23> n2 = sqrt ( (x-y2)*(x-y2)') n2 = 1.8880' は転置です. (x-y2) などは横ベクトルです. 距離を求めるには 自分自身との内積のルートを取ればよかったのですから, sqrt( (x-y2) * (x-y2)' ) で計算できます. ベクトル演算なので for 文を書かなくてもできます.
octave:24> n3 = sqrt ( (x-y3)*(x-y3)') n3 = 1.0580 octave:25> n4 = sqrt ( (x-y4)*(x-y4)') n4 = 0.69488確かに,誤差が少なくなっていっています. 正確にはこれらを定数(2πなど)で割る必要があると思うが, 距離の大きさだけを比較するので,ここでは省略.
octave:26> T = 2*pi; % 周期 octave:27> A = 0.7; % 振幅 octave:28> s=abs(t./T-floor(t./T+0.5)); octave:29> x=0.5.*A.*(sign(0.25-s)+1.0); octave:30> axis ([-10,10, -0.3, 1]); octave:31> plot(t,x)
octave:32> (-1)^7 ans = -1y2 ベクトルの第2番目の成分の値
octave:33> y2(2) ans = 2.9854y2 ベクトルのサイズ
octave:34> size(y2) ans = 1 1885これは 1行 1885 列の行列, つまり 1885次元ベクトルということです. 中身は y2 と入力すれば(;を付けなければ)全部表示されます.