% Test ao/sineParams
%
% M Nofrarias 10-02-10
%
% $Id: test_ao_sineParams.m,v 1.2 2010/02/18 11:17:59 miquel Exp $

clear

% general parameters
fs  = 5;
nsecs =500;


%% one sinusoid 

% signal parameters
A1 =  1.23456;
f1 = 0.1;
sigma = 0.01;

% AOs
pl = plist('waveform','sine','f',f1,'A',A1,'fs',fs,'nsecs',nsecs);
a = ao(pl);
a.setName('signal 1')
pl = plist('waveform','noise','sigma',sigma,'fs',fs,'nsecs',nsecs);
c = ao(pl);
c.setName('noise')

% sum contributions
d = a  + c;
d.iplot

% get parameters
pl = plist('N',1);
p = sineParams(d,pl);

% compare with contributions
% bb = eval(p,plist('index',1));
% iplot(b,bb)

aa = eval(p,plist('index',1));
residual = aa-a;
residual.setName
iplot(aa,a,residual)


%% one sinusoid montecarlo

% signal parameters
A1 =  1.5;
f1 = 0.12;

for i = 1:100
fprintf('iteration %d\n',i);
sigma(i) = 0.01;

% AOs
pl = plist('waveform','sine','f',f1,'A',A1,'fs',fs,'nsecs',nsecs);
a = ao(pl);
a.setName('signal 1');
pl = plist('waveform','noise','sigma',sigma(i),'fs',fs,'nsecs',nsecs);
c = ao(pl);
c.setName('noise');

% sum contributions
d = a  + c;
% d.iplot

% get parameters
pl = plist('N',1);
p(i) = sineParams(d,pl);
A(i) = p(i).y(1);
dA(i) = p(i).dy(1);
f(i) = p(i).y(2);
df(i) = p(i).dy(2);
% compare with contributions
% bb = eval(p,plist('index',1));
% iplot(b,bb)

aa = eval(p(i),plist('index',1));
residual = aa-a;
residual.setName;
end

% compare mean(MSE) with computed Var + bias^2 
predicted_A = sqrt(std(A)^2 + (mean(A)-A1)^2)
computed_A = mean(dA)

predicted_f = sqrt(std(f)^2 + (mean(f)-f1)^2)
computed_f = mean(df)

%% two sinusoid

% signal parameters
A1 =  2.2;
f1 = 0.1;
A2 = 4.5;
f2 = 0.33;
sigma = 0.02;

% AOs
pl = plist('waveform','sine','f',f1,'A',A1,'fs',fs,'nsecs',nsecs);
a = ao(pl);
a.setName('signal 1')
pl = plist('waveform','sine','f',f2,'A',A2,'fs',fs,'nsecs',nsecs);
b = ao(pl);
b.setName('signal 2')
pl = plist('waveform','noise','sigma',sigma,'fs',fs,'nsecs',nsecs);
c = ao(pl);
c.setName('noise')

% sum contributions
d = a  + b + c;
% d.iplot

% get parameters
pl = plist('N',2);
p = sineParams(d,pl);

aa = eval(p,plist('index',1));


% compare with contributions
bb = eval(p,plist('index',1));
residual = bb-b;
residual.setName
iplot(bb,b,residual)

aa = eval(p,plist('index',2));
residual = aa-a;
residual.setName
iplot(aa,a,residual)