% HFO_Detection_BENESCOClass
% Kaspar Schindler, 6.September 2015
%
% This script illustrates the use of the JSPECT algorithm with test
% signals.
%
% Reference: Niederhauser JJ, Esteller R, Echauz J, Vachtsevanos G, Litt B. Detection of
% seizure precursors from depth-EEG using a sign periodogram transform. IEEE Trans
% Biomed Eng. 2003 Apr;50(4):449-58. PubMed PMID: 12723056.
% Homework:
%%%%%%%%%%%%%
% 1. Go through the code and try to understand each step
% 2. Play with the amplitude of the fast oscillation...how large do
% you have to set it to detect it with normal FFT?
% 3. Try different test signals: for example add a "spike" at 14.5s
% (Ansatz: set some values for example to amplitude=4 and check
% if this affects JSPECT...
% 4. Replace the high frequency sine wave with a chirp signal (check
% Matlabs chirp command)
% 5. As an important general question: Why are monochromatic colormaps to be preferred? for answer see: http://www.cs.odu.edu/~mweigle/cs725s15/presentations/nam-presentation.pdf
%
% for questions, just email me: kaspar.schindler@gmail.com
clear all;clc % clear all variables and the command window
fs = 512; % define a sampling rate typical for EEG signals (though for the test signals used here, this does not really matter).
time = 1/fs:1/fs:30; % define a time vector of length 30s
signalA = randn(1,length(time)); % this creates a random TestSignal drawn form a normal distribution.
signalA(14*fs:15*fs) = 0.1*sin(2*pi*85*time(14*fs:15*fs)); % this gives us a small amplitude (amplitude=0.1) fast oscillation
% lets try to detect the fast oscillations by standard FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = fs; % this is the length of the moving window. By using 1s, the frequency resolution will then be 1Hz, which is convenient
dT = 1; % shift of the moving window in sampling points, dT = 1 is smalles possible advance of moving window, or put differently: maximal overlap for consecutive moving windows
nSteps = floor((length(signalA)-T)/dT); % how many times the moving window can be moved along TestSignal C
FFTcoll = zeros(101,nSteps); % in this variable the amplitude of the FFT will be saved for each time step and for 0-100Hz.
for iStep = 1:nSteps
dataInMovingWindow = signalA((iStep-1)*dT+1:(iStep-1)*dT+T);
ampFFTmovingWindow = abs(fft(dataInMovingWindow));
FFTcoll(:,iStep) = ampFFTmovingWindow(1:101);
end % end for iStep = 1:nSteps
% now we use the JSPECT algorithm: CAVE: carefully note the
% differences
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
JSPECTcoll = zeros(101,nSteps); % in this variable the amplitude of the FFT will be saved for each time step and for 0-100Hz.
for iStep = 1:nSteps
dataInMovingWindow = sign(diff(signalA((iStep-1)*dT+1:(iStep-1)*dT+T))); % CAVE: have you seen the sign(diff()) operator
ampFFTmovingWindow = abs(fft(dataInMovingWindow));
JSPECTcoll(:,iStep) = ampFFTmovingWindow(1:101);
end % end for iStep = 1:nSteps
% now we plot the results to compare
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(100)
set(gcf,'color','w') % this sets figure background to white, which I prefer because of high contrast
subplot(311) % first plot the TestSignal
plot(time(1:size(JSPECTcoll,2)),signalA(1:size(JSPECTcoll,2)),'k');
axis tight % automatically adjusts the axis
title('test signal')
set(gca,'TickDir','out','box','off','XTickLabel','')
subplot(312) % now plot the spectrogram using the powerful imagesc command (cave: you need the image processing toolbox)
imagesc(FFTcoll,[0 100]) % note that the automatic scaling is overrided by defining the range. Otherwise the comparison with JSPECT would not be fair!
set(gca,'XTick',5*fs:5*fs:25*fs,'XTickLabel',{''},'TickDir','out','box','off')
axis xy % this will define the origin (0,0) in the left lower corner
title('FFT') % give the plot a title
ylabel('frequency [Hz]') % label the axis (as should be done ALWAYS), if you are minimalistic you may restrict yourself to the units (i.e. if it is Hz then it has to be frequency...)
subplot(313)
imagesc(JSPECTcoll,[0 100])
axis xy
title('JSPECT')
xlabel('time [s]')
ylabel('frequency [Hz]')
set(gca,'XTick',5*fs:5*fs:25*fs,'XTickLabel',{'5','10','15','20','25'},'TickDir','out','box','off')
colormap(bone) % using monochromatic color map (see question above)