close all
clear all
ch =[];



[fname,pathname]=uigetfile('*.bin','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)=40;
        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
% scale(1)=scale(4);   % nur fuer Tascin experiment
% scale(2)=scale(5);

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

% delay filter definition
nz=1;
sgn=1;
%eta=input('enter eta for delay filter (=0 no filter, 1=nonrecursive')
eta=0;
if eta > 0
    nz =input('enter order of delay + 1 (fs/f0+1)')
    sgn=input('enter sign: +1=delay1, -1=delay2, no DC')
end
a(1:nz) = 0;
a(1) = 1;
a(nz) = sgn*1;
b(1:nz) = 0;
b(1)=1;
b(nz) = sgn*eta;


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


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

ndat = n - nskip;
fname

% %estimate period in signal
% for ich=1:nch
%     [nper,per_time,per_count,per_length]=TH_perest1(data(ich,1+nskip:n))
% end

nfft=input('enter nfft=')
%f0 = input('enter frequency for notch filter f0 in Hz = ')

if isempty(nfft) == 0    % spectra
    
    if nfft==0 
        error ('break with nfft=0')
    end


    for ich=1:nch
        ch = data(ich,1+nskip:n)';
        ch=ch./scale(ich);
 
        %ch = detrend(ch);
        if eta > 0
            if eta < 1
                ch=filter(a,b,ch);
            else
                ch=filter(a,1,ch);
            end
        end
        ch = ch-mean(ch);


        nf = nfft/2 +1;
        noverlap=nfft/2;

% rms for each time series
        rms(ich)=sqrt(ch(1:ndat)'*ch(1:ndat)/ndat);

% power density spectrum
         [Pxx F]=pwelch(ch,hanning(nfft),noverlap,nfft,fs,'onesided');

        pow_f(ich)=sqrt(sum(Pxx)/nfft/dt);
        P11(:,ich) =Pxx;
    end
    rms
    pow_f
    
    figure  % plot spectra
    for ich=1:nch
        switch (ich)
            case 1
                loglog (F(2:nf),sqrt(P11(2:nf,1)))
                legend ('CH1')
            case 2
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)))
                legend ('CH1','CH2')
            case 3
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)))
                legend ('CH1','CH2','CH3')
            case 4
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)))
                legend ('CH1','CH2','CH3','CH4')
            case 5
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)),F(2:nf),sqrt(P11(2:nf,5)))
                legend ('CH1','CH2','CH3','CH4','CH5')
            case 6
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)),F(2:nf),sqrt(P11(2:nf,5)),F(2:nf),sqrt(P11(2:nf,6)))
                legend ('CH1','CH2','CH3','CH4','CH5','CH6')
            otherwise
                error ('exceeding the maximum number of 6 channels')
        end
        
        xlabel 'Frequency  (Hz)'
        ylabel 'Amplitude  (V/sqrt(Hz))'
        ht=title ([fname,',  N=',int2str(nfft)]);
        set (ht,'Interpreter','none')
    end
