% A test script for utils.math.blwhitenoise
% 
% L. Ferraioli 18-11-08
% 
% $Id: test_eigcsd.m,v 1.2 2009/06/10 16:12:06 luigi Exp $
% 

%%

clear all

%%

% sampling frequency
fs = 10;
% higher bandwhidth frequency
fh = fs/2;
% lower bandwidth frequency
fl = 1e-3;
% number of data points
npts = 1e6;

% noise generation
wn = utils.math.blwhitenoise(npts,fs,fl,fh);

wnao = ao(tsdata(wn,fs));
% wnao.setNane('BLWN');

%% plot and statistics


iplot(wnao)

% Check histogram
wnao_hist = hist(wnao,plist('N',500,'Norm','true'));
wnao_norm = normdist(wnao,plist('N',500,'Norm','true'));

plpl = plist( ...
  'Yscales', {'All','log'}, ...
  'Legends', {'Histogram', 'Gaussian'});

iplot(wnao_hist,wnao_norm,plpl)

%% fft

WNf = fft(wnao);
iplot(WNf)

%% psd

% expected white noise at fs
f = [logspace(-5,fs/2,300)]';
expec = ao(plist('fsfcn', 'ones(size(f))', 'f', f));
expec = expec.*(2/fs);
WNxx = psd(wnao,plist('order',-1,'nfft',1e4));
iplot(WNxx,expec)

%% compare randn

aa = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));

iplot(wnao)

% Check histogram
aa_hist = hist(aa,plist('N',500,'Norm','true'));
aa_norm = normdist(aa,plist('N',500,'Norm','true'));

plpl = plist( ...
  'Yscales', {'All','log'}, ...
  'Legends', {'Histogram', 'Gaussian'});

iplot(aa_hist,aa_norm,plpl)