%% 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