function [P,C,F,ghat,gerr] = gest_simple(x,y,nfft,Fsamp,Fminmax,flag,pltflg)
%function [P,C,F,ghat,gerr] = gest_simple(x,y,nfft,Fsamp,Fminmax,flag,pltflg)
%
% estimate transfer function from x to y
% D. MacMynowski 10/28/10 (edited from gest.m)
%
% OUTPUT:
%       P, C, F: transfer function, coherence, at frequency vector F
%   Optional:
%       ghat, gerr: estimate for average value of transfer function, and error
%
% INPUT:
%       x, y are data series (estimates response from x to y
%       nfft is number of points to use in fft
%           (more means more resolution, less averaging)
%           Default should work out to ~7 averages
%       Fsamp is sample frequency (defaults to one)
%       Fminmax is a two-element vector with the minimum and maximum
%           frequency to use in estimating ghat.  Defaults to [0,inf]
%       flag determines the type of fit...
%           flag=0: fit a constant
%           flag=1: fit an integrator (and Fmin automatically changed to
%                   eps if it was zero)
%           flag=2 then fit mu/(s+epsilon) (and output ghat has two entries for mu, epsilon)
%       pltflg = 1 to plot results
%
%   Note overlap of 50% and hamming windows are assumed (these could be
%   modified below)

% Original version D. MacMynowski 9/25/08

if nargin<7,pltflg=1;end
if nargin<6,flag=0;end
if nargin<5,Fminmax=[0 inf];end
if nargin<4,Fsamp=1;end
if nargin<3,nfft=floor(length(x)/4);end
if length(Fsamp)==0,Fsamp=1;end
if length(Fminmax)==0,Fminmax=[0 inf];end
if length(Fminmax)==2,
    Fmin=Fminmax(1);
    Fmax=Fminmax(2);
else,
    Fmin=0;
    Fmax=Fminmax;
end

% Note here: fix length of time (frequency) rather than number of averages
% compute number of averages (for estimating error in g)
M=length(x);
L=nfft;
noverlap = fix(0.5.*L); % number of overlap
LminusOverlap = L-noverlap;
k = fix((M-noverlap)/LminusOverlap);    % or, k=1+(M-L)/(L-O) = (M-O)/(L-O)

% Next 3 lines are the actual content of the code:
win=hamming(nfft);
[P,F]=tfestimate(x,y,win,[],nfft,Fsamp);
[C,F]=mscohere(x,y,win,[],nfft,Fsamp);


if nargout>3,
    if flag==1,P=P.*(2*pi*F);Fmin=max(eps,Fmin);end
    
    If=find(F<=Fmax&F>=Fmin);
    
    P2=abs(P).^2;               % magnitude-squared; fit this rather than magnitude
    err2n=(1-C)./C/(2*k);       % variance of normalized error
    sig1=sqrt(err2n).*abs(P);   % standard deviation of magnitude error
    
    % added DGM 12/9/09 to fit gamma/(s+epsilon) to data.
    if flag==2,
        X=[ones(length(If),1) -P2(If)];
        Y=(2*pi*F(If).*abs(P(If))).^2;
        epsest=0;
        % iterate to converge on fit...
        for kii=1:5,
            sig=(epsest^2+(2*pi*F(If)).^2)*2.*sig1(If);  %  *abs(P)
            W=diag(sig.^(-2));
            gerr=inv(X'*W*X);
            ghat=gerr*X'*W*Y;
            ghat=sqrt(ghat);
            epsest=ghat(2);
        end
        ghat=flipud(ghat);  % for consistency with previous.
    else,
        ghat=sum(C(If).*abs(P(If)))/sum(C(If));
        % Only fit magnitude, not complex value, since
        % with time delay, then phase varies with frequency.
        X=ones(length(If),1);
        W=diag(sig1(If).^(-2));
        Y=abs(P(If));
        gerr=inv(X'*W*X);
        ghat=gerr*X'*W*Y;
        gerr=sqrt(gerr);  % report std
    end
end


if pltflg>0,
    
    if flag==1,P=P./(2*pi*F);end
    
    %set(gca,'FontSize',fnsize)
    err2=(1-C)./C/(2*k);   % variance of normalized error
    he=errorbar(F,abs(P),2*abs(P).*sqrt(err2));
    set(he,'Marker','.')
    %set(he,'MarkerSize',mrksize)
    set(he,'Color','k')
    set(gca,'XScale','log')
    set(gca,'YScale','log')
    if nargout>3,
    hold on
    if flag==1,
        loglog(F,abs(ghat)*(2*pi*F).^(-1),'r--')
    elseif flag==0,
        loglog(F,abs(ghat)*ones(size(F)),'r--')
    else,
        loglog(F,abs(ghat(2)*(ghat(1)+2*pi*sqrt(-1)*F).^(-1)),'r--')
    end
    end
    if 0,
        set(gca,'XTick',[1/24 1/12 1/6 1/3 1/2 1]*12)
        set(gca,'XTickLabel',int2str([24;12;6;3;2;1]))
        xlabel('Period (months)')
        ylabel(['Transfer function'])
    end
    grid on
    set(gca,'XMinorGrid','off')
    set(gca,'YMinorGrid','off')
    
 end

