Mercurial > hg > ltpda
view m-toolbox/classes/+utils/@math/pfresp.m @ 11:9174aadb93a5 database-connection-manager
Add LTPDA Repository utility functions into utils.repository
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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