B11 IRの切り出しと窓関数

室長:今回はインパルスレスポンスの切り出しをやってみよう。44kで測定したデータ’test.wav’を使って適当にやってみてくれ。

助手:は~い。いつも通りoctaveを立ち上げて…


a=wavread ("test.wav");
a=a(:,1);
b=pcmread("inverse44k2sec.pcm");
c=fftconv(a,b);
plot(c);


反射の前を切り出せば良いので

c2=c(132900:133176);
filterplot(c2,44100);

切り出したサンプル数が276なのですがインパルス以降のデータは211サンプルなので44100/211≒210Hz以下は信頼できないけど高域は部屋の影響を除いたデータのはずです。

室長:まあ大体そんな感じだな。だが正確にいうと周期性のない長いデータから切り出した時には端の影響を避けるために窓関数を使った方が良い。とくにノイズなどの影響で端の部分の値が0から離れている場合にはな。素人測定だからあまり気にしなくても良いが一応窓関数を使ってみよう。適当な窓関数をでっち上げてやってみてくれ。インパルスを中心に非対称にしてくれ。

助手:分かりました。

mado.m

#fuction a=mado(n,t)
#切り出したインパルスレスポンス用の窓関数
#長さ:n、インパルスの位置:t
function a=mado(n,t)
  a1=[0:-1:-t+1];
  a2=[t-n+1:0];
  a1=1-exp(a1/t*6);
  a2=1-exp(a2/(n-t)*6);
  a=[a1 a2];
  a=a';
endfunction

 

plot(mado(277,69));

こんなんでどうですか?少ないサンプル数で周波数分解能を確保するためにかなり無理やりなカーブになってます。

ついでだからインパルスのピークを自動判定する関数も作って見ました。

mado2.m


#function b=mado2(a);
#切り出したインパルスレスポンス:a
#に窓関数madoをかける
function b = mado2 (a)
[m1 m2]=max(abs(a));
b = a .* mado(length(a),m2);
endfunction

実際に試してみます。


c3=mado2(c2);
clf;
hold on;
filterplot(c2,44100);
filterplot(c3,44100,'1');


この窓関数だと低域と高域で少し違うだけであまり差がありませんね。

室長:まあ切り出したデータに窓関数を掛けたほうがなんとなく格好良いではないか。では次回は測定データから周波数特性の補正フィルタを作るぞ。

2013年1月7日

次ページへ

目次へ