Previous: 音生成 / Synthesizers Up: Matlab M-Files for Audio Next: 周波数応答 / Frequency Response
function [B,A] = biquad_lpf(fc, q, fs) %BIQUAD_LPF biquad low-pass filter % [B,A] = biquad_lpf(fc, q, fs) % % Input: % fc - cutoff frequency (Hz) % q - Q-factor % fs - sampling frequency (Hz) % % 2005-02-14 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. % 2010-01-18 updated to conform better to Bristow-Johnson's definition omega = 2 * pi * fc / fs; alpha = sin(omega)/(2*q); B = [ (1 - cos(omega)) / 2 1 - cos(omega) (1 - cos(omega)) / 2 ]; A = [ 1 + alpha -2 * cos(omega) 1 - alpha ];
function [B,A] = biquad_hpf(fc, q, fs) %BIQUAD_HPF biquad high-pass filter % [B,A] = biquad_hpf(fc, q, fs) % % Input: % fc - cutoff frequency (Hz) % q - Q-factor % fs - sampling frequency (Hz) % % 2009-11-30 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. % 2010-01-18 updated to conform better to Bristow-Johnson's definition omega = 2 * pi * fc / fs; alpha = sin(omega)/(2*q); B = [ (1 + cos(omega)) / 2 -(1 + cos(omega)) (1 + cos(omega)) / 2 ]; A = [ 1 + alpha -2 * cos(omega) 1 - alpha ];
function [B,A] = biquad_bpf(fc, bw, fs) % BIQUAD_BPF biquad band-pass filter % [B,A] = biquad_bpf(fc, bw, fs) % % Input: % fc - center frequency (Hz) % bw - bandwidth at -3dB point in octaves % fs - sampling frequency (Hz) % % 2005-12-16 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. omega = 2 * pi * fc / fs; alpha = sin(omega) * sinh(log(2) / 2 * bw * omega / sin(omega)); B = [ alpha 0 -alpha ]; A = [ 1 + alpha -2 * cos(omega) 1 - alpha ];
function [B,A] = biquad_lowshelf(fc, gain, s, fs)
%BIQUAD_LOWSHELF Biquad low-shelf filter
% [B,A] = biquad_lowshelf(fc, gain, s, fs)
%
% Input:
% fc - shelf midpoint frequency (Hz)
% gain - gain at the midpoint (dB)
% s - shelf slope (S=1 for steepest slope)
% fs - sampling frequency (Hz)
%
% 2009-09-13 by MARUI Atsushi
% based on "Cookbook formulae for audio EQ biquad filter
% coefficients" by Robert Bristow-Johnson.
omega = 2 * pi * fc / fs;
a = 10 ^ (gain/40);
alpha = sin(omega) / 2 * sqrt((a + 1/a) * (1/s - 1) + 2);
B = [
a * ((a+1) - (a-1)*cos(omega) + 2*sqrt(a)*alpha)
2*a * ((a-1) - (a+1)*cos(omega) )
a * ((a+1) - (a-1)*cos(omega) - 2*sqrt(a)*alpha)
]';
A = [
((a+1) + (a-1)*cos(omega) + 2*sqrt(a)*alpha)
-2 * ((a-1) + (a+1)*cos(omega) )
((a+1) + (a-1)*cos(omega) - 2*sqrt(a)*alpha)
]';
function [B,A] = biquad_highshelf(fc, gain, s, fs)
%BIQUAD_HIGHSHELF Biquad high-shelf filter
% [B,A] = biquad_highshelf(fc, gain, s, fs)
%
% Input:
% fc - shelf midpoint frequency (Hz)
% gain - gain at the midpoint (dB)
% s - shelf slope (S=1 for steepest slope)
% fs - sampling frequency (Hz)
%
% 2010-01-18 by MARUI Atsushi
% based on "Cookbook formulae for audio EQ biquad filter
% coefficients" by Robert Bristow-Johnson.
omega = 2 * pi * fc / fs;
a = 10 ^ (gain/40);
alpha = sin(omega) / 2 * sqrt((a + 1/a) * (1/s - 1) + 2);
B = [
a * ((a+1) + (a-1)*cos(omega) + 2*sqrt(a)*alpha)
-2*a * ((a-1) + (a+1)*cos(omega) )
a * ((a+1) + (a-1)*cos(omega) - 2*sqrt(a)*alpha)
]';
A = [
((a+1) - (a-1)*cos(omega) + 2*sqrt(a)*alpha)
2 * ((a-1) - (a+1)*cos(omega) )
((a+1) - (a-1)*cos(omega) - 2*sqrt(a)*alpha)
]';
function [B,A] = biquad_peak(fc, gain, q, fs) %BIQUAD_PEAK Biquad peak/dip filter % [B,A] = biquad_peak(fc, gain, q, fs) % % Input: % fc - cutoff frequency (Hz) % gain - gain at the peak (dB) % q - Q factor % fs - sampling frequency (Hz) % % 2005-10-13 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. omega = 2 * pi * fc / fs; alpha = sin(omega) / q; a = 10 ^ (gain/40); B = [ 1 + alpha * a -2 * cos(omega) 1 - alpha * a ]'; A = [ 1 + alpha / a -2 * cos(omega) 1 - alpha / a ]';
function [B,A] = biquad_notch(fc, bw, fs) %BIQUAD_NOTCH Biquad notch filter % [B,A] = biquad_notch(fc, bw, fs) % % Input: % fc - cutoff frequency (Hz) % bw - bandwidth in octaves (between -3dB frequencies) % fs - sampling frequency (Hz) % % 2010-01-18 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. omega = 2 * pi * fc / fs; alpha = sin(omega) * sinh(log(2) / 2 * bw * omega / sin(omega)); B = [ 1 -2 * cos(omega) 1 ]'; A = [ 1 + alpha -2 * cos(omega) 1 - alpha ]';
function [B,A] = biquad_apf(fc, q, fs) % BIQUAD_APF biquad all-pass filter % [B,A] = biquad_apf(fc, q, fs) % % Input: % fc - center frequency (Hz) % q - Q-factor % fs - sampling frequency (Hz) % % 2010-01-18 by MARUI Atsushi % based on "Cookbook formulae for audio EQ biquad filter % coefficients" by Robert Bristow-Johnson. omega = 2 * pi * fc / fs; alpha = sin(omega)/(2*q); B = [ 1 - alpha -2 * cos(omega) 1 + alpha ]; A = [ 1 + alpha -2 * cos(omega) 1 - alpha ];
帯域ごとに分割し、また元の信号に戻せるようなフィルターバンクが欲しいときがあります。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