Octave では FFT が fft() 関数一発で実行できます。次の実行例でその雰囲気がわかります。

フーリエ変換と逆フーリェ変換

x = 0 から x = 5 までの区間の 256 等分点をサンプル点として y = exp(-x) のサンプリングデータ a を作成します。

octave:1> t = linspace (0, 5, 256);
octave:2> a = exp (-t);
octave:3> plot (a)

それをフーリエ変換して b に入れます。fft() の出力は複素数なので絶対値を計算する abs() 関数を利用して、plot(abs(b)) で振幅スペクトルを表示します。(厳密にはパワースペクトラム密度(PSD)の計算が必要です。)

octave:4> b = fft (a)
octave:5> plot ( abs(b) )

b を ifft() で逆フーリエ変換して c に入れます。ifft の出力は複素数ですが、虚数部は0なので実数部のみ real() で取りだします。

octave:6> c = real(ifft(b));
octave:7> plot (c) 

雑音信号の周波数分析

5 Hz と 50Hz の信号成分を含む雑音信号を作成してフーリエ変換すると信号成分の周波数を抽出することができます。b(1:128) として範囲を指定すると振幅スペクトルの折り返しの部分を表示しないようにできます。

octave:1> t = linspace (0, 1, 256);
octave:2> a = sin(10*pi*t) + 0.7*cos(100*pi*t) + randn(size(t));
octave:3> plot (a)
octave:4> b = fft(a);
octave:5> plot ( abs(b(1:128)) )

フィルタリング

フーリエ変換の応用で、雑音にまぎれた信号を取りだすことができます。まず元の信号 a を作成します。

octave:1> t = linspace (0, 1, 256);
octave:2> a = sin(4*pi*t) + cos(6*pi*t);
octave:3> plot (a)

次に a に雑音を混入したサンプル b を作成します。

octave:4> b = a + randn(size(t));
octave:5> plot (b)

b をフーリエ変換したデータ c を作成します。

octave:6> c = fft (b);
octave:7> plot (abs(c(1:128)) 

c のうち雑音成分をカットします

octave:8> for i = (1:256) 
> if abs(c(i)) < 50
> c(i) = 0
> endif
> endfor
octave:9> plot (abs(c(1:128)))

c を逆フーリエ変換してもとのデータを取りだします。

octave:10> d = real(ifft(c));
octave:11> plot (a, "-", d, "-")

白色雑音の入出力分析による入出力関数の推定

入出力関数の不明な装置に白色雑音を入力し、出力とあわせて周波数分析をすることによって入出力関数を推定することができます。そのシミュレーションをやってみましょう。

まず指数関数的に減衰するインパルス応答を作成します。

octave:> t = linspace (0, 100, 1000);
octave:> a = exp(-0.1*t);
octave:> plot (a)

次に入力となる白色雑音を作成します。

octave:> b = randn(size(t));
octave:> plot (b)

b を装置に入力してその出力を取りだします。シミュレーションでは a と b の convolution を計算します。

octave:> c = conv (a, b);
octave:> plot (c)

c のうち 入力 b の直接の出力の部分 (1:1000) の値を取りだします。

octave:> c1 = c(1:1000);
octave:> plot (b, "-", c1, "-")

b と c1 をフーリエ変換して c1 と b の商 df を計算します。これが求めるインパルス応答関数の周波数特性になります。そして、df を逆フーリエ変換した d が求めるインパルス応答関数の推定値になります。乱数が完全な白色雑音ではないので多少でこぼこしていますが、グラフを見ると元の入出力関数の推定値になっていることが分かります。

octave:> bf = fft(b);
octave:> cf = fft(c1);
octave:> df = cf ./ bf;
octave:> d = real(ifft(df));
octave:> plot (d, "-", a, "-")