function [b,a] = afd_elip(Wp,Ws,Rp,As); % Analog Lowpass Filter Design: Elliptic % -------------------------------------- % [b,a] = afd_elip(Wp,Ws,Rp,As);
% b = Numerator coefficients of Ha(s) % a = Denominator coefficients of Ha(s)
% Wp = Passband edge frequency in rad/sec; Wp > 0
% Ws = Stopband edge frequency in rad/sec; Ws > Wp > 0 % Rp = Passband ripple in +dB; (Rp > 0)
% As = Stopband attenuation in +dB; (As > 0) %
if Wp <= 0
error('Passband edge must be larger than 0') end
if Ws <= Wp
error('Stopband edge must be larger than Passband edge') end
if (Rp <= 0) | (As < 0)
error('PB ripple and/or SB attenuation ust be larger than 0') end
ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); OmegaC = Wp; k = Wp/Ws;
k1 = ep/sqrt(A*A-1);
capk = ellipke([k.^2 1-k.^2]); % Version 4.0 code
capk1 = ellipke([(k1 .^2) 1-(k1 .^2)]); % Version 4.0 code N = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));
fprintf('\\n*** Elliptic Filter Order = %2.0f \\n',N) [b,a]=u_elipap(N,Rp,As,OmegaC); function w_black = Blackman(M); % M-point Blackman window % ----------------------- % w_black = Blackman(M); %
M1 = M-1; m = [0:1:M1];
w_black = abs(0.42 - 0.5*cos(2*pi*m'/(M1)) + 0.08*cos(4*pi*m'/(M1)));
function [b,a] = cas2dir(b0,B,A); % CASCADE-to-DIRECT form conversion % ---------------------------------
26
% [b,a] = cas2dir(b0,B,A)
% b = numerator polynomial coefficients of DIRECT form % a = denominator polynomial coefficients of DIRECT form % b0 = gain coefficient
% B = K by 3 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's %
[K,L] = size(B); b = [1]; a = [1]; for i=1:1:K
b=conv(b,B(i,:)); a=conv(a,A(i,:)); end
b = b*b0;
function y = casfiltr(b0,B,A,x);
% CASCADE form realization of IIR and FIR filters % ----------------------------------------------- % y = casfiltr(b0,B,A,x); % y = output sequence
% b0 = gain coefficient of CASCADE form
% B = K by 3 matrix of real coefficients containing bk's % A = K by 3 matrix of real coefficients containing ak's % x = input sequence %
[K,L] = size(B); N = length(x); w = zeros(K+1,N); w(1,:) = x; for i = 1:1:K
w(i+1,:) = filter(B(i,:),A(i,:),w(i,:)); end
y = b0*w(K+1,:);
function [b,a] = cheb1hpf(wp,ws,Rp,As)
% IIR Highpass filter design using Chebyshev-1 prototype % [b,a] = cheb1hpf(wp,ws,Rp,As)
% b = Numerator polynomial of the highpass filter % a = Denominator polynomial of the highpass filter % wp = Passband frequency in radians % ws = Stopband frequency in radians % Rp = Passband ripple in dB
27
% As = Stopband attenuation in dB %
% Determine the digital lowpass cutoff frequecies: wplp = 0.2*pi;
alpha = -(cos((wplp+wp)/2))/(cos((wplp-wp)/2));
wslp = angle(-(exp(-j*ws)+alpha)/(1+alpha*exp(-j*ws))); %
% Compute Analog lowpass Prototype Specifications: T = 1; Fs = 1/T;
OmegaP = (2/T)*tan(wplp/2); OmegaS = (2/T)*tan(wslp/2);
% Design Analog Chebyshev Prototype Lowpass Filter: [cs,ds] = afd_chb1(OmegaP,OmegaS,Rp,As);
% Perform Bilinear transformation to obtain digital lowpass [blp,alp] = bilinear(cs,ds,Fs);
% Transform digital lowpass into highpass filter Nz = -[alpha,1]; Dz = [1,alpha]; [b,a] = zmapping(blp,alp,Nz,Dz); function [xec, xoc] = circevod(x)
% signal decomposition into circular-even and circular-odd parts% --------------------------------------------------------------% [xec, xoc] = circecod(x) %
if any(imag(x) ~= 0)
error('x is not a real sequence') end
N = length(x); n = 0:(N-1); xec = 0.5*(x + x(mod(-n,N)+1)); xoc = 0.5*(x - x(mod(-n,N)+1));
function y = circonvt(x1,x2,N)
% N-point circular convolution between x1 and x2: (time-domain) % ------------------------------------------------------------- % [y] = circonvt(x1,x2,N)
% y = output sequence containing the circular convolution % x1 = input sequence of length N1 <= N % x2 = input sequence of length N2 <= N % N = size of circular buffer
% Method: y(n) = sum (x1(m)*x2((n-m) mod N))
28
% Check for length of x1 if length(x1) > N
error('N must be >= the length of x1') end
% Check for length of x2 if length(x2) > N
error('N must be >= the length of x2') end
x1=[x1 zeros(1,N-length(x1))]; x2=[x2 zeros(1,N-length(x2))];
m = [0:1:N-1];
x2 = x2(mod(-m,N)+1); H = zeros(N,N);
for n = 1:1:N
H(n,:) = cirshftt(x2,n-1,N); end
y = x1*H';
function y = cirshftt(x,m,N)
% Circular shift of m samples wrt size N in sequence x: (time domain) % ------------------------------------------------------------------- if length(x) > N
error('N must be >= the length of x') end
x = [x zeros(1,N-length(x))]; n = [0:1:N-1]; n = mod(n-m,N); y = x(n+1);
function [y,ny] = conv_m(x,nx,h,nh)
% Modified convolution routine for signal processing % -------------------------------------------------- % [y,ny] = conv_m(x,nx,h,nh) % y = convolution result % ny = support of y
% x = first signal on support nx % nx = support of x
% h = second signal on support nh
29
% nh = support of h %
nyb = nx(1)+nh(1); nye = nx(length(x)) + nh(length(h)); ny = [nyb:nye]; y = conv(x,h);
function [y,H]=conv_tp(h,x)
% Linear Convolution using Toeplitz Matrix % ---------------------------------------- % [y,H] = conv_tp(h,x)
% y = output sequence in column vector form
% H = Toeplitz matrix corresponding to sequence h so that y = Hx % h = Impulse response sequence in column vector form % x = input sequence in column vector form %
Nx = length(x); Nh = length(h); hc=[h; zeros(Nx-1, 1)]; hr=[h(1),zeros(1,Nx-1)]; H=toeplitz(hc,hr); y=H*x;
function I = cplxcomp(p1,p2) % I = cplxcomp(p1,p2)
% Compares two complex pairs which contain the same scalar elements % but (possibly) at differrent indices. This routine should be % used after CPLXPAIR routine for rearranging pole vector and its % corresponding residue vector. % p2 = cplxpair(p1) % I=[];
for j=1:1:length(p2) for i=1:1:length(p1)
if (abs(p1(i)-p2(j)) < 0.0001) I=[I,i]; end end end I=I';
function [d1,d2] = db2delta(Rp,As)
% Conversion from Relative dB specs to Absolute delta specs. % ---------------------------------------------------------- % [d1,d2] = db2delta(Rp,As)
30