バンドパスフィルター（双二次） / Biquad Bandpass Filter

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
];