clear all

%% Make test AO
fs    = 10;
Nsecs = 10;
a = ao(plist('waveform', 'sine wave', 'A', 1, 'f', .3, 'phi', 0, 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
b = ao(plist('waveform', 'sine wave', 'A', 2*pi*0.3, 'f', .3, 'phi', 90, 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm')); % Theoretical first derivative
b2 = -1*((2*pi*0.3)^2).*a; % Theoretical second derivative

%% Test first derivative
c = diff(a, plist('method', '2POINT'));
d = diff(a, plist('method', 'ORDER2'));
e = diff(a, plist('method', 'ORDER2SMOOTH'));
f = diff(a, plist('method', '3POINT'));
g = diff(a, plist('method', '5POINT'));
h = diff(a, plist('method', 'FPS', 'ORDER', 'FIRST', 'COEFF', -1/5));

pl = plist('Legends', {'Original','2POINT','ORDER2','ORDER2SMOOTH','3POINT','5POINT','FPS PAR','Theory'},...
  'Markers', {'x', 's', '^', 'p', '+', '<', 'o',''});
iplot(a,c,d,e,f,g,h,b, pl)

%% Use noise and measure TFs for First Derivatives

Nsecs = 1000;

a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
b = diff(a, plist('method', '2POINT'));
c = diff(a, plist('method', 'ORDER2'));
d = diff(a, plist('method', 'ORDER2SMOOTH'));
e = diff(a, plist('method', '3POINT'));
f = diff(a, plist('method', '5POINT'));
g = diff(a, plist('method', 'FPS', 'ORDER', 'FIRST', 'COEFF', -1/5));
h = diff(a, plist('method', 'FPS', 'ORDER', 'FIRST', 'COEFF', 1/12));

pl  = plist('Nfft', fs*100);
tf1 = tfe(a,b,pl);
tf2 = tfe(a,c,pl);
tf3 = tfe(a,d,pl);
tf4 = tfe(a,e,pl);
tf5 = tfe(a,f,pl);
tf6 = tfe(a,g,pl);
tf7 = tfe(a,h,pl);

freq = tf1(1,2).data.x;
tfthy = 2.*pi.*1i.*freq;
tfth = ao(fsdata(freq,tfthy));

pl = plist('Legends', {'2POINT','ORDER2','ORDER2SMOOTH','3POINT',...
  '5POINT', 'FPS PAR', 'FPS SER','Theory',});
iplot(tf1(1,2), tf2(1,2), tf3(1,2), tf4(1,2), tf5(1,2), tf6(1,2), tf7(1,2), tfth, pl)

%% Make test AO for second derivative

fs    = 10;
Nsecs = 10;
a = ao(plist('waveform', 'sine wave', 'A', 1, 'f', .3, 'phi', 0, 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
b = ao(plist('waveform', 'sine wave', 'A', 2*pi*0.3, 'f', .3, 'phi', 90, 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm')); % Theoretical first derivative
b2 = -1*((2*pi*0.3)^2).*a; % Theoretical second derivative

%% Test second derivatives

c = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', 2/7));
d = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', -1/12));
e = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', 1/4));

pl = plist('Legends', {'Original','FPS PAR','FPS SER','FPS PI','Theory'},...
  'Markers', {'x', 'o', 's', '^', ''});
iplot(a,c,d,e,b2, pl)

%% Use noise and measure TFs for SECOND Derivatives

Nsecs = 1000;

a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
c = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', 2/7));
d = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', -1/12));
e = diff(a, plist('method', 'FPS', 'ORDER', 'SECOND', 'COEFF', 1/4));

pl  = plist('Nfft', fs*100);

tf1 = tfe(a,c,pl);
tf2 = tfe(a,d,pl);
tf3 = tfe(a,e,pl);

freq = tf1(1,2).data.x;
tfthy = abs((2.*pi.*1i.*freq).^2);
tfth = ao(fsdata(freq,tfthy));

pl = plist('Legends', {'Theory','FPS PAR', 'FPS PI', 'FPS SER',});
iplot(tfth, tf1(1,2), tf3(1,2), tf2(1,2), pl)