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