line source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST Fitting procedure in 2dim z-domain VDFIT
% At the end you have 4 stable transfer functions in partial fractions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HISTORY: 02-10-2008 L Ferraioli
% Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VERSION
'$Id: test_2dim_vdfit.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);
csd11 = mpsd(:,2);
csd12 = mpsd(:,3);
csd21 = [];
csd22 = mpsd(:,4);
fs = 10;
%% Eigendecomposition
[tf11,tf12,tf21,tf22] = utils.math.eigcsd(csd11,csd12,csd21,csd22,'USESYM',0,'DIG',50,'OTP','TF');
%% Constructing vector
f1 = [tf11 tf21];
f2 = [tf12 tf22];
%% VDFIT
tic;
clear mlr1 mlr2
N = 32; %Order of approximation
% Max Iteration
Nmaxiter = 40;
% mlr1 = zeros(Nmaxiter,1);
% mlr2 = zeros(Nmaxiter,1);
% Fitting params
fitin.stable = 0;
fitin.dterm = 0;
fitin.plot = 0;
fitin.fs = fs;
weight = utils.math.wfun(f1,2);
% weight = 1./abs(f1);
pparams = struct('spolesopt',3, 'type','DISC', 'pamp', 0.98);
poles1 = utils.math.startpoles(N,f,pparams);
for hh = 1:Nmaxiter
[res1,poles1,dterm1,mresp1,rdl1,sqe1] = utils.math.vdfit(f1,f,poles1,weight,fitin); % Fitting
disp(['Iter' num2str(hh)])
% Stop condition checking
mlr1(hh) = sqe1(:,1);
mlr2(hh) = sqe1(:,2);
% [ext1,msg1] = utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,'lrsrmse',2,15);
% [ext2,msg2] = utils.math.stopfit(f1(:,2),rdl1(:,2),mlr2,'lrsrmse',2,15);
% order = length(poles1);
% mlr1(hh) = mean(log10(abs(f1(:,1)))-log10(abs(rdl1(:,1))));
% mlr2(hh) = mean(log10(abs(f1(:,2)))-log10(abs(rdl1(:,2))));
% [ext1,msg1] = utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,order,'lrs',1,0.0001,2);
% [ext2,msg2] = utils.math.stopfit(f1(:,2),rdl1(:,2),mlr2,order,'lrs',1,0.0001,2);
% if ext1 && ext2
% disp(msg1)
% break
% end
end
elpstime = toc;
% plotting squared error
figure()
semilogy(mlr1,'-ok')
hold on
grid on
semilogy(mlr2,'-or')
% % plotting RMS error variation
% figure()
% semilogy(abs(diff(sqrt(mlr1))),'-ok')
% hold on
% grid on
% semilogy(abs(diff(sqrt(mlr2))),'-or')
% plotting squared error variation
figure()
semilogy(abs(diff(mlr1)./mlr1(1,end-1)),'-ok')
hold on
grid on
semilogy(abs(diff(mlr2)./mlr2(1,end-1)),'-or')
%% All Passing
% [np1,resp1] = pfallpz(res1,poles1,dterm1,mresp1,f,fs);
[nr1,np1,nd1,resp1] = utils.math.pfallpsymz(res1,poles1,dterm1,mresp1,f,fs);
figure()
subplot(2,1,1);
p1 = loglog(f,abs(mresp1),'k');
hold on
p2 = loglog(f,abs(resp1),'r');
xlabel('Frequency [Hz]')
ylabel('Amplitude')
legend([p1(1) p2(1)],'VDFIT','Stabilized')
hold off
subplot(2,1,2);
p4 = semilogx(f,(180/pi).*unwrap(angle(mresp1)),'k');
hold on
p5 = semilogx(f,(180/pi).*unwrap(angle(resp1)),'r');
xlabel('Frequency [Hz]')
ylabel('Phase [Deg]')
legend([p4(1) p5(1)],'VDFIT', 'Stabilized')
hold off
%% VDFIT
clear mlr3 mlr4 ext1 ext2
N = 24; %Order of approximation
% Max Iteration
Nmaxiter = 70;
% Fitting params
fitin.stable = 0;
fitin.dterm = 0;
fitin.plot = 1;
fitin.fs = fs;
weight = utils.math.wfun(f2,2);
pparams = struct('spolesopt',2, 'type','DISC', 'pamp', 0.98);
poles2 = utils.math.startpoles(N,f,pparams);
for hh = 1:Nmaxiter
[res2,poles2,dterm2,mresp2,rdl2,rmse2] = utils.math.vdfit(f2,f,poles2,weight,fitin); % Fitting
disp(['Iter' num2str(hh)])
% Stop condition checking
mlr3(hh) = rmse2(:,1);
mlr4(hh) = rmse2(:,2);
[ext1,msg1] = utils.math.stopfit(f2(:,1),rdl2(:,1),mlr3,'lrsrmse',2,15);
[ext2,msg2] = utils.math.stopfit(f2(:,2),rdl2(:,2),mlr4,'lrsrmse',2,15);
if ext1 && ext2
disp(msg1)
break
end
end
% plotting RMS error
figure()
semilogy(mlr3,'-ok')
hold on
grid on
semilogy(mlr4,'-or')
%%
% params = struct('spolesopt',2, 'Nmaxiter',250, 'minorder',30,...
% 'maxorder',60, 'weightparam',1, 'plot',0,...
% 'ctp','lrs','lrscond',1,'lrsvarcond',1,'nrmsecond',2,...
% 'stabfit',0,'dterm',0,'spy',0);
%
% [res2,poles2,dterm2,mresp2,rdl2] = autodfit(f2,f,fs,params);