end  % end of spectra
% 
% [C12,F] = mscohere(ch1,ch2,hanning(nfft),noverlap,nfft,fs);
% [C13,F] = mscohere(ch1,ch3,hanning(nfft),noverlap,nfft,fs);
% [C23,F] = mscohere(ch2,ch3,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),F(2:nf),C13(2:nf),F(2:nf),C23(2:nf))
%         
% xlabel 'Frequency  (Hz)'
% ylabel 'Coherence'
% ht=title ([fname,',  N=',int2str(nfft)]);
% set (ht,'Interpreter','none')
% legend ('Coh12','Coh13','Coh23')
% 
% 
% %  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
str = input('enter flag for time series plot (a=all, s=single plot):','s');
str
if (str == 'a')
    figure    % plot all together

    ind = (0:ndat-1)*dt;
    for ich=1:nch
        switch (ich)
            case 1
                plot (ind(1:ndat),data(1,1+nskip:n))
                legend ('CH1')
            case 2
                plot (ind(1:ndat),data(1,1+nskip:n),ind(1:ndat),data(2,1+nskip:n))
                legend ('CH1','CH2')
            case 3
                plot (ind(1:ndat),data(1,1+nskip:n),ind(1:ndat),data(2,1+nskip:n),ind(1:ndat),data(3,1+nskip:n))
               legend ('CH1','CH2','CH3')
            case 4
                plot (ind(1:ndat),data(1,1+nskip:n),ind(1:ndat),data(2,1+nskip:n),ind(1:ndat),data(3,1+nskip:n), ...
                      ind(1:ndat),data(4,1+nskip:n))
                legend ('CH1','CH2','CH3','CH4')
            case 5
                plot (ind(1:ndat),data(1,1+nskip:n),ind(1:ndat),data(2,1+nskip:n),ind(1:ndat),data(3,1+nskip:n), ...
                      ind(1:ndat),data(4,1+nskip:n),ind(1:ndat),data(5,1+nskip:n))
                legend ('CH1','CH2','CH3','CH4','CH5')
            case 6
                plot (ind(1:ndat),data(1,1+nskip:n),ind(1:ndat),data(2,1+nskip:n),ind(1:ndat),data(3,1+nskip:n), ...
                      ind(1:ndat),data(4,1+nskip:n),ind(1:ndat),data(5,1+nskip:n),ind(1:ndat),data(6,1+nskip:n))
                legend ('CH1','CH2','CH3','CH4','CH5','CH6')
            otherwise
                error ('exceeding the maximum number of 6 channels')
        end
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', raw data']);
        set (ht,'Interpreter','none')
    end
end     % end of time series plot, all together
      
if (str == 's')
% plot each time series in a single plot
    for ich=1:nch
        figure
        ind = (0:ndat-1)*dt;
        plot (ind(1:ndat),data(ich,1+nskip:n)')
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', raw data']);
        set (ht,'Interpreter','none')
        legend (['CH' int2str(ich)])
    end
end
% STACKING over a window, alpha-timmed mean

%nwlen=input('enter number points for periodicy for stacking');
nwlen = 0;
if isempty(nwlen)
    nwlen = 0;
end
if nwlen > 0
    alpha = input('enter alpha-trimmed mean value (0=mean,0.5=median): ');
    if isempty(alpha)
        alpha = 0;
    end
    for ich=1:nch
        ch = data(ich,1+nskip:n);
        for k=1:nwlen
            aux = ch(k:nwlen:length(ch));
            nalpha = round(alpha*length(aux));
            istart=1+nalpha;
            istop = length(aux)-nalpha;
            if (istart > istop)
                isw = istart;
                istart = istop;
                istop = isw;
            end
            if (nalpha >= 1)
                aux = sort(aux);
            end
            x(ich,k)=mean(aux(istart:istop));
        end
        length(aux)
%        nalpha
%        istart
%        istop
%     k=(1:nwlen);
%     x(ich,k)=mean(ch(k:nwlen:length(ch)));
    
        figure  % plot stacked time series
        ind = (0:nwlen-1)*dt;
        plot (ind(1:nwlen),x(ich,1:nwlen))
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', stacked data']);
        set (ht,'Interpreter','none')
        legend (['CH' int2str(ich)])
    end
if (str == 'a')
    figure    % plot all together

    ind = (0:nwlen-1)*dt;
    for ich=1:nch
        switch (ich)
            case 1
                plot (ind(1:nwlen),x(1,1:nwlen))
                legend ('CH1')
            case 2
                plot (ind(1:nwlen),x(1,1:nwlen),ind(1:nwlen),x(2,1:nwlen))
                legend ('CH1','CH2')
            case 3
                plot (ind(1:nwlen),x(1,1:nwlen),ind(1:nwlen),x(2,1:nwlen),ind(1:nwlen),x(3,1:nwlen))
               legend ('CH1','CH2','CH3')
            case 4
                plot (ind(1:nwlen),x(1,1:nwlen),ind(1:ndat),x(2,1:nwlen),ind(1:nwlen),x(3,1:nwlen), ...
                      ind(1:nwlen),x(4,1:nwlen))
                legend ('CH1','CH2','CH3','CH4')
            case 5
                plot (ind(1:nwlen),x(1,1:nwlen),ind(1:nwlen),x(2,1:nwlen),ind(1:nwlen),x(3,1:nwlen), ...
                      ind(1:nwlen),x(4,1:nwlen),ind(1:nwlen),x(5,1:nwlen))
                legend ('CH1','CH2','CH3','CH4','CH5')
            case 6
                plot (ind(1:nwlen),x(1,1:nwlen),ind(1:nwlen),x(2,1:nwlen),ind(1:nwlen),x(3,1:nwlen), ...
                      ind(1:nwlen),x(4,1:nwlen),ind(1:nwlen),x(5,1:nwlen),ind(1:nwlen),x(6,1:nwlen))
                legend ('CH1','CH2','CH3','CH4','CH5','CH6')
            otherwise
                error ('exceeding the maximum number of 6 channels')
        end
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', stacked data']);
        set (ht,'Interpreter','none')
    end
