A07 Octaveの関数ファイルを作る

室長：octaveの基本文法は他のサイトをあたってくれ。とりあえず前回の解説をするぞ。最初の行はaに行ベクトルを格納している。セミコロン(;)は出力を抑制するときにつける。つぎの行はfreqzという関数でフィルタを計算している。freqzの説明はoctaveのコマンドラインで

help freqz

とすれば表示される。freqzが定義されているファイルの場所も表示されるはずじゃ。必要とあればfreqz.mの内容をみてどういう計算がされているのか知ることも出来る。最後の3行目だがこれはsemilogxという関数でX軸を対数表示にしてグラフを描いている。普通のグラフはplot、y軸が対数の場合にはsemilogyじゃ。

20*log10()というのは何をしているのかわかるな？

助手：それくらいは私でも解ります。dBに変換しているんですよ。n倍をdBにするには
20*log10(n)
n dBが何倍かの計算は
10^(n/20)
です。

室長：それではdB計算を関数にしてみてくれ。

助手：えーと、functionとendfunctionで挟めばいいんでしたね。

function db = a2db (a)
  db=20*log10(a);
endfunction

function a = db2a (db)
  a=10.^(db/20);
endfunction

こんな感じでどうですか。aやdbはローカル変数となって他には影響しないんですよね。

室長：よいだろう。ただ毎回これを入力するのは手間なので関数ファイルというのをつくる。1個の関数をコメントをつけて１つのファイルにする。”ファイル名は間数名.m”とする。ファイルの場所はパスが通っている場所にする。カレントディレクトリならOKじゃ。べつのディレクトリにするばあいには ~/.octaverc ファイルを作ってaddpath (“パスをとおすディレクトリの絶対パス”);とでもしとけば良い。

a2db.m

#function db = a2db (a)
#増幅率をdBに変換
#
function db = a2db (a)
  db=20*log10(a);
endfunction

db2a.m

#function a = db2a (db)
#dBを増幅率に変換する
#
function a = db2a (db)
  a=10.^(db/20);
endfunction

この2つのファイルを作っておけばいちいち入力しなくても最初から組み込まれている関数のように使えるのだ。他の関数から呼び出すことも可能じゃ。help a2dbとすればコメントが表示されるぞ。

ではフィルタを表示するfilterplotという関数を作ってみてくれ。コメントも忘れずにな。今作った関数a2dbもそのまま使用できるぞ。

助手：出来ました。

filterplot.m

#function filterplot (a,fs[,n][,fmt])
#フィルタaのグラフをdB表示、x軸対数で表示
#nを省略すれば4096

function filterplot (a,fs,n,fmt)

  n1=4096;
  fmt1='3';

  if (exist("n") && ischar(n))
    fmt1=n;
  elseif (exist("n"))
    n1=n;
    if (exist("fmt"))
      fmt1=fmt;
    endif
  endif

  [h,w]=freqz(a,1,n1,fs);
  semilogx(w,a2db(abs(h)),fmt1);
  grid on;

endfunction

室長：センスのない関数だな。まあ動けばよい。これでfilterplot(a,44100);と入力するだけでグラフが表示できる。この関数ファイルはこれからも使うので作って保存しておいてくれ。次回はoctave-signalにあるfir1という関数を使ってフィルタを設計してみよう。

2012年12月11日