%% 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{BT AT } =} tftrans (@var{BP}, @var{AP}, @var{W}, @var{STOP})
%%
%% Transform band edges of a generic lowpass filter (cutoff at W=1)
%% represented in splane transfer function form. W is the edge of the
%% target filter (or edges if band pass or band stop). Stop is true for
%% high pass and band stop filters or false for low pass and band pass
%% filters. Filter edges are specified in radians/second.
%%
%% @seealso{lp2lp, lp2hp, lp2bp, lp2bs}
%% @end deftypefn
%% Author: Tony Richardson
%% Created: 2021-04-10
function [BT AT] = tftrans (BP, AP, W, stop)
if (nargin != 4)
print_usage;
endif
BP = BP(:)'; AP = AP(:)';
N = length(AP); M = length(BP);
if length(W)==2
if stop
%% Band Stop
%% S*BW
%% S -> --------
%% S^2+FhFl
NBP = 2*(N - 1); MBP = 2*(M - 1); % Polynomial orders
v1 = [ abs(W(2)-W(1)) 0 ]; v2 = [1 0 W(2)*W(1)];
AT = 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
AT = AT + [zeros(1,NBP+1-length(vT)) AP(n)*vT];
endfor
BT = 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
BT = BT + [zeros(1,MBP+1-length(vT)) BP(n)*vT];
endfor
if (M < N)
for k = 1:(N-M)
BT = conv(BT, v2);
end
elseif (N < M)
for k = 1:(M-N)
AT = conv(AT, v2);
end
endif
else
%% Band Pass
%% S^2+FhFl
%% S -> --------
%% S*BW
NBP = 2*(N - 1); MBP = 2*(M - 1); % Polynomial orders
v1 = [1 0 W(2)*W(1)]; v2 = [abs(W(2)-W(1)) 0 ];
AT = 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
AT = AT + [zeros(1,NBP+1-length(vT)) AP(n)*vT];
endfor
BT = 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
BT = BT + [zeros(1,MBP+1-length(vT)) BP(n)*vT];
endfor
if (M < N)
for k = 1:(N-M)
BT = conv(BT, v2);
end
elseif (N < M)
for k = 1:(M-N)
AT = conv(AT, v2);
end
endif
endif
else
Fc = W;
if stop
%% High Pass
%% S -> Fc/S
% The adjustment is done in such a manner so that if AP(N) = 1
% then AT(1) = 1. This should help to avoid underflow/overflow problems.
n = 0:(M-1);
BT = [fliplr(BP).*Fc.^(n) zeros(1,N-M)];
% Adjust denominator coefficients
n = 0:(N-1);
AT = [fliplr(AP).*Fc.^(n) zeros(1,M-N)];
else
%% Low Pass
%% S -> S/Fc
% The adjustment is done in such a manner so that if AP(1) = 1
% then AT(1) = 1. This should help to avoid underflow/overflow problems.
% Adjust numerator coefficients
n = 0:(M-1);
BT = Fc^(N-M)*(BP .* Fc.^(n));
% Adjust denominator coefficients
n = 0:(N-1);
AT = AP .* Fc.^(n);
endif
endif
% Normalize the coefficients
BT = BT/AT(1); AT = AT/AT(1);
endfunction