function [obsgood, amplsf, pha,error] = obscllsf(data, del_t,fr, amp)
% OBSCLLSF  uses least-squares fit to find the amplitude and phase of  
% sine wave interference (due to multiple OBSs) for single bursts of 
% OBS data, using a frequency and amplitude input by the user (for instance, 
% determined from ADVOBS). Then the burst is cleaned using the 
% pre-determined frequency and the burst-specific amplitude and phase.
%
%[obsgood, amplsf, pha, error] = obscl(data, del_t, fr, amp)
%  input data = raw data
%        del_t = sampling interval, seconds
%        fr = frequency of signal to be removed (determined by
%        findfreq.m)
%        amp = amplitude of signal to be removed 
%            (determined by findfreq.m), in counts
%
%  output  obsgood = clean data    
%        amplsf = amplitude of signal for this burst determined lsf.
%           least-squares fit, in counts.
%	 pha = phase lag of removed signal for each burst, in radians
%    error = mse error of data vs. least-squares fit (I think)
%
% Jessie Lacy, USGS
% November 1, 2001
% with help from tide_fit.m by Rich Signell

n = length(data);
mn = mean(data);
datadem = data-mn;
t = [del_t:del_t:n*del_t]';
fr_rad = 2*pi*fr;

 D=zeros(n,3); 
       D(:,1)=ones(n,1);
           D(:,2)=cos(fr_rad*t);
           D(:,3)=sin(fr_rad*t);
%-------solve [D coeff = data] -----------------
           coef=D\data;
%-------generate solution components-----
       mean=coef(1);
       bad=coef(2)*cos(fr_rad*t)+coef(3)*sin(fr_rad*t);
       error=std(D*coef-data)/sqrt(n-3);
		
amplsf = sqrt(coef(2)^2+coef(3)^2);
pha = atan2(coef(3),coef(2));
badtest = amplsf*cos(2*pi*fr*t -pha);
bad_amp = amp*cos(2*pi*fr*t - pha);
% plot
%
figure(1) 
clf
subplot (3,1,1)
plot (t, bad, t, badtest,'+') 
%
%clean data by subtracting out interfering signal using predetermined amplitude
%
obsgood = data-bad_amp;
%
%plot results
%
offset = max(bad - datadem);
imin = find((bad - datadem)==offset);
bad2 = bad - (bad(imin(1))-datadem(imin(1)))+ mn;
%bad2 = bad+mn;
%figure 
subplot (3,1,2)
plot (t, data, t, bad2,'r')
title ('Raw data (blue) and identified interference(red)')
axis ([0 n*del_t 0 2000])
subplot (3,1,3)
plot (t, obsgood)
title ('Clean data')
axis ([0 n*del_t 0 2000])



