function [x1]=FitandRemove60Hz(x,t)
% Remove 60 Hz from waveform x versus t
% From Paul Krehbiel -- adapted by R.Sonnenfeld
% Sonnenfeld tested it on periods from 0.2 sec to
% 20 sec and it seems to work nicely.  Probably
% best to use it for periods of 1-5 seconds. 
% DOES NOT WORK for files with a lot of lightning in the
% time period -- for this reason -- I have created Fit60Hz
% which merely returns fit coefficients for a subset of data
% (Whatever you give it).
% Then there is another file -- Remove 60Hz -- which applied the
% coefficients to the data.
%
% For really long times the 60 Hz phase drifts, making
% the data not fittable with this simple method.

% Solve overdetermined matrix equation
%       |cos(wt_1)  sin(wt_1)||a0|    | x_1 |
%       |cos(wt_2)  sin(wt_2)||a1|    | x_2 |
%       |   .           .    |      = |  .  |
%       |   .           .    |        |  .  |
%       |cos(wt_n)  sin(wt_n)|        | x_n |
%
% or               A     *     a    =    x
%
%         | a0 |
% for a = |    | or in matlab notation a = [a0 a1]' 
%         | a1 | 
% where a0,a1 are coefficients of least square fit
%
%       x_fit(t) = a_0*cos(wt) + a_1*sin(wt)
%
% of 60 Hz sine wave to actual data values x(t).
%
% In matlab, a = A^{-1}x == A\x  gives the least square solution for a.
% (\ = left matrix division - see 'help mldivide')
%
% t and x are both expected to be column vectors
% so size(t) or size(x) would return (N, 1)

% w = 2*pi*60;  %If you change from 60, you can remove 50 Hz -- or anything else sinusoidal!
% A = [cos(w*t) sin(w*t) ];
% %%A = [cos(w*t) sin(w*t) cos(3*w*t) sin(3*w*t)]; %This will generate
% % [a0 a1 a2 a3] and take out third harmonics too  -- but it's not needed
% disp('Fitting 60Hz to data')
% a = A\x;           % This is all the action.  Now a0 and a1 are determined.
% x_fit = A*a;       % Use a0 and a1 to generate x_fit per equation (1)
% disp('Removing 60Hz from data')
% x1 = x - x_fit;
%% Alternate Filter  Lapierre 02/29/2012
Fs=50000; % sampling frequency
f0=60; % frequency to notch
df=1; % notch width
nyq=Fs/2; % Nyquist frequency

[Bb, Ab]=butter(1, [f0-df/2 f0+df/2]/nyq, 'stop');
x1=filtfilt(Bb,Ab,x);
% plot(t,x)