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


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

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

%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 
%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'
            ch_gain1(i)=152;
        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 (firmware_version >= 4.2)
    scale = scale*0.4
end
scale

fs=header.sampling_frequency
dt = 1/fs;
[nch n]=size(data)

nskip=input('enter number of points to skip at beginning');
ndat=input('enter sample number of last point (=0 take all data points)= ');
nfft=input('enter nfft (default 10 000) =')

if isempty(nfft)
    nfft = 10000
end


if isempty(nskip)
    nskip = 0;
end
if isempty(ndat)
    ndat = 0;
end
if ndat > 0 
    n=ndat;
end

ndat = n - nskip;
fname

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

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

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

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

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

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

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)),F(2:nf),sqrt(P66(2:nf)))
xlabel 'Frequency  (Hz)'
ylabel 'Amplitude  (V/sqrt(Hz))'
ht=title ([fname,',  N=',int2str(nfft)]);
set (ht,'Interpreter','none')
legend ('CH1','CH2','CH3','CH4','CH5','CH6')


[C16,F] = mscohere(ch1,ch6,hanning(nfft),noverlap,nfft,fs);
[C26,F] = mscohere(ch2,ch6,hanning(nfft),noverlap,nfft,fs);
[C36,F] = mscohere(ch3,ch6,hanning(nfft),noverlap,nfft,fs);
[C46,F] = mscohere(ch4,ch6,hanning(nfft),noverlap,nfft,fs);
[C56,F] = mscohere(ch5,ch6,hanning(nfft),noverlap,nfft,fs);
figure
semilogx (F(2:nf),C16(2:nf),F(2:nf),C26(2:nf),F(2:nf),C36(2:nf), ...
        F(2:nf),C46(2:nf),F(2:nf),C56(2:nf))
xlabel 'Frequency  (Hz)'
ylabel 'Coherence'
ht=title ([fname,',  N=',int2str(nfft)]);
set (ht,'Interpreter','none')
legend ('Coh16','Coh26','Coh36','Coh46','Coh56')


%  figure
%  %semilogy (F(2:nf),sqrt(P11(2:nf)))
%  loglog (F(2:nf),sqrt(P44(2:nf)),F(2:nf),sqrt(P55(2:nf)),F(2:nf),sqrt(P66(2:nf)))
%          
%  xlabel 'Frequency  (Hz)'
%  ylabel 'Amplitude  (mV/sqrt(Hz))'
%  ht=title ([fname,',  N=',int2str(nfft)]);
%  set (ht,'Interpreter','none')
%  legend ('CH4','CH5','CH6')




% time series
figure
ind = (0:ndat-1)*dt;
plot (ind(1:ndat),ch1(1:ndat),ind(1:ndat),ch2(1:ndat),ind(1:ndat),ch3(1:ndat), ...
      ind(1:ndat),ch4(1:ndat),ind(1:ndat),ch5(1:ndat),ind(1:ndat),ch6(1:ndat))
xlabel 'Time (s)'
ylabel 'Amplitude  (V)'
ht=title ([fname, ', after DC-leveling']);
set (ht,'Interpreter','none')
legend ('CH1','CH2','CH3','CH4','CH5','CH6')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch2(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (mV)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH2')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch3(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (mV)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH3')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch4(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (mV)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH4')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch5(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (mV)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH5')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch6(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (mV)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH6')


% histogram
%figure
%hist(ch1,101)
