% A test script for ao/noisegen1D
% 
% L. Ferraioli 10-11-08
% 
% $Id: test_ao_noisegen1D.m,v 1.7 2010/05/04 13:30:35 luigi Exp $
% 

%% Make white noise and compare results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Description:
% 1) Generates a random series of data (white)
% 2) Build a pzm
% 3) Calculate equivalent one sided data psd from pzm response
% 4) Generates a miir object from pzm
% 5) Filer white data with miir to generate reference colored noise
% 6) Generates colored noise with noisegen1D
% 7) calculated psd of generated and reference data
% 8) check result by plotting

% 1)
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));

% PSD model
% 2)
pzm = pzmodel(1, {0.01}, {0.1});
% 3)
mod = (abs(pzm.resp).^2).*(2/fs); % this corresponds to the theoretical one sided psd of the data
% 4)
fpzm = miir(pzm,plist('fs',fs));
% 5)
acr = filter(a,fpzm);

% 6) Noise generation
pl = plist(...
    'model', mod, ...
    'MaxIter', 10, ...
    'PoleType', 3, ...
    'MinOrder', 2, ...
    'MaxOrder', 9, ...
    'Weights', 2, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-3,...
    'FITTOL', 1e-5);

ac = noisegen1D(a, pl);

% 7)
acxx = ac.psd;
acrxx = acr.psd;
% 8)
iplot(acxx,acrxx,mod);

%% Noise generation from fsdata model object - output time series %%%%%%%%

% Description:
% 1) Generates a fsdata object to be used as psd model
% 2) Generates a random series of data (white)
% 3) Generates colored noise with noisegen1D
% 4) calculated psd of generated data
% 5) check result by plotting

% 1)
fs = 10; % sampling frequency
f = logspace(-6,log10(5),1000);
pl_mod1 = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
mod1 = ao(pl_mod1); % fsdata model object

% 2)
% generating white noise
a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000));

% 3) Noise generation

pl1 = plist(...
    'model', mod1, ...
    'MaxIter', 30, ...
    'PoleType', 2, ...
    'MinOrder', 3, ...
    'MaxOrder', 20, ...
    'Weights', 3, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-3,...
    'FITTOL', 1e-5);

ac1 = noisegen1D(a1, pl1);

% 4)
acxx1 = ac1.psd;
% 5)
iplot(acxx1, mod1);

%% Noise generation from fsdata model object - output filter %%%%%%%%%%%%%

% Description:
% 1) Generates a fsdata object to be used as psd model
% 2) Generates coloring filter with noisegen1D
% 3) Generates white noise data
% 4) filter data
% 4) calculated psd of filtered data
% 5) check result by plotting

% 1) Generates a fsdata object to be used as psd model
fs = 10; % sampling frequency
f = logspace(-6,log10(5),1000);
pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
mod = ao(pl_mod); % fsdata model object

% 2) noise coloring filter generation

pl1 = plist(...
    'fs', fs, ...
    'Iunits', '', ...
    'Ounits', 'm', ....
    'MaxIter', 30, ...
    'PoleType', 2, ...
    'MinOrder', 3, ...
    'MaxOrder', 20, ...
    'Weights', 3, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-3,...
    'FITTOL', 1e-5);

fil = noisegen1D(mod, pl1);

% 3) generating white noise
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));

% 4) filter white noise
b = filter(a,fil);

% 4) psd estimation
bxx = b.psd(plist('navs',8,'olap',50,'order',1));
% 5)
iplot(bxx, mod);

%% Noise generation from fsdata model object - output filter 2 %%%%%%%%%%%%%

% Description:
% 1) Generates a fsdata object to be used as psd model
% 2) Generates coloring filter with noisegen1D
% 3) Generates white noise data
% 4) filter data
% 4) calculated psd of filtered data
% 5) check result by plotting

% 1) Generates a fsdata object to be used as psd model
fs = 10; % sampling frequency
f = logspace(-6,log10(5),1000);
pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
mod = ao(pl_mod); % fsdata model object

% 2) noise coloring filter generation

pl1 = plist(...
    'fs', fs, ...
    'Iunits', '', ...
    'Ounits', 'm', ....
    'MaxIter', 30, ...
    'PoleType', 2, ...
    'MinOrder', 3, ...
    'MaxOrder', 20, ...
    'Weights', 3, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-3,...
    'FITTOL', 1e-5);

fil = noisegen1D(mod, mod, pl1);

% 3) generating white noise
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));

% 4) filter white noise
b = filter(a,fil(1));

% 4) psd estimation
bxx = b.psd(plist('navs',8,'olap',50,'order',1));
% 5)
iplot(bxx, mod);

%% Noise generation from fsdata model object - output filter 3 %%%%%%%%%%%%%

% Description:
% 1) Generates a fsdata object to be used as psd model
% 2) Generates coloring filter with noisegen1D
% 3) Generates white noise data
% 4) filter data
% 4) calculated psd of filtered data
% 5) check result by plotting

% 1) Generates a fsdata object to be used as psd model
fs = 10; % sampling frequency
f = logspace(-6,log10(5),1000);
pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
mod = ao(pl_mod); % fsdata model object

% 2) noise coloring filter generation

pl1 = plist(...
    'fs', fs, ...
    'Iunits', '', ...
    'Ounits', 'm', ....
    'MaxIter', 30, ...
    'PoleType', 2, ...
    'MinOrder', 3, ...
    'MaxOrder', 20, ...
    'Weights', 3, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-3,...
    'FITTOL', 1e-5);

[fil1, fil2] = noisegen1D(mod, mod, pl1);

% 3) generating white noise
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));

% 4) filter white noise
b = filter(a,fil1);

% 4) psd estimation
bxx = b.psd(plist('navs',8,'olap',50,'order',1));
% 5)
iplot(bxx, mod);

% END