Global Index (short | long) | Local contents | Local Index (short | long)
[zf,thf]=idfilt(z,n,Wn,hs)
IDFILT Filters data.
ZF = IDFILT(Z,N,Wn)
Z : The output-input data matrix Z = [Y U]. Multi-variable
data are allowed
ZF : The filtered data. The columns of ZF correspond to those of Z.
N : The order of the filter to be used.
N = a positive integer: gives an acausal, zero-phase filter of
order 2N, determined from a Butterworth filter.
N = a negative integer: gives a causal Butterworth filter of
order abs(N).
Wn : The cut-off frequency(ies) in fractions of the Nyquist frequency.
If Wn is a scalar a Low-Pass filter will be used.
If Wn = [Wl Wh] a Band-Pass filter with pass-band Wn will be used
With ZF = IDFILT(Z,N,Wn,'high') or ZF = IDFILT(Z,N,[Wl Wh],'stop')
High-Pass and Band-Stop filters will be used instead.
With [ZF,THFILT] = IDFILT(Z,N,Wn), also the filter will be returned
in the THETA-format as THFILT. (See also THETA).
See also DTREND and IDRESAMP.
| This function calls | This function is called by |
|---|---|
function [zf,thf]=idfilt(z,n,Wn,hs)
% L. Ljung 10-1-89, 3-3-95.
% The code is adopted from several routines in the Signal Processing
% Toolbox.
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 3.3 $ $Date: 1997/12/02 03:42:26 $
if nargin<3
disp('Usage: zf = IDFILT(DATA,ORDER,BAND);')
disp(' zf = IDFILT(DATA,ORDER,BAND,''HIGH'');')
disp(' with ORDER = an integer (positive for acausal filtering).')
return
end
if nargin<4,hs=[];end
if length(Wn)>1,
if Wn(1)>=Wn(2)
error('The band is ill-defined. Should be Wn = [a b] with a < b.')
end
if abs(Wn(1))<eps
Wn=Wn(2);
elseif abs(Wn(2)-1)<eps,
Wn=Wn(1);
if isempty(hs),hs='high';else hs=[];end
end
end
[ndat,nyu]=size(z);
if ndat<nyu
error('Data should be organized columnwise.')
end
if ndat<3*n
error('The data record length must be at least three times the filter order')
end
if n<0,n=abs(n);causal=1;else causal=0;end
btype = 1;
if ~isempty(hs),btype=3;end
if length(Wn) == 2
btype = btype + 1;
end
% step 1: get analog, pre-warped frequencies
fs = 2;
u = 2*fs*tan(pi*Wn/fs);
% step 2: convert to low-pass prototype estimate
if btype == 1 % lowpass
Wn = u;
elseif btype == 2 % bandpass
Bw = u(2) - u(1);
Wn = sqrt(u(1)*u(2)); % center frequency
elseif btype == 3 % highpass
Wn = u;
elseif btype == 4 % bandstop
Bw = u(2) - u(1);
Wn = sqrt(u(1)*u(2)); % center frequency
end
% step 3: Get N-th order Butterworth analog lowpass prototype
p = exp(sqrt(-1)*(pi*(1:2:2*n-1)/(2*n) + pi/2)).';
k = real(prod(-p));
% Transform to state-space
% Strip infinities and throw away.
p = p(finite(p));
% Group into complex pairs
np = length(p);
nz= 0;%= length(z);
p = cplxpair(p,1e6*np*norm(p)*eps + eps);
% Initialize state-space matrices for running series
a=[]; b=[]; c=zeros(1,0); d=1;
% If odd number of poles only, convert the pole at the
% end into state-space.
% H(s) = 1/(s-p1) = 1/(s + den(2))
if rem(np,2)
a = p(np);
b = 1;
c = 1;
d = 0;
np = np - 1;
end
% Now we have an even number of poles and zeros, although not
% necessarily the same number - there may be more poles%.
% H(s) = (s^2+num(2)s+num(3))/(s^2+den(2)s+den(3))
% Loop thru rest of pairs, connecting in series to build the model.
i = 1;
% Take care of any left over unmatched pole pairs.
% H(s) = 1/(s^2+den(2)s+den(3))
while i < np
den = real(poly(p(i:i+1)));
wns = sqrt(prod(abs(p(i:i+1))));
if wns == 0, wns = 1; end
t = diag([1 1/wns]); % Balancing transformation
a1 = t\[-den(2) -den(3); 1 0]*t;
b1 = t\[1; 0];
c1 = [0 1]*t;
d1 = 0;
% [a,b,c,d] = series(a,b,c,d,a1,b1,c1,d1);
% Next lines perform series connection
[ma1,na1] = size(a);
[ma2,na2] = size(a1);
a = [a zeros(ma1,na2); b1*c a1];
b = [b; b1*d];
c = [d1*c c1];
d = d1*d;
i = i + 2;
end
% Apply gain k:
c = c*k;
d = d*k;
% step 4: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
wo = Wn;
if btype == 1 % Lowpass
at = wo*a;
bt = wo*b;
ct = c;
dt = d;
elseif btype == 2 % Bandpass
[ma,nb] = size(b);
[mc,ma] = size(c);
% Transform lowpass to bandpass
q = wo/Bw;
at = wo*[a/q eye(ma); -eye(ma) zeros(ma)];
bt = wo*[b/q; zeros(ma,nb)];
ct = [c zeros(mc,ma)];
dt = d;
elseif btype == 3 % Highpass
at = wo*inv(a);
bt = -wo*(a\b);
ct = c/a;
dt = d - c/a*b;
elseif btype == 4 % Bandstop
[ma,nb] = size(b);
[mc,ma] = size(c);
% Transform lowpass to bandstop
q = wo/Bw;
at = [wo/q*inv(a) wo*eye(ma); -wo*eye(ma) zeros(ma)];
bt = -[wo/q*(a\b); zeros(ma,nb)];
ct = [c/a zeros(mc,ma)];
dt = d - c/a*b;
end
a=at;b=bt;c=ct;d=dt;
% step 5: Use Bilinear transformation to find discrete equivalent:
t = 1/fs;
r = sqrt(t);
t1 = eye(size(a)) + a*t/2;
t2 = eye(size(a)) - a*t/2;
ad = t2\t1;
bd = t/r*(t2\b);
cd = r*c/t2;
dd = c/t2*b*t/2 + d;
a = ad;b = bd; c = cd; d = dd;
den = poly(a);
num = poly(a-b*c)+(d-1)*den;
if nargout>1, thf=mktheta(den,num);end
if causal
for k=1:nyu
x=ltitr(a,b,z(:,k));
zf(:,k)=x*c'+d*z(:,k);
end
else
for k = 1:nyu
x=z(:,k);
len = size(x,1); % length of input
b = num(:).';
a = den(:).';
nb = length(b);
na = length(a);
nfilt = max(nb,na);
nfact = 3*(nfilt-1); % length of edge transients
if nb < nfilt, b(nfilt)=0; end % zero-pad if necessary
if na < nfilt, a(nfilt)=0; end
zi = ( eye(nfilt-1) - [-a(2:nfilt).' ...
[eye(nfilt-2); zeros(1,nfilt-2)]] ) \ ...
( b(2:nfilt).' - a(2:nfilt).'*b(1) );
y = [2*x(1)-x((nfact+1):-1:2);x;...
2*x(len)-x((len-1):-1:len-nfact)];
% filter, reverse data, filter again, and reverse data again
y = filter(b,a,y,[zi*y(1)]);
y = y(length(y):-1:1);
y = filter(b,a,y,[zi*y(1)]);
y = y(length(y):-1:1);
% remove extrapolated pieces of y
y([1:nfact len+nfact+(1:nfact)]) = [];
zf(:,k)=y;
end
end