Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test_psdzfit.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST Fitting procedure in 1dim PSDZFIT
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 02-10-2008 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test_psdzfit.m,v 1.1 2009/05/08 13:48:02 luigi Exp $';
11
12 %% Clear
13
14 clear all
15
16 %% Loading spectra
17
18 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
19
20 f = mpsd(:,1);
21 psd = mpsd(:,2);
22 % csd12 = mpsd(:,3);
23 % csd21 = [];
24 % csd22 = mpsd(:,4);
25
26 fs = 10;
27
28 %% PSDFIT
29 tic;
30
31 clear mlr1 mlr2
32
33 N = 11; %Order of approximation
34
35
36 % Max Iteration
37 Nmaxiter = 40;
38
39 % mlr1 = zeros(Nmaxiter,1);
40 % mlr2 = zeros(Nmaxiter,1);
41
42 % Fitting params
43 fitin.plot = 1;
44 fitin.fs = fs;
45
46 weight = utils.math.wfun(psd,2);
47 % weight = 1./abs(f1);
48
49 clear res poles dterm mresp rdl sqe mlr
50
51 pparams = struct('spolesopt',2, 'type','DISC', 'pamp', 0.98);
52 poles = utils.math.startpoles(N,f,pparams);
53
54 for hh = 1:Nmaxiter
55 [res,poles,fullpoles,mresp,rdl,sqe] = utils.math.psdzfit(psd,f,poles,weight,fitin); % Fitting
56 % disp(['Iter' num2str(hh)])
57
58 % Stop condition checking
59 mlr(hh) = sqe(:,1);
60 % [ext1,msg1] =utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,'lrsrmse',2,15);
61
62 % if ext1 && ext2
63 % disp(msg1)
64 % break
65 % end
66
67 end
68 elpstime = toc;
69
70 % remove NaN from mlr
71 idnan = isnan(mlr);
72 mlr(idnan) = 0;
73
74 % plotting squared error
75 figure()
76 semilogy(mlr,'-ok')
77 legend('MSE')
78
79 % plotting squared error variation
80 figure()
81 semilogy(abs(diff(mlr)./mlr(1,end-1)),'-ok')
82 legend('MSE Variation')
83
84
85
86 %%
87
88 [num,den] = residue(res,[poles;1./poles],0);
89 zrs = roots(num);
90