From e45e8ccd236aa4f3c7abc377f730cd79b65c06b4 Mon Sep 17 00:00:00 2001 From: Jorge Date: Sun, 5 Mar 2017 10:45:03 +0100 Subject: [PATCH] VVOR 1.0 --- analize.m | 154 +++++++++++++++++++++++++++++++++++++++++++++++++++++ fourier.m | 9 ++++ isOctave.m | 12 +++++ read.m | 95 +++++++++++++++++++++++++++++++++ saccades.m | 88 ++++++++++++++++++++++++++++++ vvor.m | 15 ++++++ 6 files changed, 373 insertions(+) create mode 100644 analize.m create mode 100644 fourier.m create mode 100644 isOctave.m create mode 100644 read.m create mode 100644 saccades.m create mode 100644 vvor.m diff --git a/analize.m b/analize.m new file mode 100644 index 0000000..c99b7c1 --- /dev/null +++ b/analize.m @@ -0,0 +1,154 @@ +% Jorge Rey Martinez 2017 +% Inputs +% t = time array +% e = eye velocity array +% h = head velocity array +% s = Boolean, true if register is supresed VVOR false if VVOR +function analize(t,e,h,s) +if isOctave + scrsz = get(0,'ScreenSize'); +else + scrsz = get(groot,'ScreenSize'); +end +if s == 1 + figure1 = figure('Name','VORS PLOTS','NumberTitle','off','Position',[1 scrsz(4)/1.2 scrsz(3)/1.01 scrsz(4)/1.2]); +else + figure1 = figure('Name','VVOR PLOTS','NumberTitle','off','Position',[1 scrsz(4)/1.2 scrsz(3)/1.01 scrsz(4)/1.2]); +end +figure(figure1); +%Original plot +subplot(3,2,1) +plot(t,h,'b',t,e,'r','LineWidth',1.25) +title('Test Results - RAW data') +xlabel('Time in secs') +ylabel('Velocity in deg/sec') +ylim([-400 +400]) +legend ('Head velocity','Eye velocity') +%Desaccaded plot +subplot(3,2,2) +if s == 1 + desacE = medfilt1(e,35); +else + desacE = medfilt1(e,30); +end +plot(t,h,'b',t,desacE,'r','LineWidth',1.25) +title('Test Results - Desaccaded data') +xlabel('Time in secs') +ylabel('Velocity in deg/sec') +ylim([-400 +400]) +legend ('Head velocity','Eye velocity') +%get positive/neagtive data +limit = size(t); +posH = []; +posE = []; +negH = []; +negE = []; +for n = 1:limit + if h(n) > 0 + posH = vertcat(posH,h(n)); + if desacE(n) > 0 + posE = vertcat(posE,desacE(n)); + else + posE = vertcat(posE,0); + end + end + if h(n) < 0 + negH = vertcat(negH,h(n)); + if desacE(n) < 0 + negE = vertcat(negE,desacE(n)); + else + negE = vertcat(negE,0); + end + end +end +%Positive-Negative plots +subplot(3,2,5) +plot([posH,posE],'LineWidth',1.25) +title('Positive vs Negative data - Desaccaded data') +xlabel('Time in samples') +ylabel('Velocity in deg/sec') +hold on +plot([negH,negE],'LineWidth',1.25) +hold off +legend ('Head velocity','Eye velocity','Head velocity','Eye velocity') +%XY ploy & regresion line +subplot(3,2,6) +hold on +if isOctave + scatter(negH,negE,'c','.'); + scatter(posH,posE,'b','.'); +else + scatter(negH,negE,'.b'); + scatter(posH,posE,'.','MarkerEdgeColor',[0 .7 .7]); +end +title('Head vs Eye plot - Desaccaded data') +xlabel('Head Velocity in deg/sec') +ylabel('Eye Velocity in deg/sec') +negB = negH\negE; +posB = posH\posE; +calcNegE = negB*negH; +calcPosE = posB*posH; +plot(negH,calcNegE,posH,calcPosE,'LineWidth',5) +if s == 1 + plot(negH,negH,'r',posH,posH,'r','LineWidth',1.5) + legend ('Negative data','Positive data','Negative regresion','Positive regresion','No Suppression line','Location','eastoutside') +else + plot(negH,negH,'g',posH,posH,'g','LineWidth',1.5) + legend ('Negative data','Positive data','Negative regresion','Positive regresion','Normality line','Location','eastoutside') +end +hold off +%axis square +%Fourier plots +subplot(3,2,3) +[fHead,P1Head] = fourier(h); +[fEye,P1Eye] = fourier(e); +hold on +stem(fHead,P1Head,'b'); +title('Single-Side Amplitude Spectrum of Head and Eye - RAW data') +xlabel('f (Hz)') +ylabel('|P1(f)|') +xlim([0 5]) +stem(fEye,P1Eye,'r'); +legend('Head','Eye') +hold off +%Desacade Fourier +subplot(3,2,4) +[fHead,P1Head] = fourier(h); +[fEye,P1Eye] = fourier(desacE); +hold on +stem(fHead,P1Head,'b'); +title('Single-Side Amplitude Spectrum of Head and Eye - Desaccaded data') +xlabel('f (Hz)') +ylabel('|P1(f)|') +xlim([0 5]) +stem(fEye,P1Eye,'r'); +legend('Head','Eye') +hold off +%Numerical data +aucPosEye = trapz(posE); +aucPosHead = trapz(posH); +aucNegEye = trapz(negE); +aucNegHead = trapz(negH); +gainPos = (aucPosEye/aucPosHead); +gainNeg = (aucNegEye/aucNegHead); +velGainPos = mean(findpeaks(posE))/mean(findpeaks(posH)); +velGainNeg = mean(findpeaks(abs(negE)))/mean(findpeaks(abs(negH))); +prePosPeakE = mean(findpeaks(posE)); +preNegPeakE = mean(findpeaks(abs(negE))); +if prePosPeakE > preNegPeakE + peakE = prePosPeakE; + peakH = mean(findpeaks(posH)); +else + peakE = -preNegPeakE; + peakH = -mean(findpeaks(abs(negH))); +end +result = strcat('+Area gain: ',num2str(gainPos),' -Area gain: ',num2str(gainNeg),' || +Slope: ',num2str(posB),' -Slope: ',num2str(negB),' || +Velocity gain: ',num2str(velGainPos),' -Velocity gain: ',num2str(velGainNeg),' || Head Maxs(deg/s): ', num2str(peakH),' Eye Maxs: ',num2str(peakE),' || DCoff.: ', num2str(velGainPos-velGainNeg)); +display(result) +if ~isOctave + mTextBox = uicontrol('style','text'); + set(mTextBox,'String',result) + set(mTextBox,'FontSize',12) + set(mTextBox,'Position',[20 5 1300 25]) + set(figure1,'MenuBar','figure') +end + diff --git a/fourier.m b/fourier.m new file mode 100644 index 0000000..38f1419 --- /dev/null +++ b/fourier.m @@ -0,0 +1,9 @@ +function [f,P1] = fourier(data) +Fs = 250; % Sampling frequency +[L,~] = size(data); % Length of signal +X = data; % Data +Y = fft(X); +P2 = abs(Y/L); +P1 = P2(1:fix(L/2)+1); +P1(2:end-1) = 2*P1(2:end-1); +f = Fs*(0:(L/2))/L; diff --git a/isOctave.m b/isOctave.m new file mode 100644 index 0000000..307ad77 --- /dev/null +++ b/isOctave.m @@ -0,0 +1,12 @@ +%% +%% Return: true if the environment is Octave. +%% +function retval = isOctave + persistent cacheval; % speeds up repeated calls + + if isempty (cacheval) + cacheval = (exist ('OCTAVE_VERSION', 'builtin') > 0); + end + + retval = cacheval; +end \ No newline at end of file diff --git a/read.m b/read.m new file mode 100644 index 0000000..6d3253c --- /dev/null +++ b/read.m @@ -0,0 +1,95 @@ +function [t,e,h] = read(s) +t = []; +e = []; +h = []; +[file,directory] = uigetfile('.csv','Select "ASCII results File" to import'); +if file == 0 + return +end +fullArchive = fullfile(directory,file); +readFile = fopen(fullArchive); +tline = fgetl(readFile); +readnow = 0; +timeData = ''; +eyeData = ''; +headData = ''; +while ischar(tline) + + if readnow == 3 + timeData = tline; + end + if readnow == 4 + eyeData = tline; + end + if readnow == 5 + headData = tline; + break + end + if s == 1 + if strcmp(tline,'Test Type;VORS — Horizontal')||strcmp(tline,'Test Type;SRVO — Horizontal')||strcmp(tline,'Tipo de Prueba;SRVO — Horizontal') + readnow = 1; + end + else + if strcmp(tline,'Test Type;VVOR — Horizontal')||strcmp(tline,'Tipo de Prueba;RVVO — Horizontal')||strcmp(tline,'Test Type;RVVO — Horizontal') + readnow = 1; + end + end + + if readnow > 0 + readnow = readnow + 1; + end + tline = fgetl(readFile); +end +fclose(readFile); +timeData = strrep(timeData, ',', '.'); +timeData = strsplit(timeData,';'); +timeData = timeData(2:end); +eyeData = strrep(eyeData, ',', '.'); +eyeData = strsplit(eyeData,';'); +eyeData = eyeData(2:end); +headData = strrep(headData, ',', '.'); +headData = strsplit(headData,';'); +headData = headData(2:end); +[~,b] = size(timeData); +for n = 1:b-1 + t = vertcat(t,str2double(timeData{1,n})); + e = vertcat(e,str2double(eyeData{1,n})); + h = vertcat(h,str2double(headData{1,n})); +end +question = questdlg('Do you want to select a time period or do you want to use entrie data?','Select time window','SET TIME RANGE','USE ALL DATA','USE ALL DATA'); +if strcmp(question,'SET TIME RANGE') + prompt = {'Enter min value of seconds:','Enter max value of seconds:'}; + dlg_title = 'Input range'; + num_lines = 1; + defaultans = {'',''}; + answer = inputdlg(prompt,dlg_title,num_lines,defaultans); + validation = 1; + [a,~] = size(answer); + if a < 2 + display('Max and/or Min are/is missing') + validation = 0; + end + minim = str2double(answer{1,1}); + maxim = str2double(answer{2,1}); + if isnan(minim)||isnan(maxim) + display('Data must to be numeric. Decimal separator is "."') + validation = 0; + end + if maxim < minim + display('Max is < than Min') + validation = 0; + end + if validation == 1 + [~,minpos] = min(abs(t-minim)); + [~,maxpos] = min(abs(t-maxim)); + t = t(minpos:maxpos); + e = e(minpos:maxpos); + h = h(minpos:maxpos); + end +end + + + + + + diff --git a/saccades.m b/saccades.m new file mode 100644 index 0000000..e78f901 --- /dev/null +++ b/saccades.m @@ -0,0 +1,88 @@ +% Jorge Rey Martinez 2017 +% Inputs +% t = time array +% e = eye velocity array +% h = head velocity array +% s = Boolean, true if register is supresed VVOR false if VVOR +function saccades(t,e,h,s) +if isOctave + scrsz = get(0,'ScreenSize'); +else + scrsz = get(groot,'ScreenSize'); +end +if s == 1 + figure2 = figure('Name','SACCADES VORS PLOT','NumberTitle','off','Position',[1 scrsz(4)/1.2 scrsz(3)/1.55 scrsz(4)/1.7]); +else + figure2 = figure('Name','SACCADES VVOR PLOT','NumberTitle','off','Position',[1 scrsz(4)/1.2 scrsz(3)/1.55 scrsz(4)/1.7]); +end +figure(figure2); +%SACCADES plot +subplot(2,1,1) +if s == 1 + desacE = medfilt1(e,35); +else + desacE = medfilt1(e,30); +end +sacE = e-desacE; +plot(t,h,'b',t,sacE,'r','LineWidth',1.25) +%Calculate Sacades per Second +[peaks,locs] = findpeaks(abs(sacE),'MinPeakHeight',100,'MinPeakDistance',10); +[npeaks,~] = size(peaks); +[last,~]= size(t); +sacSec = (npeaks)/(t(last)-t(1)); +[f,P1] = fourier(h); +[~,ind] = max(P1); +headFrec = f(ind); +headPeriod = 1/headFrec; +sacCyc = (npeaks)/((t(last)-t(1))/headPeriod); +%PlotData +string = ['Test Results - Saccades data (Saccades per second = ' num2str(sacSec) ' - Saccades per cycle = ' num2str(sacCyc) ' at ' num2str(headFrec) 'Hz. )']; +title(string) +xlabel('Time in secs') +ylabel('Velocity in deg/sec') +ylim([-400 +400]) +legend ('Head velocity','Eye velocity') +%Plot Detectes Saccades +%subplot(3,1,2) +%findpeaks(abs(sacE),t,'MinPeakHeight',100,'MinPeakDistance',10); +%ylim([0 +400]); +%title('Detected Saccades Plot') +%xlabel('Time in secs') +%ylabel('Eye Velocity in deg/sec') +%Get sacades on time +headSacc = abs(h(locs)); +[sizeHead,~]= size (headSacc); +pos = 1; +sac0 = 0; +sac50 = 0; +sac100 = 0; +sac150 = 0; +sac200 = 0; +sac250 = 0; +while pos < sizeHead + if headSacc(pos) > 125 + sac250 = sac250 + 1; + elseif headSacc(pos) > 100 + sac200 = sac200 + 1; + elseif headSacc(pos) > 75 + sac150 = sac150 + 1; + elseif headSacc(pos) > 50 + sac100 = sac100 + 1; + elseif headSacc(pos) > 25 + sac50 = sac50 + 1; + elseif headSacc(pos) > 0 + sac0 = sac0 + 1; + end + pos = pos + 1; +end +subplot(2,1,2) +total = sac0+sac50+sac100+sac150+sac200+sac250; +dataX = [(sac0/total)*100,(sac50/total)*100,(sac100/total)*100,(sac150/total)*100,(sac200/total)*100,(sac250/total)*100]; +bar(dataX) +Labels = {'0-25', '25-50', '50-75', '75-100', '100-125', '>125'}; +set(gca, 'XTick', 1:6, 'XTickLabel', Labels); +string = ['Saccade count by head velocity (Total Saccades = ' num2str(npeaks) ')']; +title(string) +xlabel('Head velocity intervals in deg/sec') +ylabel('% of sacaddes') + diff --git a/vvor.m b/vvor.m new file mode 100644 index 0000000..90548d7 --- /dev/null +++ b/vvor.m @@ -0,0 +1,15 @@ +s = 0; +question = questdlg('Do you want to read the data from VVOR or from VVORS (suppressed) test?','Select test to read','VVOR','VORS (suppressed)','VVOR'); +if strcmp(question,'VORS (suppressed)') + s = 1; +end +[t,e,h] = read(s); +if isempty(t) + display('No data was loaded') + return +end +analize(t,e,h,s); +question = questdlg('Do you want to analize saccades?','Extended analysis','YES','NO','NO'); +if strcmp(question,'YES') +saccades(t,e,h,s) +end