%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST csd2tf 1dim identification
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HISTORY:     23-04-2009 L Ferraioli
%                 Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VERSION

'$Id: test1_csd2tf.m,v 1.1 2009/04/23 10:11:26 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);

%% identification in z-domain

% identification params
params = struct('Nmaxiter',50, 'minorder',10,...
    'maxorder',25, 'spolesopt',2, 'weightparam',2, 'plot',0,...
    'ctp','chival','lrscond',3,'msevar',2,...
    'fs',10,'dterm',0, 'spy',0, 'fullauto',1,...
    'extweights', []);

ostruct = utils.math.csd2tf(psd,f,params);

%% Model response

res = ostruct.res;
poles = ostruct.poles;
dterm = ostruct.dterm;

pfparams.type = 'disc';
pfparams.freq = f;
pfparams.fs = 10;
pfparams.res = res;
pfparams.pol = poles;
pfparams.dterm = dterm;

% response of the model filter
pfr = utils.math.pfresp(pfparams);
y = pfr.resp;

%% Checking result

figure()
loglog(f,psd,'k')
hold on
loglog(f,y.*conj(y),'r')
xlabel('Frequency [Hz]')
ylabel('psd')
legend('Original', 'Model')

figure()
loglog(f,abs(y),'k')
hold on
xlabel('Frequency [Hz]')
ylabel('Amplitude')
legend('TF Model')

figure()
semilogx(f,(180/pi).*unwrap(angle(y)),'k')
hold on
xlabel('Frequency [Hz]')
ylabel('Phase Angle [Deg]')
legend('TF Model')

figure()
semilogy(ostruct.mse)
xlabel('Iteration')
ylabel('MSE')
legend('Mean Squared Error Progression')