%% Copyright (C) 2021 Tony Richardson %% %% This program is free software: you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published by %% the Free Software Foundation, either version 3 of the License, or %% (at your option) any later version. %% %% This program is distributed in the hope that it will be useful, %% but WITHOUT ANY WARRANTY; without even the implied warranty of %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the %% GNU General Public License for more details. %% %% You should have received a copy of the GNU General Public License %% along with this program. If not, see . %% -*- texinfo -*- %% @deftypefn {} {@var{[B, A]} =} lp2bp (@var{Bp}, @var{Ap}, @var{w0}, @var{W}) %% Transform a lowpass s-plane filter into an s-plane bandpass filter. %% Bp and Ap are numerator and denominator coefficients of the lowpass %% prototype filter. w0 and W are the center frequency and bandwidth %% (rad/s) of the desired bandpass filter. The numerator and denominator %% coefficients of the desired bandpass filter are returned. %% %% @seealso{lp2lp, lp2hp, lp2bs} %% @end deftypefn %% Author: Tony Richardson %% Created: 2021-04-09 function [B, A] = lp2bp(Bp, Ap, w0, W) N = length(Ap); M = length(Bp); NBP = 2*(N - 1); MBP = 2*(M - 1); v1 = [1 0 w0^2]; v2 = [ W 0 ]; A = zeros(1, NBP + 1); for n = 1:N vT = 1; for k = 1:(N-n) vT = conv(vT, v1); endfor for k = 1:(n-1) vT = conv(vT, v2); endfor A = A + [zeros(1,NBP+1-length(vT)) Ap(n)*vT]; endfor B = zeros(1, MBP + 1); for n = 1:M vT = 1; for k = 1:(M-n) vT = conv(vT, v1); endfor for k = 1:(n-1) vT = conv(vT, v2); endfor B = B + [zeros(1,MBP+1-length(vT)) Bp(n)*vT]; endfor while (M < N) B = conv(B, v2); M = M + 1; end NORM = A(1); A = A/NORM; B = B/NORM; end %!demo %! n = 14; %! [z p k] = buttap(n); %! [b a] = zp2tf(z, p, real(k)); %! %! w = logspace(-1, 1, 501); %! H = freqs(b, a, w); %! figure(1) %! clf %! subplot(2, 1, 1) %! loglog(w, abs(H)) %! ylim([1e-15 1]); %! yticks(logspace(-15, 0, 4)) %! grid on %! subplot(2, 1, 2) %! semilogx(w, rad2deg(angle(H))) %! grid on %! %! fl = 30; fh = 100; %! Wo = 2*pi*sqrt(fl*fh); %! Bw = 2*pi*(fh - fl); %! [bt, at] = lp2bp(b, a, Wo, Bw); %! %! w = logspace(0, 4, 501); %! H = freqs(bt, at, w); %! %! figure(2) %! clf %! subplot(2, 1, 1) %! loglog(w, abs(H)) %! ylim([1e-20 1]); %! yticks(logspace(-20, 0, 3)) %! grid on %! subplot(2, 1, 2) %! semilogx(w, rad2deg(angle(H))) %! grid on