butterworth bandpass filter coefficient calculatinon (algorithm)

[EDIT: 20110720 22:53 CDT - convert character set, clarify - WDR]

I am trying to code my own butterworth bandpass filter to be implemented in C, from discrete 2th order lowpass butterworth filter by lowpass to bandpass spectral transformation, but the filter coefficient is not calculated correctly compared with the Matlab butter function. Following is the complete code of my program. Just copy to the Matlab editor and you can run it.

The algorithm is based on the following steps:

1) 2th order continuous Laplace transform is: 1/(s^2+ sqrt(2)*s + 1)

2) then transform the above filter to digital from by bilinear transform; s = 1/k*(z-1)/(z+1), where k = tan( pi * fc / fs); %fc is the cut-off frequency of the lowpass filter

3) lowpass to bandpass transform z--> (z^2 - r *z)/(r*z-1), where r = cos(w0) = cos( 2 * pi * ( fc1 + fc2 ) / 2 / fs ) / cos( 2 * pi * fc / 2 / fs); fc = fc2-fc1

 based on the following paper
Constantinides, A.G., "Design of bandpass digital filters,' IEEE Proceedings, vol. 1, pp. 1129-1231, June 1969.

From 1-3), the coefficients of this IIR filter can be determined.

Are there any one can help check where I am I wrong? Thanks.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function test_bp1

fc1 = 1e3; fc2 = 10e3; fs = 48e3;

[b, a] =  butter2_bp_direct(fc1, fc2,  fs)
figure, freqz(b, a)

[b, a] =  butter(2, [fc1, fc2]*2/fs)
figure,freqz(b, a)

function [b1, a1] =  butter2_bp_direct(fc1, fc2,  fs)

q1 = sqrt(2); 

fc = abs( fc2 - fc1 );% equivalent cut-off frequency of lowpass filter is (fc2-fc1)
k = tan( pi * fc / fs);

r = cos( 2 * pi * ( fc1 + fc2 ) / 2 / fs ) / cos( 2 * pi * fc / 2 / fs);

b1 = k^2 * [1, 0, -2, 0, 1];

a10 =   k^4 + k*q1  + 1;
a11 =   -2*k *q1*r    - 4*r;
a12 =   -2 *k*k  +4*r*r  + 2;
a13 =   2*k *q1*r       -    4*r;
a14 =  k^2-k*q1+1;
a1 = [a10, a11, a12, a13, a14]/a10;
b1 = b1/a10;

[Walter Roberson]

Please do not use unicode characters such as "ｉ" (unicode FF49) when posting text. I had to type your text in by hand to correct it here, as my editors would not convert it for me.

[Jan]

BTW. "close all; clear all;" is a waste of time when called at the start of a function. I strongly recommend to avoid wild clearing of all functions loaded into memory, because it can reduce the runtime by a factor of 10.000.

[Victor Lee]

Thanks, Jan!