--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/utils/test_psdzfit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,90 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST Fitting procedure in 1dim PSDZFIT
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     02-10-2008 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% VERSION
+
+'$Id: test_psdzfit.m,v 1.1 2009/05/08 13:48:02 luigi Exp $';
+
+%% Clear
+
+clear all
+
+%% Loading spectra
+
+load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
+
+f = mpsd(:,1);
+psd = mpsd(:,2);
+% csd12 = mpsd(:,3);
+% csd21 = [];
+% csd22 = mpsd(:,4);
+
+fs = 10;
+
+%% PSDFIT
+tic;
+
+clear mlr1 mlr2
+
+N = 11; %Order of approximation 
+
+
+% Max Iteration
+Nmaxiter = 40;
+
+% mlr1 = zeros(Nmaxiter,1);
+% mlr2 = zeros(Nmaxiter,1);
+
+% Fitting params
+fitin.plot = 1;
+fitin.fs = fs;
+
+weight = utils.math.wfun(psd,2);
+% weight = 1./abs(f1);
+
+clear res poles dterm mresp rdl sqe mlr
+
+pparams = struct('spolesopt',2, 'type','DISC', 'pamp', 0.98);
+poles = utils.math.startpoles(N,f,pparams);
+
+for hh = 1:Nmaxiter
+  [res,poles,fullpoles,mresp,rdl,sqe] = utils.math.psdzfit(psd,f,poles,weight,fitin); % Fitting
+%   disp(['Iter' num2str(hh)])
+
+  % Stop condition checking
+  mlr(hh) = sqe(:,1);
+%   [ext1,msg1] =utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,'lrsrmse',2,15);
+
+%   if ext1 && ext2
+%     disp(msg1)
+%     break
+%   end
+
+end
+elpstime = toc;
+
+% remove NaN from mlr
+idnan = isnan(mlr);
+mlr(idnan) = 0;
+
+% plotting squared error
+figure()
+semilogy(mlr,'-ok')
+legend('MSE')
+
+% plotting squared error variation
+figure()
+semilogy(abs(diff(mlr)./mlr(1,end-1)),'-ok')
+legend('MSE Variation')
+
+
+
+%%
+
+[num,den] = residue(res,[poles;1./poles],0);
+zrs = roots(num);
+