Mercurial > hg > ltpda

comparison m-toolbox/test/diagnostics/ltpda_arma_freq.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 function varargout = ltpda_arma_freq(varargin)
2 % LTPDA_ARMA_FREQ estimates the ARMA parameters of the transfer function
3 % relating an output to an input
4 %
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %
7 % DESCRIPTION: LTPDA_ARMA_FREQ estimates the ARMA parameters of the
8 % transfer function relating an output to an input time series. The
9 % algorithm uses the IV method to find a set of initial parameters
10 % and then performs an iterative search based on the Optimization
11 % Toolbox.
12 %
13 % CALL: b = ltpda_arma_freq(ax,ay,pl)
14 %
15 % INPUTS: ax - analysis object containig the input time series
16 % ay - analysis object containig the output time series
17 % pl - parameters list
18 %
19 % OUTPUTS: b - cdata type analysis object containing the ARMA
20 % parameters and the residual of the fit
21 %
22 % PARAMETERS:
23 % MaxIter - Maximum number of iterations to be performed by the
24 % iterative search (default 4e2)
25 % MaxFunEvals - Maximum number of function evaluations (default 1e2)
26 % TolX - Tolerance on the estimated parameter (default 1e-6)
27 % TolFun - Tolerance on the evaluated function (default 1e-6)
28 % UpBound - Array of parameters upper bounds for the iterative search
29 % (default is 1 for each parameter)
30 % LowBound - Array of parameters lower bounds for the iterative search
31 % (default is 1 for each parameter)
32 % ARMANum - MA order of the ARMA filter (default 2)
33 % ARMADen - AR order of the ARMA filter (default 1)
34 %
35 % VERSION: $Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
36 %
37 % HISTORY: 25-02-2008 M Nofrarias
38 % Creation
39 %
40 %
41 % TODO: - Add parameters errors
42 % - Pass from univariate to multivariate
43 %
44 %
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46
47
48
49 %% Standard history variables
50
51 ALGONAME = mfilename;
52 VERSION = '$Id: ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $';
53 CATEGORY = 'SIGPROC';
54
55
56 %% Check if this is a call for parameters
57 if nargin == 1 && ischar(varargin{1})
58 in = char(varargin{1});
59 if strcmpi(in, 'Params')
60 varargout{1} = getDefaultPL();
61 return
62 elseif strcmpi(in, 'Version')
63 varargout{1} = VERSION;
64 return
65 elseif strcmpi(in, 'Category')
66 varargout{1} = CATEGORY;
67 return
68 end
69 end
70
71 %% Capture input variables names
72 invars = {};
73 as = [];
74 ps = [];
75 for j=1:nargin
76 invars = [invars cellstr(inputname(j))];
77 if isa(varargin{j}, 'ao')
78 as = [as varargin{j}];
79 end
80 if isa(varargin{j}, 'plist')
81 ps = [ps varargin{j}];
82 end
83 end
84
85
86 %% Check plist
87 if isempty(ps)
88 pl = getDefaultPL();
89 else
90 pl = combine(ps, getDefaultPL);
91 end
92
93
94 %% Fourier transform input varialbles
95
96 Fs=0.65;
97 L=length(xin)
98 NFFT = 2^nextpow2(L); % Next power of 2 from length of y
99 NFFT=length(xin);
100
101 x_ = fft(xin,NFFT)/L;
102 y_ = fft(yout',NFFT)/L;
103
104 Px = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
105 Py = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.
106
107 f = Fs/2*linspace(0,1,NFFT/2)';
108
109
110
111
112 %% Iterative search using the Optimization Toolbox
113
114
115 % Options for the lsqcurvefit function
116 opt = optimset(...
117 'MaxIter',find(pl, 'MaxIter'),...
118 'TolX',find(pl, 'TolXm'),...
119 'TolFun',find(pl, 'TolFun'),...
120 'MaxFunEvals',find(pl, 'MaxFunEvals'),...
121 'Display','off');
122
123 % Upper and Lower Bounds for the parameters
124 ub=find(pl,'UpBound');
125 lb=find(pl,'LowBound');
126 if isempty(ub)
127 ub=ones(p+q,1);
128 pl = append(pl, param('UpBound', ub));
129 disp('! Parameters Upper Bound set to 1')
130 end
131 if isempty(lb)
132 lb=-ones(p+q,1);
133 pl = append(pl, param('LowBound', lb));
134 disp('! Parameters Lower Bound set to -1')
135 end
136
137 % Call to the Optimization Toolbox function
138 [par,res]=lsqcurvefit(@(par,xdata)filter([par(1:p)],[1 par(p+1:p+q)],xdata),...
139 par0,x,y,lb,ub,opt);
140
141
142 %% Build output ao
143
144 % New output history
145 h = history(ALGONAME, VERSION, pl,[as(1).hist as(2).hist]);
146 h = set(h,'invars', invars);
147
148 % Make output analysis object
149 b = ao([par; res]);
150
151 % Name, mfilename, description for this object
152 b = setnh(b,...
153 'name', ['ARMA param. Input: ' sprintf('%s ', char(invars{1})) ' Output: ' char(invars{2})],...
154 'mfilename', ALGONAME, ...
155 'description', find(pl,'description'),...
156 'hist', h);
157
158 varargout{1} = b;
159
160
161 %% --- FUNCTIONS ---
162
163
164 % Get default parameters
165 function plo = getDefaultPL();
166 disp('* creating default plist...');
167 plo = plist();
168 plo = append(plo, param('MaxIter',4e2));
169 plo = append(plo, param('MaxFunEvals',1e2));
170 plo = append(plo, param('TolX',1e-6));
171 plo = append(plo, param('TolFun',1e-6));
172 plo = append(plo, param('ARMANum',2));
173 plo = append(plo, param('ARMADen',1));
174 plo = append(plo, param('fs', 1));
175 disp('* done.');
176
177
178
179 % siz=length(Px);
180 % % siz=200;
181 %
182 % xdataf=Px(1:siz)*1e5;
183 % ydataf=Py(1:siz)*1e5;
184 % % xdataf=Px.data.y;
185 % % ydataf=Py.data.y;
186 % global freq
187 % % freq=Px.data.x;
188 % freq =f(1:siz);
189
190
191 tic;
192 [x1f,resnorm1f]=...
193 lsqcurvefit(@myfunf,p0,xdataf,ydataf,[-1 -1 -1],[1 1 1],options);
194 elaps=toc
195
196 % %
197 % %%
198 % % Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
199 % % F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
200 % %
201 % % figure(1)
202 % % hold off
203 % % loglog(Px.data.x,ydataf./xdataf)
204 % % hold
205 % % loglog(Px.data.x,abs(Ff),'k')
206 % % loglog(Px.data.x,abs(F),'m')
207 % % loglog(f,abs((39.6e-3-39.5e-3*exp(-i*2*pi*f/0.65))./(1-0.996*exp(-i*2*pi*f/0.65))),'r')
208 % %
209 % %
210 %
211 % %%
212 %
213 % Fs=18;
214 % Lw=2;
215 %
216 % Ff = ((x1f(1)+x1f(2)*exp(-i*2*pi*f/0.65))./(1+x1f(3)*exp(-i*2*pi*f/0.65)));
217 % F = ((x1(1)+x1(2)*exp(-i*2*pi*f/0.65))./(1+x1(3)*exp(-i*2*pi*f/0.65)));
218 %
219 % figure(1)
220 % set(gca,'Fontsize',Fs)
221 % hold off
222 % loglog(f,Py./Px,'LineWidth',Lw)
223 % hold
224 % % % loglog(f,ydataf,'m')
225 % loglog(f,abs(Ff),'k','LineWidth',Lw)
226 % loglog(f,abs(F),'m','LineWidth',Lw)
227 %
228 % % loglog(f,(x1(1)+x1(2)*exp(-i*2*pi*f/1.32))./(1+x1(3)*exp(-i*2*pi*f/1.32)),'r')
229 % loglog(f,abs((39.6e-3-39.5e-3*exp(-i*2*pi*f/0.65))./(1-0.996*exp(-i*2*pi*f/0.65))),'r')
230 %
231 % legend('Empirical','Freq.domain','Time domain','Analytic')
232 %