end     % end of time series plot, all together
      
if (str == 's')
% plot each time series in a single plot
    for ich=1:nch
        figure
        ind = (0:nwlen-1)*dt;
        plot (ind(1:nwlen),x(ich,1:nwlen))
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', stacked data']);
        set (ht,'Interpreter','none')
        legend (['CH' int2str(ich)])
    end
end
    
    
end

% delay hanning window (3-point)
% ndel=nwleng/2;
% x(1:ndat-ndel)=[0.5*(x(1:ndat-ndel)-x(1+ndel:ndat))];
% x(ndat:-1:ndel+1)=0.5*(x(ndat:-1:ndel+1)-x(ndat-ndel:-1:1));
% x(1:1:n)=-x(n+1:+1:2*n);
% x(ndat-n+1:ndat)=-x(ndat-2*n+1:1:ndat-n);


%nwlen=input('enter # sample of periodicy for DELAY HANNING (3-point)');
nwlen = 0;
if isempty(nwlen)
    nwlen = 0;
end
if nwlen > 0 
    for ich=1:nch
        ch = data(ich,1+nskip:n);
        ndel=nwlen/2;
        ch(1:ndat-ndel)=[0.5*(ch(1:ndat-ndel)-ch(1+ndel:ndat))];
        ch(ndat:-1:ndel+1)=0.5*(ch(ndat:-1:ndel+1)-ch(ndat-ndel:-1:1));
        ch(1:1:ndel)=0; %-ch(ndel+1:+1:2*ndel);
        ch(ndat-ndel+1:ndat)=0;  %-ch(ndat-2*ndel+1:1:ndat-ndel);
        figure
        ind = (0:ndat-1)*dt;
        plot (ind(1:ndat),ch)
        xlabel 'Time (s)'
        ylabel 'Amplitude  (mV)'
        ht=title ([fname, ', filtered data']);
        set (ht,'Interpreter','none')
        legend (['CH' int2str(ich)])
%      spectra of filtered data
        [Pxx F]=pwelch(ch,hanning(nfft),noverlap,nfft,fs,'onesided');

        pow_f(ich)=sqrt(sum(Pxx)/nfft/dt);
        P11(:,ich) =Pxx;
        
        
    end
end
if ((nfft > 0) & (eta > 0))
    figure  % plot spectra
    for ich=1:nch
        switch (ich)
            case 1
                loglog (F(2:nf),sqrt(P11(2:nf,1)))
                legend ('CH1')
            case 2
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)))
                legend ('CH1','CH2')
            case 3
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)))
                legend ('CH1','CH2','CH3')
            case 4
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)))
                legend ('CH1','CH2','CH3','CH4')
            case 5
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)),F(2:nf),sqrt(P11(2:nf,5)))
                legend ('CH1','CH2','CH3','CH4','CH5')
            case 6
                loglog (F(2:nf),sqrt(P11(2:nf,1)),F(2:nf),sqrt(P11(2:nf,2)),F(2:nf),sqrt(P11(2:nf,3)), ...
                        F(2:nf),sqrt(P11(2:nf,4)),F(2:nf),sqrt(P11(2:nf,5)),F(2:nf),sqrt(P11(2:nf,6)))
                legend ('CH1','CH2','CH3','CH4','CH5','CH6')
            otherwise
                error ('exceeding the maximum number of 6 channels')
        end
        
        xlabel 'Frequency  (Hz)'
        ylabel 'Amplitude  (V/sqrt(Hz))'
        ht=title ([fname,',  N=',int2str(nfft), 'filtered']);
        set (ht,'Interpreter','none')
    end
end
% end of spectra
%
% % histogram
% figure
% hist(ch1,101)