function [ghat,gerr,P,C,F] = gest(x,y,nfft,Fsamp,Fmax,flag,figinfo,fitval)
%function [ghat,gerr,P,C,F] = gest(x,y,nfft,Fsamp,Fminmax,flag,figinfo,fitval)
%
% estimate DC value of a transfer function from x to y
% if x is complex, then interpret the imaginary part as indicating the data
% quality (from TAO array); same for y.
%
% nfft is number of points to use in fft (more means more resolution, less
%       averaging)
% Fsamp is sample frequency
% Fmax is maximum real frequency to use in estimate (e.g. 6/year)
%      if length(Fmax)=2, then first element is minimum freq
% flag=1 then multiply xfer by s before averaging (i.e. fit an integrator)
% flag=2 then fit mu/(s+epsilon) (and ghat has two entries for mu, epsilon)
% %%% now part of figinfo!! pltstr=0 then no plot is generated, otherwise pltstr is the string label
% optional fitval over-rides estimates for fit to use on plot
% figinfo provides additional arguments to use to control plotting
% of multiple transfer functions... figure # for magnitude, phase, # of
% subplot columns & rows, subplot number, and font size
%
% Outputs are estimate and error, 
% plus xfer fn P and coherence C and freq F
%
% D. MacMynowski 9/25/08

% unfortunately, the following arguments aren't passed in:
newfig=0;
%fnsize=figinfo(4);  % for subplot; otherwise 18 is better
mrksize=8; % for subplot, otherwise 12 is better
if nargin<7,figinfo.pflg=0;end
if ~isstruct(figinfo), % old calling
    pltstr=figinfo;
    clear figinfo
    if pltstr(1)==0,figinfo.pflg=0;else,figinfo.pflg=1;end
    % NOTE missing new functionality?
    if pltstr(1)=='x',newfig=0;pltstr=pltstr(2:end);end
    figinfo.Magfig=1;
    figinfo.Phfig=2;
    figinfo.ncol=1;figinfo.nrow=1;figinfo.subplt=1;
    figinfo.fnsize=12;
    figinfo.pltstr=pltstr;
    figinfo.colour='k';
end
if nargin<6,flag=0;end
if length(Fmax)==2,
    Fmin=Fmax(1);
    Fmax=Fmax(2);
else,
    Fmin=eps;
end

if isreal(x),
    % 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)
    
    win=hamming(nfft);
    [P,F]=tfestimate(x,y,win,[],nfft,Fsamp);
    [C,F]=mscohere(x,y,win,[],nfft,Fsamp);
else,
    [P,F,C,k] = tfnan(real(x),real(y),imag(x),imag(y),-nfft,Fsamp);
    F=F';
end

if flag==1,P=P.*(2*pi*F);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
sig2=err2n.*P2;        % variance of measurement

% added DGM 12/9/09 to fit gamma/(s+epsilon) to data.
if flag==2,
    % maybe phase info is terrible?
    if 0,
        P=abs(P);
        M=[-P(If).*sqrt(C(If)) sqrt(C(If))];
        V=P(If).*(2*pi*sqrt(-1)*F(If)).*sqrt(C(If));
        ghatc=pinv(M)*V;
        ghatb=sum(C(If).*abs(2*pi*F(If).*P(If)))/sum(C(If));  % was abs(P), DGM 6/2/09
        % try fitting magnitude squared directly:
        P2=abs(P).^2;
        M=[-P2(If).*C(If) C(If)];
        V=(P2(If).*(2*pi*F(If)).^2).*C(If);
        ghat=pinv(M)*V;
        ghat=sqrt(ghat);
    end
    % updated...
    X=[ones(length(If),1) -P2(If)];
    Y=(2*pi*F(If).*abs(P(If))).^2;
    epsest=1;  % has to start nonzero if zero freq point included
    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.
    gerr=fliplr(flipud(sqrt(gerr)));
else,
    ghat=sum(C(If).*abs(P(If)))/sum(C(If));  % was abs(P), DGM 6/2/09
    % changed back to abs... note that with time delay, then phase
    % varies with frequency. (8/4/09)
    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
%gerr=sqrt(1/sum(C(If))/k);

if nargin==8,ghat=fitval;end

