Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/pfresp.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % PFRESP returns frequency response of a partial fraction TF.
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % DESCRIPTION:
6 %
7 % Returns frequency response of a partial fraction expanded function
8 % (continuous or discrete).
9 %
10 % Continuous case
11 % The expected model is:
12 %
13 % r1 rN
14 % f(s) = ------ + ... + ------ + d1 + d2*s + ... + dK*s^{K-1}
15 % s - p1 s - p1
16 %
17 % Discrete case
18 % The expected model is:
19 %
20 % z*r1 z*rN
21 % f(z) = ------ + ... + ------ + d1 + d2*z^{-1} + ... + dK*z^{-(K-1)}
22 % z - p1 z - p1
23 %
24 % NOTE: The function cannot handle poles multiplicity higher than 1 in
25 % z domain. Multiple poles in s-domain are accepted.
26 %
27 % CALL: pfr = pfresp(pfparams)
28 %
29 % INPUT:
30 %
31 % pfparams is a struct containing input parameters
32 % pfparams.type = 'cont' Assumes a continuous model
33 % pfparams.type = 'disc' Assumes a discrete model
34 % pfparams.freq set the frequencies vector in Hz
35 % pfparams.res - set the vector of residues
36 % pfparams.pol - set the vector of poles
37 % pfparams.pmul - set the vectr flag with poles multiplicity (this option
38 % is used only for continuous models)
39 % pfparams.dterm - set the vector of direct terms
40 % pfparams.fs - set the sampling frequency (Necessary for the discrete
41 % case)
42 %
43 % OUTPUT:
44 %
45 % pfr is a struct containing output data and parameters
46 % pfr.type = 'cont' if the model is continuous
47 % pfr.type = 'disc' if the model is discrete
48 % pfr.freq - frequencies vector
49 % pfr.nfreq - normalized frequencies vector (Discrete case)
50 % pfr.angfreq - angular frequencies vector
51 % pfr.resp - frequency response data
52 %
53 %
54 % VERSION: $Id: pfresp.m,v 1.5 2011/02/03 15:09:01 luigi Exp $
55 %
56 % HISTORY: 12-09-2008 L Ferraioli
57 % Creation
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59
60 function pfr = pfresp(pfparams)
61
62 %%% switching between continuous and discrete
63
64 switch pfparams.type
65 case 'cont'
66 % collecting input parameters
67 f = pfparams.freq;
68 % willing to work with row
69 if size(f,1)>size(f,2)
70 f = f.';
71 end
72 r = pfparams.res;
73 p = pfparams.pol;
74 pmul = pfparams.pmul;
75 d = pfparams.dterm;
76 % willing to work with row
77 if ~isempty(d) && size(d,1)>size(d,2)
78 d = d.';
79 end
80
81 N = length(p);
82
83 % substituted by faster code 03-Feb-2011
84 % Nf = length(f);
85 % rsp = zeros(Nf,1);
86 % indx = (0:length(d)-1).';
87 % for ii = 1:Nf
88 % for jj = 1:N
89 % m = pmul(jj);
90 % rsptemp = r(jj)/(1i*2*pi*f(ii)-p(jj))^m;
91 % rsp(ii) = rsp(ii) + rsptemp;
92 % end
93 % % Direct terms response
94 % rsp(ii) = rsp(ii) + sum((((1i*2*pi*f(ii))*ones(length(d),1)).^indx).*d);
95 % end
96
97 % new code for a faster response calculation 03-Feb-2011
98 rsp = zeros(size(f));
99 for jj = 1:N
100 m = pmul(jj);
101 rsptemp = r(jj)./(1i*2*pi.*f-p(jj)).^m;
102 rsp = rsp + rsptemp;
103 end
104 % get direct term response
105 if ~isempty(d)
106 Z = ones(numel(d),numel(f));
107 ss = 2.*pi.*1i.*f;
108 for jj=2:size(Z,1)
109 Z(jj,:) = ss.^(jj-1);
110 end
111 rdtemp = d*Z;
112 rsp = rsp + sum(rdtemp,1);
113 end
114
115 % Output
116 pfr.type = 'cont';
117 pfr.freq = f;
118 pfr.angfreq = 2*pi*f;
119 pfr.resp = rsp;
120
121 case 'disc'
122 % collecting input parameters
123 f = pfparams.freq;
124 fs = pfparams.fs;
125 r = pfparams.res;
126 p = pfparams.pol;
127 d = pfparams.dterm;
128
129 Nf = length(f);
130 N = length(p);
131
132 % Defining normalized frequencies
133 fn = f./fs;
134
135 rsp = zeros(Nf,1);
136 indx = 0:length(d)-1;
137 for ii = 1:Nf
138 for jj = 1:N
139 rsptemp = exp(1i*2*pi*fn(ii))*r(jj)/(exp(1i*2*pi*fn(ii))-p(jj));
140 rsp(ii) = rsp(ii) + rsptemp;
141 end
142 % Direct terms response
143 rsp(ii) = rsp(ii) + sum(((exp((1i*2*pi*f(ii))*ones(length(d),1))).^(-1.*indx)).*d);
144 end
145
146 % Output
147 pfr.type = 'disc';
148 pfr.freq = f;
149 pfr.nfreq = fn;
150 pfr.angfreq = 2*pi*f;
151 pfr.resp = rsp;
152 end
153 end