フィルターバンク / Reconstructing Filterbank

帯域ごとに分割し、また元の信号に戻せるようなフィルターバンクが欲しいときがあります。JIS規格にのっとったフィルターではないので他者の研究と比べるときには不都合があるかもしれませんが、グラフィック・イコライザなど様々な用途に使えそうです。

原信号とフィルターバンクから再合成した信号との誤差の大きさを調べるために、以下のような計算をしてみました。

fs = 44100;
x = randn(fs*10, 1);
z = filterbank(x, fs);
zz = sum(z, 2);
err = x - zz;
10*log10(sqrt(mean(err.^2))/sqrt(mean(x.^2)))
corr(x, zz)

10秒の白色雑音を生成し、再合成した信号との誤差の大きさをデシベル値で計算します（筆者環境では誤差は25〜35dB程度）。また、最後の行では原信号と再合成信号とのピアソン相関係数を計算します（筆者環境では1.0000）。

function [out, cfreqs] = filterbank(x, fs, oct)
%FILTERBANK   octave filterbank using raised cosine frequency window
%   [out, cfreqs] = filterbank(x, fs) create filtered signals `out'
%   from x.  Columns of `out' is a filtered signal with
%   corresponding center frequency in `cfreqs'.  By default,
%   1/1-octave band filters are used.
%
%   Filtered signals at the lowest frequency is processed with
%   low-pass, the highest frequency is processed with high-pass,
%   and all frequencies in between are processed with band-pass
%   filters.  This yields to (almost) perfect reconstruction of the
%   original signal, such that:
%      out = filterbank(x, fs);
%      z = sum(out, 2);
%   creates z almost the same as x.
%
%   [cfreqs, out] = filterbank(x, fs, oct) create and filters using
%   1/oct-octave band filters.
%
%   2011-07-13 MARUI Atsushi

if nargin<3
  oct = 1;
end
cfreqs = 1000 * 2.^((-15:12)/oct);
cfreqs = cfreqs(20<=cfreqs & cfreqs<=fs/2);



%% create octave-band bandpass filter
out = zeros(length(x), length(cfreqs));
for n = 1:length(cfreqs)
  fc = cfreqs(n);
  fl = fc * 2.0^(-1/oct);
  fh = fc * 2.0^(+1/oct);
  
  W = zeros(length(x)/2, 1);
  f = linspace(0, fs/2, length(W))';
  ff = fl<=f & f<=fh;
  cc = (cos(log2(f(ff)/fc)*pi*oct) + 1) / 2;
  W(ff) = cc;
  if n==1
    ff = f < fc;
    W(ff) = 1;
  elseif n==length(cfreqs)
    ff = fc < f;
    W(ff) = 1;
  end
  W = [W ; flipud(W)];
  W(1) = 0;
  W(end/2+1) = 0;
  
  % apply filter
  y = real(ifft(fft(x) .* W));
  out(:,n) = y;
end