if figinfo.pflg>0,
    if ~isfield(figinfo,'colour'),figinfo.colour='k';end
    pltstr=figinfo.pltstr;
    fnsize=figinfo.fnsize;
    if flag==1,P=P./(2*pi*F);end
    %if newfig,figure;else,figure(1);end
    figure(figinfo.Magfig);
    subplot(figinfo.nrow,figinfo.ncol,figinfo.subplt)
    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',figinfo.colour)
    set(gca,'XScale','log')
    set(gca,'YScale','log')
    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--')
        %loglog(F,abs(ghatb)*(2*pi*F).^(-1),'b--')
        %ghat=[ghat;ghatb];
    end
    if flag==0,xlim([.4 15]);else,xlim([0.2 15]);end
    xlim([0.2 15])
    if F(2)>0.4,xlim([0.4 15]);end
    if ~flag,ylim([abs(ghat)/10 abs(ghat)*10]);end
    %xlabel('Frequency (yr^{-1})')
    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 T_{' pltstr '}'])
    %legend('TF estimate','TF fit')
    grid on
    set(gca,'XMinorGrid','off')
    set(gca,'YMinorGrid','off')
    if flag,scl=2*pi*2;else,scl=1;end
    if flag<2,
        ht=text(F(2)*1.2,abs(ghat)*5/scl,['T_{' pltstr '} = ' num2str(abs(ghat),2)]);
    else,
        ht=text(F(2)*1.2,abs(ghat(2))*5/scl,['T_{' pltstr '}: \mu = ' num2str(abs(ghat(2)),2) ', \epsilon = ' num2str(abs(ghat(1)),2)]);
    end
    set(ht,'FontSize',fnsize);set(ht,'BackGroundColor','w')
    drawnow
    
    if figinfo.nrow>1&figinfo.ncol>1,
        zz=get(gca,'OuterPosition');
    %textboxloc=[zz(1) zz(2)+zz(4) .025 .025];
    textboxloc=[zz(1)-.01 zz(2)+zz(4) .025 .025];
    hi=annotation('textbox',textboxloc,'String',['(' char('a'+figinfo.subplt-1), ')']);
    % the above need to be moved...
    set(hi,'LineStyle','none');
    set(hi,'FontSize',fnsize);
    end
    if figinfo.pflg>=2,  % plot phase too
        figure(figinfo.Phfig);
        subplot(figinfo.nrow,figinfo.ncol,figinfo.subplt)
        set(gca,'FontSize',figinfo.fnsize)
        he=errorbar(F,180/pi*(phase(P)),180/pi*atan(sqrt(err2)));
        set(gca,'FontSize',fnsize)
        xlim([0.2 15])
        if F(2)>0.4,xlim([0.4 15]);end
        set(he,'Marker','.')
        set(he,'MarkerSize',mrksize)
        set(he,'Color',figinfo.colour)
        set(gca,'XScale','log')
        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(['Phase (deg), T_{' pltstr '}'])
        grid on
        set(gca,'XMinorGrid','off')
        set(gca,'YMinorGrid','off')
        hold on
        if flag==1,
            semilogx(F,180/pi*phase(ghat)*ones(size(F))-90,'r--')
        elseif flag==0,
            semilogx(F,180/pi*phase(ghat)*ones(size(F)),'r--')
        else,
            semilogx(F,180/pi*phase(ghat(2)*(ghat(1)+2*pi*sqrt(-1)*F).^(-1)),'r--')
        end
    end
    if figinfo.pflg==3,  % plot coherence too
        if ~isfield(figinfo,'Cfig'),figure(3);else,
            figure(figinfo.Cfig);
        end
        subplot(figinfo.nrow,figinfo.ncol,figinfo.subplt)
        set(gca,'FontSize',figinfo.fnsize)
        he=semilogx(F,C);
        set(gca,'FontSize',fnsize)
        xlim([0.2 15])
        if F(2)>0.4,xlim([0.4 15]);end
        set(he,'Marker','.')
        set(he,'MarkerSize',mrksize)
        set(he,'Color',figinfo.colour)
        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(['Coherence, T_{' pltstr '}'])
        grid on
        set(gca,'XMinorGrid','off')
        set(gca,'YMinorGrid','off')
    end
end

