No BSD License  

Impulse response invariant discretization of fractional order low-pass filters



07 Sep 2008 (Updated )

Discretize [1/(\tau s +1)]^r with "r" a real number

% Impulse response invariant discretization of fractional order 
% low-pass filters
% irid_folpf function is prepared to compute a discrete-time finite 
% dimensional (z) transfer function to approximate a continuous-time 
% fractional order low-pass filter (LPF) [1/(\tau s +1)]^r, where "s" is 
% the Laplace transform variable, and "r" is a real number in the range of 
% (0,1), \tau is the time constant of LPF [1/(\tau s +1)]. 
% The proposed approximation keeps the impulse response "invariant"
% IN: 
%       tau: the time constant of (the first order) LPF
%       r: the fractional order \in (0,1)
%       Ts: the sampling period
%       norder: the finite order of the approximate z-transfer function 
%       (the orders of denominator and numerator z-polynomial are the same)
% OUT: 
%       sr: returns the LTI object that approximates the [1/(\tau s +1)]^r
%           in the sense of invariant impulse response. 
% dfod=irid_folpf(.01,0.5,.001,5);figure;pzmap(dfod)
% Reference: YangQuan Chen. "Impulse-invariant discretization of fractional
% order low-pass filters".
% Sept. 2008. CSOIS AFC (Applied Fractional Calculus) Seminar.
% --------------------------------------------------------------------
% YangQuan Chen, Ph.D, Associate Professor and Graduate Coordinator
% Department of Electrical and Computer Engineering,
% Director, Center for Self-Organizing and Intelligent Systems (CSOIS)
% Utah State University, 4120 Old Main Hill, Logan, UT 84322-4120, USA
% E: or, T/F: 1(435)797-0148/3054; 
% W: or 
% --------------------------------------------------------------------
% 9/7/2009 
% Only supports when r in (0,1). That is fractional order low pass filter.
% HOWEVER, if r is in (-1,0), we call this is a "fractional order
% (proportional and derivative controller)" - we call it FO(PD).
% Note: it may be needed to make FO-LPF discretization minimum phase first.
% See also irid_fod.m 
%          at
% See also srid_fod.m 
%          (See how the nonminimum phase zeros are handled)
% See also gml_fun.m 
%          at

function [sr]=irid_folpf(tau,r,Ts,norder)

if nargin<4; norder=5; end
if tau < 0 , sprintf('%s','Time constant has to be positive'),     return, end
if Ts < 0 , sprintf('%s','Sampling period has to be positive'),     return, end
if r>=1 | r<= 0, sprintf('%s','The fractional order should be in (0,1)'), return, end
if norder<2, sprintf('%s','The order of the approximate transfer function has to be greater than 1'), return, end
wmax0=2*pi/Ts/2; % rad./sec. Nyquist frequency
L=abs(tau)*4/Ts; % decides the number of points of the impulse response function h(n)
n=1:L-1; n=n*Ts ;  
h1=gml_fun(1,r,r,-1/tau *n); %
h0= ha0*(1/(tau+ (7.0*Ts/8)))^r;   
h=[h0,h2*Ts]; %% [ha0, (Ts^r)*(n.^(r-1))/gamma(r)]; 
% --  visualize the approximation
% approximated h()
hold on;
hold on; 
plot(Taxis,ioir*Ts,'g:'); % compare the r=1 case - exponential function)
xlabel('time');ylabel('impulse response');
legend(['approximated  for 1/(',num2str(tau), '* s +1 )^{',num2str(abs(r)),'}'],'true','first order low pass filter')

wmax=floor(1+ log10(wmax0) ); wmin=wmax-5;
srfr=(1./(tau*j*w +1)).^(r);

%sc=tf(1,[tau,1]); sd=c2d(sc,Ts,'tustin');
semilogx(w,20*log10(abs(srfr)),'r');grid on
hold on;
semilogx(w,20*log10(abs(srfr)),'b:');grid on
semilogx(w,20*log10(abs(reshape(srfrhat, 1000, 1))),'k');grid on
xlabel('frequency in Hz');ylabel('dB');
legend('true mag. Bode - continuous time','true mag. Bode - discrete time','approximated mag. Bode')

semilogx(w,(180/pi) * (angle(srfr)),'r');grid on;hold on
semilogx(w,(180/pi) * (angle(srfr1)),'b:');grid on;hold on
semilogx(w,(180/pi) * (angle(reshape(srfrhat, 1000, 1))),'k');grid on
xlabel('frequency in Hz');ylabel('phase (degrees)');
legend('true phase Bode - continuous time','true phase Bode - discrete time','approximated phase Bode')

% % Want to discretize FOPD (1+\tau s)^r ??
% % use sr1=1/sr, with danger.
% % Try this - get stable, minimum phase approximation.
% [zz,pp,kk]=zpkdata(sr,'v');
% for i=1:norder;
%     if abs(zz(i)) > 1
%         kk=kk*(-zz(i)); zz(i)=1/zz(i); 
%         sprintf('%s','nonminimum phase approximation - forced minimum phase!!'), 
%     end
%     if abs(pp(i)) > 1
%         kk=kk/(-pp(i)); pp(i)=1/pp(i); 
%         sprintf('%s','unstable approximation - forced stable!!'), 
%     end
% end
% sr1=zpk(zz,pp,kk,Ts);

Contact us