view m-toolbox/classes/+utils/@math/pfresp.m @ 24:056f8e1e995e
database-connection-manager
Properly record history in fromRepository constructors
author
Daniele Nicolodi <nicolodi@science.unitn.it>
date
Mon, 05 Dec 2011 16:20:06 +0100 (2011-12-05)
parents
f0afece42f48
children
line source
+ − % PFRESP returns frequency response of a partial fraction TF.
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION:
+ − %
+ − % Returns frequency response of a partial fraction expanded function
+ − % (continuous or discrete).
+ − %
+ − % Continuous case
+ − % The expected model is:
+ − %
+ − % r1 rN
+ − % f(s) = ------ + ... + ------ + d1 + d2*s + ... + dK*s^{K-1}
+ − % s - p1 s - p1
+ − %
+ − % Discrete case
+ − % The expected model is:
+ − %
+ − % z*r1 z*rN
+ − % f(z) = ------ + ... + ------ + d1 + d2*z^{-1} + ... + dK*z^{-(K-1)}
+ − % z - p1 z - p1
+ − %
+ − % NOTE: The function cannot handle poles multiplicity higher than 1 in
+ − % z domain. Multiple poles in s-domain are accepted.
+ − %
+ − % CALL: pfr = pfresp(pfparams)
+ − %
+ − % INPUT:
+ − %
+ − % pfparams is a struct containing input parameters
+ − % pfparams.type = 'cont' Assumes a continuous model
+ − % pfparams.type = 'disc' Assumes a discrete model
+ − % pfparams.freq set the frequencies vector in Hz
+ − % pfparams.res - set the vector of residues
+ − % pfparams.pol - set the vector of poles
+ − % pfparams.pmul - set the vectr flag with poles multiplicity (this option
+ − % is used only for continuous models)
+ − % pfparams.dterm - set the vector of direct terms
+ − % pfparams.fs - set the sampling frequency (Necessary for the discrete
+ − % case)
+ − %
+ − % OUTPUT:
+ − %
+ − % pfr is a struct containing output data and parameters
+ − % pfr.type = 'cont' if the model is continuous
+ − % pfr.type = 'disc' if the model is discrete
+ − % pfr.freq - frequencies vector
+ − % pfr.nfreq - normalized frequencies vector (Discrete case)
+ − % pfr.angfreq - angular frequencies vector
+ − % pfr.resp - frequency response data
+ − %
+ − %
+ − % VERSION: $Id: pfresp.m,v 1.5 2011/02/03 15:09:01 luigi Exp $
+ − %
+ − % HISTORY: 12-09-2008 L Ferraioli
+ − % Creation
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function pfr = pfresp(pfparams)
+ −
+ − %%% switching between continuous and discrete
+ −
+ − switch pfparams.type
+ − case 'cont'
+ − % collecting input parameters
+ − f = pfparams.freq;
+ − % willing to work with row
+ − if size(f,1)>size(f,2)
+ − f = f.';
+ − end
+ − r = pfparams.res;
+ − p = pfparams.pol;
+ − pmul = pfparams.pmul;
+ − d = pfparams.dterm;
+ − % willing to work with row
+ − if ~isempty(d) && size(d,1)>size(d,2)
+ − d = d.';
+ − end
+ −
+ − N = length(p);
+ −
+ − % substituted by faster code 03-Feb-2011
+ − % Nf = length(f);
+ − % rsp = zeros(Nf,1);
+ − % indx = (0:length(d)-1).';
+ − % for ii = 1:Nf
+ − % for jj = 1:N
+ − % m = pmul(jj);
+ − % rsptemp = r(jj)/(1i*2*pi*f(ii)-p(jj))^m;
+ − % rsp(ii) = rsp(ii) + rsptemp;
+ − % end
+ − % % Direct terms response
+ − % rsp(ii) = rsp(ii) + sum((((1i*2*pi*f(ii))*ones(length(d),1)).^indx).*d);
+ − % end
+ −
+ − % new code for a faster response calculation 03-Feb-2011
+ − rsp = zeros(size(f));
+ − for jj = 1:N
+ − m = pmul(jj);
+ − rsptemp = r(jj)./(1i*2*pi.*f-p(jj)).^m;
+ − rsp = rsp + rsptemp;
+ − end
+ − % get direct term response
+ − if ~isempty(d)
+ − Z = ones(numel(d),numel(f));
+ − ss = 2.*pi.*1i.*f;
+ − for jj=2:size(Z,1)
+ − Z(jj,:) = ss.^(jj-1);
+ − end
+ − rdtemp = d*Z;
+ − rsp = rsp + sum(rdtemp,1);
+ − end
+ −
+ − % Output
+ − pfr.type = 'cont';
+ − pfr.freq = f;
+ − pfr.angfreq = 2*pi*f;
+ − pfr.resp = rsp;
+ −
+ − case 'disc'
+ − % collecting input parameters
+ − f = pfparams.freq;
+ − fs = pfparams.fs;
+ − r = pfparams.res;
+ − p = pfparams.pol;
+ − d = pfparams.dterm;
+ −
+ − Nf = length(f);
+ − N = length(p);
+ −
+ − % Defining normalized frequencies
+ − fn = f./fs;
+ −
+ − rsp = zeros(Nf,1);
+ − indx = 0:length(d)-1;
+ − for ii = 1:Nf
+ − for jj = 1:N
+ − rsptemp = exp(1i*2*pi*fn(ii))*r(jj)/(exp(1i*2*pi*fn(ii))-p(jj));
+ − rsp(ii) = rsp(ii) + rsptemp;
+ − end
+ − % Direct terms response
+ − rsp(ii) = rsp(ii) + sum(((exp((1i*2*pi*f(ii))*ones(length(d),1))).^(-1.*indx)).*d);
+ − end
+ −
+ − % Output
+ − pfr.type = 'disc';
+ − pfr.freq = f;
+ − pfr.nfreq = fn;
+ − pfr.angfreq = 2*pi*f;
+ − pfr.resp = rsp;
+ − end
+ − end