% PROTO2DBP
% prototype to bandpass disrete transformation
% [Bz Az] = proto2dbp(Zs, Ps, Ks, f0, fBWpm, fsamp)
% Z,P,K - analog/continous normalized lowpass prototype
% f0 - center freq [Hz]
% fBWpm - bandwidth [Hz] e.g. fBWpm=10 => +/-5Hz
% fsamp - sampling freq [Hz]
%
% auth: Petr Sladek CTU FEE, id sladep1

%function [Bs As Bz Az] = proto2dbp(Zs, Ps, Ks, f0, fBWpm, fsamp);

%[Zs Ps Ks] = ellipap(2,0.22,30);
wBWpm = fBWpm*2*pi;
w0= f0*2*pi;
% zeros
Bs = real(poly(Zs))*Ks;
% poles
As = real(poly(Ps));
%bandpass from prototype plus freq shift:
[Bs As]=lp2bp(Bs,As,w0,wBWpm);


%residue:

[R P K] = residue(Bs,As); %decomposition

Tsamp = 1/fsamp;

% Ra contains numerator of Ri*1/(1-exp(PkT)*z^-1)
% Pa ----------------------------------^
% Ka

%manual impulse invariance:
Rz = R.*Tsamp; %("1/T" scale)
Pz = exp(P.*Tsamp);
Kz = K;

[Bz Az]=residuez(Rz,Pz,Kz); %back to poly
Bz = real(Bz); %trash small im parts, potential problem!
Az = real(Az); %trash small im parts, potential problem!