B13 f特の補正2

室長:では周波数特性の逆特性を作るぞ。補正関数を作ってくれ。

助手:ほい。

fhosei.m

# hc = fhosei (h,y1,y2)
#y1補正する上限値
#y2補正する目標値
#
function hc=fhosei(h,y1,y2)
  hdb=a2db(abs(h));
  hc=y2-hdb;
  for n=1:length(h)
    if (hdb(n) > y1)
      hc(n)=y2-y1;
    endif
    if (hdb(n) < y2)
      hc(n)=0;
    endif
  endfor
  hc=db2a(hc);
endfunction

動作確認してみます。切り出したインパルスレスポンスをc2とします。

c3=mado2(c2);
[h w]=freqz(c3,1,4096,44100);
hs=fsmooth(h,w,3);
hs2=fsmooth(hs,w,6);
hc=fhosei(hs2,0,-18);
clf;hold on;
semilogx(w,a2db(hs2),'1');
semilogx(w,a2db(hc));
figure;
semilogx(w,a2db(hs2)+a2db(hc));


うまくいきました。ピークを下げて-18dBでフラットにするカーブができました。

室長:次はいま計算したhcからfir2関数でフィルタを計算する。

助手:ふぁーい。

fir2m.m


#a = fir2m(n,h,w,fs)
function a = fir2m (n,h,w,fs)
  f=w/(fs/2);
  f=[f;1];
  m=[abs(h);abs(h(length(h)))];
  a=fir2(n,f,m,4096*4,kaiser(n+1,8));
endfunction

この関数も確認してみます。

a=fir2m(254,hc,w,44100);
filterplot(a,44100);
figure;
filterplot(fftconv(c2,a),44100);


うまくフラットに補正できました。

室長:では次回はフィルタの種類を変えてみよう。

2013年1月12日

次ページへ

目次へ