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


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

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

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

nfft=input('enter nfft=')
nskip=input('enter number of points to skip at beginning')
ndat = n - nskip;
fname
gain1 = input('enter gain1=')
gain2 = input('enter gain2=')
gain1*gain2
scale = gain1*gain2*1000;

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

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

ch3 = ch1-ch2;
%%%ch3=ch3/scale;
ch3 = ch3-mean(ch3);

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

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

%ch6 = data(6,1+nskip:n)';
%ch6=ch6/scale;
%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)))
        
xlabel 'Frequency  (Hz)'
ylabel 'Amplitude  (V/sqrt(Hz))'
ht=title ([fname,',  N=',int2str(nfft)]);
set (ht,'Interpreter','none')
legend ('CH1','CH2','CH1-CH2')


[C12,F] = mscohere(ch1,ch2,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),C12(2:nf))
        
xlabel 'Frequency  (Hz)'
ylabel 'Coherence)'
ht=title ([fname,',  N=',int2str(nfft)]);
set (ht,'Interpreter','none')
legend ('Coh12')


%  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  (V/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))
      
xlabel 'Time (s)'
ylabel 'Amplitude  (V)'
ht=title ([fname, ', after DC-leveling']);
set (ht,'Interpreter','none')
legend ('CH1','CH2','CH1-CH2')
% figure
% ind = (0:ndat-1)*dt;
% plot (ind(1:ndat),ch2(1:ndat))
% xlabel 'Time (s)'
% ylabel 'Amplitude  (V)'
% 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  (V)'
% 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  (V)'
% 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  (V)'
% 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  (V)'
% ht=title ([fname, ', after DC-leveling']);
% set (ht,'Interpreter','none')
% legend ('CH6')

% 
% % histogram
% figure
% hist(ch1,101)
