% test TF code on simple example
%
% DGM 1/6/11
%

% set parameters of model:
% NOTE: gest currently returns mu/tau and 1/tau as "mu" and "epsilon"
mu=1/4;
tau=1/4;          % time constant in years    
% ODE is \tau\dot{T}=mu*H-T+n, or \dot{T}=(mu/tau)H-(1/tau)T+(1/tau)n
% or in Laplace, T/H = (mu/tau)*(s+1/tau)^(-1) = mu/(1+tau*s)

% set parameters of disturbances: 
Namp=0.25;   % ratio of noise amplitude at low frequencies
tauH=tau/4;  % "known" input lowpass filter time constant

% set parameters of simulation
Tmax=20;       % years
Fsamp=52;       % number of samples per year

% set parameters of transfer function estimation:
nfft=2*Fsamp;   % number of samples per time-interval 
% (this determines both number of averages and minimum frequency resolved)
Fminmax=[0.0001 6];  % minimum and maximum frequency (cycles/year) to use
% in estimating mu and tau


% now generate simulation vectors...
dt=1/Fsamp;     % generate data at weekly sample rate
tvec=dt:dt:Tmax;nt=length(tvec);
% unknown disturbance:
n=Namp*randn(nt,1);
% convert continuous time model parameters to discrete time:
aH=1-dt/tauH;bH=dt/tauH;
if aH<0.5,disp('Time step too large');end
aT=1-dt/tau;bT=dt/tau;
% generate "known" disturbance:
H=filter(bH,[1 -aH],randn(nt,1));
% generate time sequence:    
T=filter(bT,[1 -aT],mu*H+n);

%% transfer function:
% set plotting controls for gest
% (for gest_simple, set figinfo=1)
figinfo.pflg=3;figinfo.Magfig=1;figinfo.Phfig=2;
figinfo.nrow=1;figinfo.ncol=1;
figinfo.fnsize=12;
figinfo.subplt=1;
figinfo.pltstr='HT';

% compare different amounts of averaging
for k=1:3,figure(k);hold off;end
figinfo.colour='k';
[Ght1,gerr_ht1,P,C,F] = gest(H,T,2*nfft,Fsamp,Fminmax,2,figinfo);
for k=1:3,figure(k);hold on;end
figinfo.colour='b';
[Ght2,gerr_ht2,P,C,F] = gest(H,T,nfft,Fsamp,Fminmax,2,figinfo);
for k=1:3,figure(k);xlim([0.2 15]);end
