GNU Octave 3.8.1のバグ

Octave 3.8.1にはバグがあってfir1が正しく動きません。 Ubuntuの場合、14.04、14.10が該当します。

Ubuntu 15.04でインストールされるOctave 3.8.2は大丈夫です。

Raspbianでインストールされる3.6.Xも大丈夫です。

fir1のバグ

signal packageに入っているfir1には長いこと報告されていたのに修正されていなかったバグがありました。

fir1の中からfir2関数を呼び出すところでgridサイズが512に決め打ちされていたのです。

そのためtap数の多いフィルタを設計しようとするとエラーが起きてしまいます。これも3.8.2でようやく修正されました。

fir1.mを見てみるとfir2を呼び出すところでgridサイズを[]としてfir2のデフォルト値を使うようになりました。これでtap数の大きなフィルタを計算してもエラーが出ないようになりました。

しかし、これでも不十分です。 gridサイズが小さいのでフィルタの精度が出ないのです。

3.8.2でクロス周波数300Hz、3000Hzの3wayのフィルタを設計してみます。

fig1_1

この様に一見正常な特性に見えます。しかし拡大してみると300Hzになるはずのクロスが 288Hzになっています。

fig1_2

同様に3000Hzも2950Hzになってます。

fir1.mを修正してしまいましょう。

~/.octaverc を編集します。なければ作成します。

format long;
addpath ("home/hoge/octave_function"); #ディレクトリは環境に合わせて修正する

これでoctaveを起動すると/home/hoge/octave_functionディレクトリ内の関数ファイルが優先して使われるようになります。 このディレクトリにfir1.mをコピーして修正します。

cd octave_function cp /usr/share/octave/packages/signal-1.3.0/fir1.m ./ #octave 3.8.2の場合

適当なエディタでfir2を呼び出しているところを編集します。 3.8.2の場合には117行目のグリッドサイズをどーーんと大きくします。

b = fir2(n, f, m, [], 2, window);

修正後

b = fir2(n, f, m, 8192*32, 2, window);

これで上書き保存します。

再びフィルタを設計してみます。 赤がgridを大きくして設計したものです。

fig2_1

拡大してみます。

fig2_2

こんどは正しくクロスが300Hzになってます。

あと自作の3wayと4way用の関数ファイルを貼っておきます。これをパスの通ったディレクトリにコピーすればfir3wayとfir4way関数が使えるようになります。

fir3way.m

## Copyright (C) 2015 Amanogawa Audio Labo
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

# function [low, mid, high]=fir3way(fs,f1,f2,n1,n2[,w])
# 3way用FIRフィルタの設計
# n : 遮断周波数の何倍のtapを使うか?
# w : 窓関数 @hammingなど 省略時はkaiser,beta=9
function [low, mid, high] = fir3way (fs,f1,f2,n1,n2,w)

  if !(exist("w"))
    w=@(x)kaiser(x,9);
  endif


  if (f1 > f2)
    tmp=f1;
    f1=f2;
    f2=tmp;
  endif
  t1=floor(fs/f1*n1);
  t2=floor(fs/f2*n2);
  t1=t1+mod(t1,2);
  t2=t2+mod(t2,2);
  low=fir1(t1,f1/fs*2,'low',w(t1+1));
  m1=fir1(t1,f1/fs*2,'high',w(t1+1));
  m2=fir1(t2,f2/fs*2,'low',w(t2+1));
  high=fir1(t2,f2/fs*2,'high',w(t2+1));
  mid=fftconv(m1,m2);
  if ( size(low)(1) > 1 )
    low=low';
    mid=mid';
    high=high';
  endif
  lm=length(mid);
  ll=length(low);
  lh=length(high);
  lz=zeros(1,(lm-ll)/2);
  hz=zeros(1,(lm-lh)/2);
  low=[lz,low,lz];
  high=[hz,high,hz];
endfunction

fir4way.m

## Copyright (C) 2015 Amanogawa Audio Labo
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

# function [l, ml,mh, h]=fir4way(fs,f1,f2,f3,n1,n2,n3[,w])
# 4way用FIRフィルタの設計
# n : 遮断周波数の何倍のtapを使うか?
# w : 窓関数 @hammingなど 省略時はkaiser,beta=9
#
function [l, ml, mh, h] = fir4way (fs,f1,f2,f3,n1,n2,n3,w)

  if !(exist("w"))
    w=@(x)kaiser(x,9);
  endif

  t1=floor(fs/f1*n1);
  t2=floor(fs/f2*n2);
  t3=floor(fs/f3*n3);
  t1=t1+mod(t1,2);
  t2=t2+mod(t2,2);
  t3=t3+mod(t3,2);
  l=fir1(t1,f1/fs*2,'low',w(t1+1));
  ml1=fir1(t1,f1/fs*2,'high',w(t1+1));
  ml2=fir1(t2,f2/fs*2,'low',w(t2+1));
  mh1=fir1(t2,f2/fs*2,'high',w(t2+1));
  mh2=fir1(t3,f3/fs*2,'low',w(t3+1));
  h=fir1(t3,f3/fs*2,'high',w(t3+1));
  ml=fftconv(ml1,ml2);
  mh=fftconv(mh1,mh2);
  if ( size(l)(1) > 1 )
    l=l';
    ml=ml';
    mh=mh';
    h=h';
  endif
  ll=length(l);
  lml=length(ml);
  lmh=length(mh);
  lh=length(h);
  m1=max([ll,lml,lmh,lh]);

  lz=zeros(1,(m1-ll)/2);
  mlz=zeros(1,(m1-lml)/2);
  mhz=zeros(1,(m1-lmh)/2);
  hz=zeros(1,(m1-lh)/2);
  l=[lz,l,lz];
  ml=[mlz,ml,mlz];
  mh=[mhz,mh,mhz];
  h=[hz,h,hz];
endfunction
2015年8月15日