% util script to generate training session data for topic 5
% L FERRAIOLI 22-02-09
%
% $Id: gendata.m,v 1.4 2009/10/02 08:31:21 mauro Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% save path
data_path = '/Users/MH/Matlab/ltpda_cvs/software/m-toolbox/test/LTPDA_Training/topic5/';

%% 1) model filter

% Model Residues and Poles
mRes = [2.44554138162509e-011 - 1.79482547894083e-011i;
2.44554138162509e-011 + 1.79482547894083e-011i;
2.66402334803101e-009 +  1.1025122049153e-009i;
2.66402334803101e-009 -  1.1025122049153e-009i;
-7.3560293387644e-009;
-1.82811618589835e-009 - 1.21803627800855e-009i;
-1.82811618589835e-009 + 1.21803627800855e-009i;
1.16258677367555e-009;
1.65216557639319e-016;                         
-1.78092396888606e-016;
-2.80420398962379e-017;
9.21305973049041e-013 - 8.24686706827269e-014i;
9.21305973049041e-013 + 8.24686706827269e-014i;
5.10730060739905e-010 - 3.76571756625722e-011i;
5.10730060739905e-010 + 3.76571756625722e-011i;
3.45893698149735e-009;
3.98139182134446e-014 - 8.25503935419059e-014i;
3.98139182134446e-014 + 8.25503935419059e-014i;
-1.40595719147164e-011];

mPoles = [0.843464045655194 - 0.0959986292915475i;
0.843464045655194 + 0.0959986292915475i;
0.953187595424927 - 0.0190043625473383i;
0.953187595424927 + 0.0190043625473383i;
0.967176277937188;
0.995012027005247 - 0.00268322602801729i;
0.995012027005247 + 0.00268322602801729i;
0.996564761885673;
0.999999366165445;
0.999981722418555;
0.999921882627659;
0.999624431675213 - 0.000813407848742761i;
0.999624431675213 + 0.000813407848742761i;
0.997312006278751 - 0.00265611346834941i;
0.997312006278751 + 0.00265611346834941i;
0.990516544257531;
0.477796923118318 - 0.311064085401834i;
0.477796923118318 + 0.311064085401834i;
0];

f = logspace(-6,log10(5),100).';
fs = 10;

pfparams.type = 'disc';
pfparams.freq = f;
pfparams.fs = fs;
pfparams.res = mRes;
pfparams.pol = mPoles;
pfparams.dterm = 0;

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

% saving AO
rfilt = ao(plist('xvals', f, 'yvals', y, 'fs', fs, 'dtype', 'fsdata'));
rfilt.setName;

%
rfilt.save(plist('filename',[data_path 'T5_Ex02_rfilt.mat']));
rfilt.save(plist('filename',[data_path 'T5_Ex02_rfilt.xml']));

%% 2) Peaked model psd

% Create a psd model with two peak resonance at 5e-2 and 7e-2 Hz
func = '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e6.*f.^2).*1e-10 + 1e-5./(1e-3+((f./1e-2).^2-1).^2) + 1e-6./(1e-3+((f./2e-2).^2-1).^2)';
pl_data = plist('fsfcn', func, 'f1', 1e-6, 'f2', .5, 'nf', 300);
mod = ao(pl_data);
iplot(mod)

%% Building white noise

a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000, 'yunits', 'm'));
a.setName;

%% Calling the noise generator

pl = plist(...
    'model', mod, ... % Multiplication by fs needed to preserve energy
    'MaxIter', 70, ...
    'PoleType', 2, ...
    'MinOrder', 10, ...
    'MaxOrder', 45, ...
    'Weights', 2, ...
    'Plot', false,...
    'Disp', false,...
    'RMSEVar', 7,...
    'FitTolerance', 2);

ac = noisegen1D(a, pl);

%% 

acxx = ac.psd(plist('Nfft',2000));
iplot(acxx,mod)

%%

ac.save(plist('filename',sprintf([data_path 'T5_Ex03_TestNoise.mat'])));
ac.save(plist('filename',sprintf([data_path 'T5_Ex03_TestNoise.xml'])));

%% 3) Generate polynomial tsdata

% Make fake AO from polyval
nsecs = 1000;
fs    = 1;
  
a1 = ao(plist('tsfcn', 'polyval([1e-11 0 -5e-6 -1e-3 -5e-1 0.5], t) + 5e2*randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));

% iplot(a1)

a1.save(plist('filename',sprintf([data_path 'T5_Ex04_TestNoise.mat'])));
a1.save(plist('filename',sprintf([data_path 'T5_Ex04_TestNoise.xml'])));

%% 4) Generate tsdata from function

nsecs = 1000;
fs = 1;
a2 = ao(plist('tsfcn','5 + 3.*sin(2.*pi.*(1e-4 + 1e-5.*t).*t + 0.3) + randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));

iplot(a2)

a2.save(plist('filename',sprintf([data_path 'T5_Ex05_TestNoise.mat'])));
a2.save(plist('filename',sprintf([data_path 'T5_Ex05_TestNoise.xml'])));

%% 5) Temp and Ifo data

% read ao
ifo_rw = ao('ifo_temp_example/ifo_training.dat');
ifo_rw.setYunits('rad');
T_rw = ao('ifo_temp_example/temp_training.dat');
T_rw.setYunits('degC');

%% 

% Consolidate data
[ifo_c,T_c] = consolidate(ifo_rw,T_rw,plist('fs',1));
ifo_c.save(plist('filename',sprintf([data_path 'ifo_Ex06.mat'])));
ifo_c.save(plist('filename',sprintf([data_path 'ifo_Ex06.xml'])));
T_c.save(plist('filename',sprintf([data_path 'T_Ex06.mat'])));
T_c.save(plist('filename',sprintf([data_path 'T_Ex06.xml'])));

%%

% do TF
% pl2 = plist('Nfft', 2000);
% tf = tfe(T_c,ifo_c,pl2);
tf = tfe(T_c,ifo_c);