% MSE curves for zoh + loudspeaker A, Volterra-CPD model up to P=4,
% memory N=100, rank R=3,5,10 and 15, 100 realizations and 25 scales of input
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

clear

% Order
P=4;
% Sampling frequency
Fs=5000;

%----------------------------------------------------------------
load inputs_r100.mat inputs

% Simulation output generated with ode133 by
% cbvkdemos/PhysicalModels/LoudspeakerA/sim_Azoh_r100_s25.m. The output was
% saved since the execution time is very long
load outputs_Azoh_r100_s25.mat outputs scales
reali=(1:100);


% Volterra-CPD models generated with Tensorlab 3.0

% R=3
load vcpdmodel_Azoh_P4_N100_R3.mat model
[mse3,times] = genmse(inputs,outputs,scales,model,'CPD',P,reali,[],0);

% R=5
load vcpdmodel_Azoh_P4_N100_R5.mat model
[mse5,times] = genmse(inputs,outputs,scales,model,'CPD',P,reali,[],0);

% R=10
load vcpdmodel_Azoh_P4_N100_R10.mat model
[mse10,times] = genmse(inputs,outputs,scales,model,'CPD',P,reali,[],0);

% R=15
load vcpdmodel_Azoh_P4_N100_R15.mat model
init15=tic;
[mse15,times15] = genmse(inputs,outputs,scales,model,'CPD',P,reali,[],0);
t15=toc(init15);

%----

%save mse_Azoh_P4_N100_R3_15_r100_s25 mse3 mse5 mse10 mse15 scales times15 t15

lsc=20*log10(scales);

for p=1:P
    lm{p}=10*log10([mse3(:,p) mse5(:,p) mse10(:,p) mse15(:,p)]);
end

plot(lsc,lm{1},'r',lsc,lm{2},'m',lsc,lm{3},'b',lsc,lm{4},'g')



