close all
clear all
ch1 =[];
ch2 =[];
ch3 =[];
ch4 =[];
ch5 =[];
ch6 =[];


[fname,pathname]=uigetfile('*.*','Select data file');

[header,data] = readKMS([pathname fname]);


header.GPS_pos_info
header.GPS_time_info

ch_sort=input('enter method 1= HxHyHzEyEy (def), 2=ExEyHxHyHz (India)');
if isempty(ch_sort);
    ch_sort = 1;   % default
end
if (ch_sort == 1)
    ch_pos = [1 2 3 4 5];
end
if (ch_sort == 2)
    ch_pos = [3 4 5 1 2];
end
ch_pos

fs=header.sampling_frequency
dt = 1/fs;
[nch n]=size(data)
nfft=input('enter nfft=')
if isempty(nfft)
    nfft = 0
end
nskip=input('enter number of points to skip at beginning');
if isempty(nskip)
    nskip =0;
end
ndat=input('enter sample number of last point (=0 take all samples) = ');
if isempty(ndat)
    ndat=0;
end
if ndat > 0 
    n=ndat
end


%gain1 = input('enter gain1=')
%gain2 = input('enter gain2=')
gain1 = 1;
gain2 = 1;
gain1*gain2
scale = gain1*gain2;

ndat = n-nskip
ch1 = data(ch_pos(1),1+nskip:n)';
ch1=ch1/scale;
%ch1 = detrend(ch1);
ch1 = ch1-mean(ch1);


ch2 = data(ch_pos(2),1+nskip:n)';
ch2=ch2/scale;
%ch2 = detrend(ch2);
ch2 = ch2-mean(ch2);

ch3 = data(ch_pos(3),1+nskip:n)';
ch3=ch3/scale;
%ch3 = detrend(ch3);
ch3 = ch3-mean(ch3);

ch4 = data(ch_pos(4),1+nskip:n)';
ch4=ch4/scale;
%ch4 = detrend(ch4);
ch4 = ch4-mean(ch4);

ch5 = data(ch_pos(5),1+nskip:n)';
ch5=ch5/scale;
%ch5 = detrend(ch5);
ch5 = ch5-mean(ch5);

%ch6 = data(6,1+nskip:n)';
%ch6=ch6/scale;
%ch6 = detrend(ch6);
%ch6 = ch6-mean(ch6);

if (nfft > 0)
    
    nf = nfft/2 +1
    noverlap=nfft/2

% rms for each time series
    rms1=sqrt(ch1(1:ndat)'*ch1(1:ndat)/ndat);
    rms2=sqrt(ch2(1:ndat)'*ch2(1:ndat)/ndat);
    rms3=sqrt(ch3(1:ndat)'*ch3(1:ndat)/ndat);
    rms4=sqrt(ch4(1:ndat)'*ch4(1:ndat)/ndat);
    rms5=sqrt(ch5(1:ndat)'*ch5(1:ndat)/ndat);
%rms6=sqrt(ch6(1:ndat)'*ch6(1:ndat)/ndat);

% power density spectrum
    [P11 F]=pwelch(ch1,hanning(nfft),noverlap,nfft,fs,'onesided');
    [P22 F]=pwelch(ch2,hanning(nfft),noverlap,nfft,fs,'onesided');
    [P33 F]=pwelch(ch3,hanning(nfft),noverlap,nfft,fs,'onesided');
    [P44 F]=pwelch(ch4,hanning(nfft),noverlap,nfft,fs,'onesided');
    [P55 F]=pwelch(ch5,hanning(nfft),noverlap,nfft,fs,'onesided');
%[P66 F]=pwelch(ch6,hanning(nfft),noverlap,nfft,fs,'onesided');

    pow_f1=sqrt(sum(P11)/nfft/dt);
    pow_f2=sqrt(sum(P22)/nfft/dt);
    pow_f3=sqrt(sum(P33)/nfft/dt);
    pow_f4=sqrt(sum(P44)/nfft/dt);
    pow_f5=sqrt(sum(P55)/nfft/dt);
%pow_f6=sqrt(sum(P66)/nfft/dt);

    pow_f1,rms1
    pow_f2,rms2
    pow_f3,rms3
    pow_f4,rms4
    pow_f5,rms5
%pow_f6,rms6

    figure
%semilogy (F(2:nf),sqrt(P11(2:nf)))
    loglog (F(2:nf),sqrt(P11(2:nf)),F(2:nf),sqrt(P22(2:nf)),F(2:nf),sqrt(P33(2:nf)), ...
        F(2:nf),sqrt(P44(2:nf)),F(2:nf),sqrt(P55(2:nf)))
    xlabel 'Frequency  (Hz)'
    ylabel 'Amplitude  (mV/sqrt(Hz))'
    ht=title ([fname,',  N=',int2str(nfft)]);
    set (ht,'Interpreter','none')
    legend ('CH1 Bx','CH2 By','CH3 Bz','CH4 Ex','CH5 Ey')


    [C51,F] = mscohere(ch5,ch1,hanning(nfft),noverlap,nfft,fs);
    [C42,F] = mscohere(ch4,ch2,hanning(nfft),noverlap,nfft,fs);
    [C52,F] = mscohere(ch5,ch2,hanning(nfft),noverlap,nfft,fs);
    [C41,F] = mscohere(ch4,ch1,hanning(nfft),noverlap,nfft,fs);

    figure
    semilogx (F(2:nf),C51(2:nf),F(2:nf),C42(2:nf), ...
          F(2:nf),C52(2:nf),F(2:nf),C41(2:nf))
    xlabel 'Frequency  (Hz)'
    ylabel 'Coherence'
    ht=title ([fname,',  N=',int2str(nfft)]);
    set (ht,'Interpreter','none')
    legend ('Coh Ey&Bx','Coh Ex&By','Coh Ey&By','Coh Ex&Bx')

    [C31,F] = mscohere(ch3,ch1,hanning(nfft),noverlap,nfft,fs);
    [C32,F] = mscohere(ch3,ch2,hanning(nfft),noverlap,nfft,fs);

    figure
    semilogx (F(2:nf),C31(2:nf),F(2:nf),C32(2:nf))
          
    xlabel 'Frequency  (Hz)'
    ylabel 'Coherence'
    ht=title ([fname,',  N=',int2str(nfft)]);
    set (ht,'Interpreter','none')
    legend ('Coh Bz&Bx','Coh Bz&By')
end
% time series
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),data(ch_pos(1),1+nskip:n)')
xlabel 'Time (s)'
ylabel 'Amplitude  (mV)'
ht=title ([fname, ', raw data']);
set (ht,'Interpreter','none')
legend ('CH1 Bx')
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),data(ch_pos(2),1+nskip:n)')
xlabel 'Time (s)'
ylabel 'Amplitude  (mV)'
ht=title ([fname, ', raw data']);
set (ht,'Interpreter','none')
legend ('CH2 By')
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),data(ch_pos(3),1+nskip:n)')
xlabel 'Time (s)'
ylabel 'Amplitude  (mV)'
ht=title ([fname, ', raw data']);
set (ht,'Interpreter','none')
legend ('CH3 Bz')
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),data(ch_pos(4),1+nskip:n)')
xlabel 'Time (s)'
ylabel 'Amplitude  (mV)'
ht=title ([fname, ', raw data']);
set (ht,'Interpreter','none')
legend ('CH4 Ex')
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),data(ch_pos(5),1+nskip:n)')
xlabel 'Time (s)'
ylabel 'Amplitude  (mV)'
ht=title ([fname, ', raw data']);
set (ht,'Interpreter','none')
legend ('CH5 Ey')

