%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);