% lhbutt.m - lowpass/highpass Butterworth digital filter design%% [A, B, P] = lhbutt(s, fs, fpass, fstop, Apass, Astop)%% s = 1, -1 = lowpass, highpass% design parameters:% P = [Wpass, Wstop, epass, estop, Nex, N, Astop, W0, f0];% A, B are Kx3 matrices of cascade second-order sectionsfunction [A, B, P] = lhbutt(s, fs, fpass, fstop, Apass, Astop)Wpass = tan(pi * fpass / fs);  Wpass = Wpass^s;         % cot() for HPWstop = tan(pi * fstop / fs);  Wstop = Wstop^s;epass = sqrt(10^(Apass/10) - 1); estop = sqrt(10^(Astop/10) - 1);Nex = log(estop/epass) / log(Wstop/Wpass);N = ceil(Nex);  r = rem(N,2);  K = (N - r) / 2;         % K = no. sectionsW0 = Wpass * (epass^(-1/N));Astop = 10 * log10(1 + (Wstop/W0)^(2*N));               % actual Astopf0 = (fs/pi) * atan(W0^s);                              % 3-dB freq. in HzP = [Wpass, Wstop, epass, estop, Nex, N, Astop, W0, f0];if r==1,                                                % N = odd    G = W0 / (1 + W0);                                  % 1st order section    B(1,:) = G * [1, s, 0];    A(1,:) = [1, s*(2*G-1), 0];else                                                    % N = even    B(1,:) = [1, 0, 0];    A(1,:) = [1, 0, 0];endfor i=1:K,    th = pi * (N - 1 + 2 * i) / (2 * N);    D = 1 - 2 * W0 * cos(th) + W0^2;    G = W0^2 / D;    a1 = 2 * (W0^2 - 1) / D;    a2 = (1 + 2 * W0 * cos(th) + W0^2) / D;    B(i+1,:) = G * [1, 2*s, 1];    A(i+1,:) = [1, s*a1, a2];end