% ANALYSE SPEECH ABR
% SLB 16/5/2018
% This is rough Matlab code used to read in data and carry out time and
% frequency bootstraps
%
% ...it is not very well documented!!.. 
%
% but may be helpful to see how to read in
% data and to see how the basic analysis works
clear;close all
load P16da1.txt % THIS NEEDS TO HAVE NO HEADER, COLUMN 1 = DATA
y=P16da1(:,2); % channel here
% P16 has clear speech ABR
x1=resample(y,1,2); % downsample to 5000Hz - shouldn't make much difference


%fs=10000;
fs=5000;
t=0.090; % epoch length in signal - THIS WILL ADD THE POLARITIES TOGETHER
% TWO EPOCHS OF 90 ms, Stimulus at 10 ms and alternating at 100 ms

epochlen=t*fs;
L=epochlen

%%B=butter(2,100/2500,'high'); % introduce 200 Hz filter into signal
%B=butter(5,[140/2500 160/2500],'stop');
%B=butter(10,[148/2500 152/2500],'stop'); % introduce 100 Hz filter into signal
%freqz(B,1)
%figure;psd(x1,4096,fs)
[B]=notch_filter(100,0.1,fs);
x=filtfilt(B,1,x1);
[B]=notch_filter(150,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(200,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(250,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(300,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(350,0.1,fs);
x=filtfilt(B,1,x1);
[B]=notch_filter(400,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(450,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(500,0.1,fs);
x=filtfilt(B,1,x);
[B]=notch_filter(550,0.1,fs);
x=filtfilt(B,1,x);
%x=filtfilt(B,1,x1); %% FILTER NOT WORKING
%figure;plot(x)


N=length(x)/(t*fs)

%figure;plot(x);title('raw data')
%figure;plot(abs(fft(x(1:10000))));title('fft of 1s of data')
%figure;psd(x,4096,fs)
array=reshape(x,epochlen,N);
scale=linspace(0,t,epochlen);
average=mean(array');
average=average/10000; % scale for amplifier gain
figure;plot(scale,average);title('ABR pre artefact rejection')
xlabel('time s')
ylabel('uV')

%ADD FSP ON RAW WAVEFORM
SIG=var(average(70:90)); % TIME WINDOW FROM 4 TO 8 ms post stimulus - for ABR
NOISE=var(array(80,:));          % Single point at 6 ms post stimulus
FSP=SIG*N/NOISE

SIG1=var(average(150:170)); % TIME WINDOW FROM 4 TO 8 ms post stimulus - for ABR
NOISE1=var(array(160,:));          % Single point at 6 ms post stimulus
FSP1=SIG1*N/NOISE1

SIG2=var(average(170:190)); % TIME WINDOW FROM 4 TO 8 ms post stimulus - for ABR
NOISE2=var(array(180,:));          % Single point at 6 ms post stimulus
FSP2=SIG2*N/NOISE2

% compare sub averages
array1=array(:,1:3000)/10000;
array2=array(:,3001:6000)/10000;
figure;plot(scale,mean(array1'))
hold on
plot(scale,mean(array2'),'c');title('two subaverages pre artefact rejection')
xlabel('time s')
ylabel('uV')
%axis([0 0.1 -0.01 0.01]) % change ABR scale here - for alternating

prestimulusnoise=std(average(1:50)) % changed 24/7 to be 10ms

%%%%%

% arterfact rejection
count=1;
%arraynew=0;
for i=1:N
    sweep=array(:,i);
    %temp=sweep(100:end); % cut off stimulus artefact
    temp=sweep;
    ptop=max(temp)-min(temp);
    if ptop<0.8 % 0.1? %0.04 % change this back?
        arraynew(:,count)=sweep;
        count=count+1;
    end
end
count
%figure;plot(scale,mean(arraynew'))
%axis([0 0.06 -0.008 0.008]) % change ABR scale here
average=mean(arraynew');
SIG1=var(average(140:180)); % for speech
%SIG=var(mean(array((40:80),:))); % TIME WINDOW FROM 4 TO 8 ms
NOISE1=var(array(160,:));          % Single point at 20 ms
FS1P=SIG1*N/NOISE1


%% time-frequency analsysis
% try windows of 24 samples. Hanning is default. 5 sample overalap
% may need to define window, number of frequency bins etc.
% would help to have same bins as Kraus et al.
%? default = hanning. Use blackman?
% using window 256 with 250 overlap is similar to adobe audition with 128
% points
window=blackman(256); % TRY USING BLACKMAN WINDOW
%lwindow=48;
noverlap=250;
[Y,F,T,P]=spectrogram(average,window,noverlap,[],fs); % P is the output array
%[Y,F,T,P]=spectrogram(average,lwindow,noverlap,[],fs); % P is the output array
%figure;surf(T,F,10*log10(abs(P)),'EdgeColor','none');view(0,90);
Timefreq=10*log10(P); % save specgram output (in dB?)

V=var(average(50:150))

%count=N % no rejection
for n=1:499
    temp=array';
    for i=1:count
        rotate=fix(rand*(L-1))+1; % avoid zero by adding 1
        temp(i,:)=[temp(i,rotate+1:L) temp(i,1:rotate)]; %% 
    end
    bootstrap(n,:)=mean(temp);
    [Y,F,T,P]=spectrogram(bootstrap(n,:),window,noverlap,[],fs); % P is the output array
    bootstrapTF(n,:,:)=10*log10(P); % put specgram of new average in array
    Bootvar(n)=var(bootstrap(n,50:150));
end

Bootvar=sort(Bootvar);
Bootvar(475)

for n=1:L
    temp2=bootstrap(:,n);
    temp1=sort(temp2);
    lower5(n)=temp1(25); %lowest 1%
    upper5(n)=temp1(475); % upper 1%    
    lower1(n)=temp1(5); %lowest 1%
    upper1(n)=temp1(495); % upper 1%
    lowerp2(n)=temp1(1); %lowest 1%
    upperp2(n)=temp1(499); % upper 1%    
end
figure;
plot(scale,average);hold on
plot(scale,lower5,'c')
plot(scale,upper5,'g')
plot(scale,lower1,'c:')
plot(scale,upper1,'g:')
plot(scale,lowerp2,'r:')
plot(scale,upperp2,'r:')

title('mean response with upper and lower 5th, 1st percentile amplitudes from bootstrap and upper and lower .2%')

%bootstrap specgram
for i=1:size(Timefreq,1);
    for n2=1:size(Timefreq,2);
        temp2=bootstrapTF(:,i,n2); % need to sort each time freq point, find
        temp3=sort(temp2);
        if Timefreq(i,n2)>temp3(495) % 1%
            Masker1pc(i,n2)=1; % put 1 if exceeded specgram at 1%
        else
            Masker1pc(i,n2)=0;
        end
        if Timefreq(i,n2)>temp3(475) % 1%
            Masker5pc(i,n2)=1; % put 1 if exceeded specgram at 5%
        else
            Masker5pc(i,n2)=0;
        end
    end
end
figure;surf(T,F,Masker1pc);view(0,90);title('spectral bootstap at 1%')
xlabel('time ms')
ylabel('freq(Hz)')
figure;surf(T,F,Masker5pc);view(0,90);title('spectral bootstap at 5%')
xlabel('time ms')
ylabel('freq(Hz)')