% MSE curves for zoh + loudspeaker A, Volterra filter model up to P=4,
% memory N=75,90,120 and 150, 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
% sim_Azoh_r100_s25.m. The output was
% saved since the execution time is very long
load outputs_Azoh_r100_s25.mat outputs scales

% Literal linear-analytic model
[f,g,h,x,names,values]=dloudAitgr_lit;      

% State vector dimension
m=length(x);

% N=75
disp('N=75');
% Kernels (diagonal form)
kernels = voltkernels(f,g,h,x,names,values,P,75,1/Fs,1,1,1,0);
% Simulations for mse curves
reali=(1:100);
[mse75,times] = genmse(inputs,outputs,scales,kernels,'VF',P,reali,[],1);

%----

% N=90
disp('N=90');
% Kernels (diagonal form)
kernels = voltkernels(f,g,h,x,names,values,P,90,1/Fs,1,1,1,0);
% Simulations for mse curves
[mse90,times] = genmse(inputs,outputs,scales,kernels,'VF',P,reali,[],1);

%----

% N=120
disp('N=120');
% Kernels (diagonal form)
kernels = voltkernels(f,g,h,x,names,values,P,120,1/Fs,1,1,1,0);
% Simulations for mse curves
ini120=tic;
[mse120,times120] = genmse(inputs,outputs,scales,kernels,'VF',P,reali,[],1);
t120=toc(ini120);

%----

% N=150
disp('N=150');
% Kernels (diagonal form)
kernels = voltkernels(f,g,h,x,names,values,P,150,1/Fs,1,1,1,0);
% Simulations for mse curves
[mse150,times] = genmse(inputs,outputs,scales,kernels,'VF',P,reali,[],1);

%----

lsc=20*log10(scales);

for p=1:P
    lm{p}=10*log10([mse75(:,p) mse90(:,p) mse120(:,p) mse150(:,p)]);
end

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

%save mse_Azoh_P4_N75_150_r100_s25 mse75 mse90 mse120 mse150 scales times120 t120


