%advobs
% ADVOBS.m estimates the frequency, amplitude and phase of a sine wave interference
% in OBS data records (or selected bursts) from the 2001 SWW deployment, by applying 
% FINDFREQ.  ADVOBS also produces some statistics on the amplitudes and phases, 
% so values can be chosen for cleaning the entire record. 
%
% J. Lacy, USGS, November 2001
%
% Pick some bursts:
nmin = 6;
nmax = 205;
brange = [nmin:nmax];
nb = length(brange);
sta = input('enter station name in single quotes ')
%fn = '/work2/sherwood/data/GH2001/d1/ND/ND1001.adr'
%fn = '/work2/sherwood/data/GH2001/d1/SD/SD1001.adr'
%fn = '/work2/sherwood/data/GH2001/d1/MD/MD1001.adr'
%fn = '/work2/sherwood/data/GH2001/d1/MIB/MIB1001.adr'
%fn = '/work2/sherwood/data/GH2001/d1/MS/MS1001.adr'
%
%!! Note: when running MS (dep 1 or 2), you need to go into findfreq.m and change
% nlow from 5 to 2!!!
%
%fn = '/work3/jlacy/GH2001/GH2001d2/ND/ADV_B45B/ND2002.adr';
fn = '/work3/jlacy/GH2001/GH2001d2/SD/ADV_B52B/SD2001.adr';
%fn = '/work3/jlacy/GH2001/GH2001d2/MD/ADV_B51B/MD2002.adr';
%fn = '/work3/jlacy/GH2001/GH2001d2/MIB/ADV_B59B/MIB2001.adr';
%fn = '/work3/jlacy/GH2001/GH2001d2/MS/ADV_B158/MS2001.adr';
%
fr1a = zeros(nb,1);
fr2a = zeros(nb,1);
amponea = zeros(nb,1);
amptwoa = zeros(nb,1);
phaonea = zeros(nb,1);
phatwoa = zeros(nb,1);
mn1 = zeros(nb,1);
mn2 = zeros(nb,1);
std1 = zeros(nb,1);
std2 = zeros(nb,1);
temp = zeros(nb,1);
%A = adr2mat(fn,[nmin,nmax]);
for i = nmin:nmax;
  ii = i - nmin+1;
A = adr2mat(fn,[i,i]);
[nr,nc] = size(A.extsensor{1});
del_t = 1/A.samprate(1);
 [obsgood1, frmax1, amp1, pha1] = findfreq(A.extsensor{1}(1:nr,1), del_t);
%  axis ([0 n 200 1600])
  title (['OBS 1, Burst ',num2str(i)])
%  pause;
 [obsgood2, frmax2, amp2, pha2] = findfreq(A.extsensor{1}(1:nr,2), del_t);
%   axis ([0 n 200 1600])
    title (['OBS 2, Burst ',num2str(i)])
%   pause;
   fr1a(ii) = frmax1/del_t;
   amponea(ii) = amp1;
   phaonea(ii) = pha1;
   fr2a(ii) = frmax2/del_t;
   amptwoa(ii) = amp2;
   phatwoa(ii) = pha2;
    mn1(ii) = mean(A.extsensor{1}(1:nr,1));
   mn2(ii) = mean(A.extsensor{1}(1:nr,2));
   std1(ii) = std(A.extsensor{1}(1:nr,1));
   std2(ii) =  std(A.extsensor{1}(1:nr,2));
   temp(ii) = mean(A.temperature{1}(:,1));
end
figure 
subplot (3,1,1)
bnum = [nmin:nmax];
plot (bnum,amponea,'bo',bnum,amptwoa,'g+')
title ('Amplitudes')
ylabel ('Counts')
subplot (3,1,2)
plot (bnum, fr1a,'bo',bnum,fr2a,'g+')
title ('Frequencies')
ylabel ('s^-1')
subplot (3,1,3)
plot (bnum, phaonea, 'bo',bnum, phatwoa,'g+')
title('Phase')
ylabel ('Radians')
xlabel ('Burst number')
orient tall

figure 
subplot (2,2,1)
plot (mn1, amponea,'o',mn2,amptwoa,'+')
xlabel ('Burst mean')
ylabel ('Amplitude')
subplot (2,2,2)
plot (std1, amponea, 'o',std2, amptwoa, '+')
xlabel ('Burst std')
ylabel ('Amplitude')
subplot (2,2,3)
plot (mn1, fr1a,'o',mn2,fr2a,'+')
ylabel ('Frequency')
xlabel ('Burst mean')
subplot (2,2,4)
plot (temp, amponea, 'o',temp,amptwoa,'+')
xlabel('Temperature')
ylabel('Amplitude')

ione = find(mn1<1000);
itwo = find(mn2<1000);

eval([sta,'mnamp1 = mean(amponea(ione))'])
eval([sta,'mnamp2 = mean(amptwoa(ione))'])

eval([sta,'medamp1 = median(amponea(ione))'])
eval([sta,'medamp2 = median(amptwoa(ione))'])

eval([sta,'mnfr1a = mean(fr1a(ione))'])
eval([sta,'mnfr2a = mean(fr2a(ione))'])

eval([sta,'medfr1a = median(fr1a(ione))'])
eval([sta,'medfr2a = median(fr2a(ione))'])

[h,sig] = ttest2(amponea(ione), amptwoa(itwo), 0.05, 0)

eval([sta,'amp1 = amponea;']);
eval([sta,'amp2 = amptwoa;']);
eval([sta,'fr1 = fr1a;']);
eval([sta, 'fr2 = fr2a;']);
eval([sta, 'pha1 = phaonea;']);
eval([sta, 'pha2 = phatwoa;']);

eval(['save ', sta, 'obsstats ', sta, '* mn1 mn2 ione itwo temp bnum'])
