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


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

text = {'enter method: 1= HxHyHzEyEy (def)'; ...
        '              2= ExEyHxHyHz (India)'; ...
        '              3= fluxgate (7 8 9 4 5)'}
%ch_sort=input('enter method 1= HxHyHzEyEy (def), 2=ExEyHxHyHz (India)');
ch_sort=input('method=');
if isempty(ch_sort);
    ch_sort = 1;   % default
end
if (ch_sort == 1)
    ch_pos = [1 2 3 4 5];
elseif (ch_sort == 2)
    ch_pos = [3 4 5 1 2];
elseif (ch_sort == 3)
    ch_pos = [7 8 9 4 5];
end
ch_pos


if (ch_sort == 3)
    [header,data] = readKMSflux([pathname fname]);
else
    [header,data] = readKMS([pathname fname]);
end
%header.GPS_pos_info
%header.GPS_time_info

%get firmware version
fw = dec2bin(header.firmware_version,8);
fw = (fw*1) - 48;
firmware_version = fw(1)*2^3+fw(2)*2^2+fw(3)*2+fw(4);
fw_dec = fw(5)*2^3+fw(6)*2^2+fw(7)*2+fw(8);
if (fw_dec > 9)
    firmware_version = firmware_version + 0.01*fw_dec;
else
    firmware_version = firmware_version + 0.1*fw_dec;
end
firmware_version 
% since firmware version has been resetted 4.2 --> 2.8
%fpga_version is now used to distinguish beteen old and new board
fpga_version = header.fpga_version
%get gain settings
s_char=dec2bin(header.ch_setting,8);
s_bit=s_char(8:-1:1,8:-1:1);
s_rev=reshape(s_bit',1,64);
%s_rev
for i=1:6
    ii=1+(i-1)*2;
    switch (s_rev(ii:ii+1));
        case '00'
            ch_gain1(i)=1;
        case '10'  %reversed bit order
            ch_gain1(i)=1;
        case '01' %reversed bit order
            ch_gain1(i)=8;
        case '11'
            if (fpga_version >= 64)
                ch_gain1(i)=152;
            else
                ch_gain1(i)=40;
            end
        otherwise
            error('This is impossible')
    end
    ii=13+(i-1)*2;
    switch (s_rev(ii:ii+1));
        case '00'
            ch_gain2(i)=1;
        case '10' %reversed bit order
            ch_gain2(i)=4;
        case '01' %reversed bit order
            ch_gain2(i)=16;
        case '11'
            ch_gain2(i)=64;
        otherwise
            error('This is impossible')
    end
    ch_gain(i)=ch_gain1(i)*ch_gain2(i);
    
end
ch_gain1
ch_gain2
ch_gain
scale = ch_gain.*1000;
if (fpga_version >= 64)
    scale = scale*0.4
end
scale

fs=header.sampling_frequency
dt = 1/fs;
[nch n]=size(data)
nfft=input('enter nfft (=0 skip spectra & coherence) = ')
if isempty(nfft)
    nfft = 0
end
nskip=input('enter number of points to skip at beginning (=0 default) = ');
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

ndat = n-nskip

if (nfft > ndat)
    nfft = ndat
end

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


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

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

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

ch5 = data(ch_pos(5),1+nskip:n)';
ch5=ch5/scale(5);
%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
    rms(1)=sqrt(ch1(1:ndat)'*ch1(1:ndat)/ndat);
    rms(2)=sqrt(ch2(1:ndat)'*ch2(1:ndat)/ndat);
    rms(3)=sqrt(ch3(1:ndat)'*ch3(1:ndat)/ndat);
    rms(4)=sqrt(ch4(1:ndat)'*ch4(1:ndat)/ndat);
    rms(5)=sqrt(ch5(1:ndat)'*ch5(1:ndat)/ndat);
%rms(6)=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_f(1)=sqrt(sum(P11)/nfft/dt);
    pow_f(2)=sqrt(sum(P22)/nfft/dt);
    pow_f(3)=sqrt(sum(P33)/nfft/dt);
    pow_f(4)=sqrt(sum(P44)/nfft/dt);
    pow_f(5)=sqrt(sum(P55)/nfft/dt);
%pow_f(6)=sqrt(sum(P66)/nfft/dt);
    pow_f
    rms

    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  (V/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')

