function vfitsinus = gppolyf(data,freq,ipic) 
% GPTRIGF, General Purpose TRIGonometric Fit
%     This fits a sinusiod (with fit phase) to a set of data,
%        and makes plots for assessing goodness of fit
%
% Inputs:
%     data  -  n x 2 matrix of input data, 
%                  first column is x's, i.e. design points
%                      (Caution: don's use integers here with integer freq)
%                  second column is y's, i.e. data values
%                             OR
%              n x 1 vector of tie series data
%                  when use x = 0, 1/n, ..., (n-1)/n
%     freq   -  frequency of sinusoid to fit, on "0-1 scale"
%                  i.e. will fit model of form:
%                      m + A * cos(2*pi*freq*x) + B * sin(2*pi*freq*x)
%                  Note: integer frequencies are the Fourier Frequencies.
%                  To control period, use  freq = 1 / period
%     ipic  -  picture control (optional)
%                  0 - don't make plots
%                  1 - (default) make plots in a new figure window
%                  2 - make plots in current figure window
%
% Output:
%     vfitsinus    -  n x 1 column vector of sinusoidal fit to y's,
%                        evaluated at x's

%    Copyright (c) J. S. Marron 1998



%  set up picture parameter
%
if nargin == 2 ;    %  only 2 arguments input, use default iipic
  iipic = 1 ;
else ;              %  xgrid was specified, use that
  iipic = ipic ;
end ;


%  organize x's and y's for least square fit
%
n = size(data,1) ;
if size(data,2) == 2 ;     %  then already have x's
  xgrid = data(:,1) ;
  vy = data(:,2) ;
else ;       %  then need to create x's
  xgrid = ((0:(n-1)) / n)' ;
  vy = data(:,1) ;
end ;


%  Do conventional least square fitting
%
mx = [ones(n,1), cos(2 * pi * freq * xgrid), sin(2 * pi * freq * xgrid)] ;
vcoeff = (mx' * mx) \ (mx' * vy) ;
vfitsinus = mx * vcoeff ;




%  get sum of squares stuff
%
totss = sum((vy - mean(vy)).^2) ;
pss = sum((vfitsinus - mean(vfitsinus)).^2) ;
percpss = 100 * pss / totss ;
percrss = 100 - percpss ;




if iipic >= 1 ;     %  Then show graphics with result

  if iipic == 1 ;    %  then make a new figure window for graphics
    figure ;
  end ;

  subplot(2,1,1) ;
    plot(xgrid,[vy,vfitsinus]) ;
      titstr = ['L. S. Fit of Sin., Freq. ' num2str(freq)] ;
      titstr = [titstr, ', % of SS = ' num2str(percpss)] ;
      title(titstr) ;
      vax = axis ;


  subplot(2,1,2) ;
    plot(xgrid,(vy-vfitsinus)) ;
      title(['Residuals, % of SS = ' num2str(percrss)]) ;
      hold on ;
        plot([vax(1);vax(2)],[0;0],'--') ;
      hold off ;
      vaxhr = (vax(4) - vax(3)) / 2 ;
          %  half range of vertical axis
      axis([vax(1) vax(2) -vaxhr vaxhr]) ;
          %  make range of this axis same as for raw data,
          %  for good visual understanding


end ;

