PR score minimum head and new gain Fourier-based callculation method

analize.m

@@ -66,15 +66,65 @@ function analize(t,e,h,s)
 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)));
 
 %Eye Head Regression Gain
 negB = negH\negE;
 posB = posH\posE;
 calcNegE = negB*negH;
 calcPosE = posB*posH;
 
+%Fourier based Gain calcullation algorithm
+headThreshold = 20; % Set minimum head velocity for frequency analysis
+[dataEyeR,dataEyeL,dataHeadR,dataHeadL] = splitTest(desacE,h,headThreshold);
+[fHeadL,P1HeadL] = fourier(dataHeadL);
+[fEyeL,P1EyeL] = fourier(dataEyeL);
+[fHeadR,P1HeadR] = fourier(dataHeadR);
+[fEyeR,P1EyeR] = fourier(dataEyeR);
+[maxHeadPwrR,HeadPwrRIndex] = max(P1HeadR);
+[maxHeadPwrL,HeadPwrLIndex] = max(P1HeadL);
+maxHeadRFreq = fHeadR(HeadPwrRIndex);
+maxHeadLFreq = fHeadL(HeadPwrLIndex);
+[preMaxEyeRFreq,maxEyeRFreqIndex]=min(abs(fEyeR-maxHeadRFreq));
+[preMaxEyeLFreq,maxEyeLFreqIndex]=min(abs(fEyeL-maxHeadLFreq));
+maxEyeRFreq = fEyeR(maxEyeRFreqIndex);
+maxEyeLFreq = fEyeL(maxEyeLFreqIndex);
+maxEyeRPwr = P1EyeR(maxEyeRFreqIndex);
+maxEyeLPwr = P1EyeL(maxEyeLFreqIndex);
+leftFouGain = (maxEyeLPwr/maxHeadPwrL);
+rightFouGain = (maxEyeRPwr/maxHeadPwrR);
+
+%Analysis of head oscillations variability:
+distanciaPicos = 60;
+lHeadPeaks = findpeaks(posH,'MinPeakDistance',distanciaPicos);
+rHeadPeaks = findpeaks(abs(negH),'MinPeakDistance',distanciaPicos);
+velocityTreshold = 25;
+lHeadInvalids = lHeadPeaks<velocityTreshold;
+rHeadInvalids = rHeadPeaks<velocityTreshold;
+lHeadPeaks(lHeadInvalids) = [];
+rHeadPeaks(rHeadInvalids) = [];
+[lPeakN, ~] = size(lHeadPeaks);
+[rPeakN, ~] = size(rHeadPeaks);
+[lPR,rPR,saccadePositions] = prScoreVVR(t,e,h,s);
+
+%PR PLOT only available in VVOR tests
+if s ~= 1
+    subplot(3,2,5)
+    plot(t,h,'b',t,e,'r','LineWidth',1.5)
+    prTitle = strcat('Saccade Recognition & PR Plot || ', '  LEFT PR:',num2str(lPR),',  RIGHT PR: ',num2str(rPR));
+    title(prTitle)
+    xlabel('Time in samples')
+    ylabel('Velocity in deg/sec')
+    ylim([-400 +400])
+    %add saccade detection to plot
+    [sP,~,~] = find(t==saccadePositions);
+    hold on
+    plot(t(sP),e(sP),'ko')
+    hold off
+    legend ('Head velocity','Eye velocity','Detected Saccade')
+end
+
+
+
 %%%%% PLOTS SECTION %%%%%
 
 %RAW plot
@@ -85,26 +135,17 @@ function analize(t,e,h,s)
 ylabel('Velocity in deg/sec')
 ylim([-400 +400])
 legend ('Head velocity','Eye velocity')
+
 %Desaccaded plot
 subplot(3,2,2)
 plot(t,h,'b',t,desacE,'r','LineWidth',1.25)
-AUCTitle = strcat('Test Output - Desaccaded data  - ',' LEFT GAIN: ',num2str(gainPos),' -RIGHT GAIN: ',num2str(gainNeg));
+AUCTitle = strcat('Test Output - Desaccaded data  || ',' LEFT GAIN: ',num2str(gainPos),' - RIGHT GAIN: ',num2str(gainNeg));
 title(AUCTitle)
 xlabel('Time in secs')
 ylabel('Velocity in deg/sec')
 ylim([-400 +400])
 legend ('Head velocity','Eye velocity')
 
-% %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)
@@ -116,7 +157,7 @@ function analize(t,e,h,s)
     scatter(negH,negE,'.b');
     scatter(posH,posE,'.','MarkerEdgeColor',[0 .7 .7]);
 end
-XYTitle = strcat('Head vs Eye plot - Desaccaded data',' LEFT GAIN: ',num2str(posB),' -RIGHT GAIN: ',num2str(negB));
+XYTitle = strcat('Head vs Eye plot - Desaccaded data ||',' LEFT GAIN: ',num2str(posB),' - RIGHT GAIN: ',num2str(negB));
 title(XYTitle)
 xlabel('Head Velocity in deg/sec')
 ylabel('Eye Velocity in deg/sec')
@@ -130,70 +171,92 @@ function analize(t,e,h,s)
 end
 hold off
 %axis square
-%Fourier plots
+
+%Fourier plot
 subplot(3,2,3)
 [fHead,P1Head] = fourier(h);
 [fEye,P1Eye] = fourier(e);
+[~,ixx] = max(P1Head);
+maxFreqHeadFour = fHead(ixx);
 hold on
 stem(fHead,P1Head,'b');
-title('Single-Side Amplitude Spectrum of Head and Eye - RAW data')
+fourierTitle = strcat('Single-Side Amplitude Spectrum of Head and Eye (RAW) || Head Freq(Hz): ',num2str(maxFreqHeadFour));
+title(fourierTitle)
 xlabel('f (Hz)')
 ylabel('|P1(f)|')
 xlim([0 5])
 stem(fEye,P1Eye,'r');
 legend('Head','Eye')
 hold off
-%Desacade Fourier
+
+%Fourier gain plot
 subplot(3,2,4)
-[fHead,P1Head] = fourier(h);
-[fEye,P1Eye] = fourier(desacE);
+axis off
 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')
+plotBar = bar([maxHeadPwrL maxEyeLPwr;maxHeadPwrR maxEyeRPwr],'hist');
+gainFtitle = strcat(['Fourier Gain || LEFT(1): ',num2str(leftFouGain),' - RIGHT(2): ',num2str(rightFouGain)]);
+title(gainFtitle)
+legend('Head','Eye','Location','eastoutside');
+if ~isOctave
+    plotBar(1).FaceColor = 'b';
+    plotBar(2).FaceColor = 'r';
+end
 hold off
 
-
-%Analysis of head oscillations variability:
-distanciaPicos = 60;
-lHeadPeaks = findpeaks(posH,'MinPeakDistance',distanciaPicos);
-rHeadPeaks = findpeaks(abs(negH),'MinPeakDistance',distanciaPicos);
-velocityTreshold = 25;
-lHeadInvalids = lHeadPeaks<velocityTreshold;
-rHeadInvalids = rHeadPeaks<velocityTreshold;
-lHeadPeaks(lHeadInvalids) = [];
-rHeadPeaks(rHeadInvalids) = [];
-[lPeakN, ~] = size(lHeadPeaks);
-[rPeakN, ~] = size(rHeadPeaks);
-[lPR,rPR,saccadePositions] = prScoreVVR(t,e,h,s);
-
-%PR PLOT only available in VVOR tests
-if s ~= 1
-    subplot(3,2,5)
-    plot(t,h,'b',t,e,'r','LineWidth',1.5)
-    prTitle = strcat('Saccade Recognition & PR Plot: ', '  Left PR:',num2str(lPR),',  Right PR: ',num2str(rPR));
-    title(prTitle)
-    xlabel('Time in samples')
-    ylabel('Velocity in deg/sec')
-    ylim([-400 +400])
-    %add saccade detection to plot
-    [sP,~,~] = find(t==saccadePositions);
-    hold on
-    plot(t(sP),e(sP),'ko')
-    hold off
-    legend ('Head velocity','Eye velocity','Detected Saccade')
-end
+% For FourierGain debug purposes only (uncoment next section)
+% fgFigure = figure;
+% subplot(2,2,3)
+% hold on
+% stem(fHeadL,P1HeadL,'b');
+% title('Single-Side Amplitude Spectrum of Head and Eye (Desaccaded) || LEFT')
+% xlabel('f (Hz)')
+% ylabel('|P1(f)|')
+% xlim([0 5])
+% stem(fEyeL,P1EyeL,'r');
+% legend('Head','Eye')
+% hold off
+% 
+% subplot(2,2,4)
+% hold on
+% stem(fHeadR,P1HeadR,'b');
+% title('Single-Side Amplitude Spectrum of Head and Eye (Desaccaded) || RIGHT')
+% xlabel('f (Hz)')
+% ylabel('|P1(f)|')
+% xlim([0 5])
+% stem(fEyeR,P1EyeR,'r');
+% legend('Head','Eye')
+% hold off
+% 
+% subplot(2,2,1)
+% hold on
+% plot(dataHeadL,'b','LineWidth',1.25)
+% plot(dataEyeL,'r','LineWidth',1.25)
+% hold off
+% title('Test Output - LEFT SIDE')
+% xlabel('Samples')
+% ylabel('Velocity in deg/sec')
+% ylim([-400 +400])
+% legend ('Head velocity','Eye velocity')
+% 
+% subplot(2,2,2)
+% hold on
+% plot(dataHeadR,'b','LineWidth',1.25)
+% plot(dataEyeR,'r','LineWidth',1.25)
+% hold off
+% title('Test Output - RIGHT SIDE')
+% xlabel('Samples')
+% ylabel('Velocity in deg/sec')
+% ylim([-400 +400])
+% legend ('Head velocity','Eye velocity')
+% 
+% set(0, 'currentfigure', figure1);
 
 
 %%%%%%%%%Output analysis results to text%%%%%%%%%%%%
 
-resultG = strcat('GAIN RESULTS: ',' Left(area): ',num2str(gainPos),' Right(area): ',num2str(gainNeg),' || Left(slope): ',num2str(posB),' Right(slope): ',num2str(negB),' || Head Max(º/s):  ', num2str(peakH),' Eye Max: ',num2str(peakE));
+resultG = strcat('GAIN RESULTS: ',' Left(area): ',num2str(gainPos),' Right(area): ',num2str(gainNeg),' || Left(slope): ',num2str(posB),' Right(slope): ',num2str(negB),' || Left(Fourier): ',num2str(leftFouGain),' Right(Fourier): ',num2str(rightFouGain),' || Head Max(º/s):  ', num2str(peakH),' Eye Max: ',num2str(peakE));
 if s~= 1
-    resultPR = strcat('PR RESULTS: ',' Left PR Score: ',num2str(lPR),' Right PR score: ',num2str(rPR),' || Left/Right peaks > 25º/s: ',num2str(lPeakN),'/',num2str(rPeakN),' || Left/Right velocity SD of peaks: ',num2str(std(lHeadPeaks)),'/',num2str(std(rHeadPeaks)));
+    resultPR = strcat('PR RESULTS: ',' Left PR Score: ',num2str(lPR),' Right PR score: ',num2str(rPR),' || Left/Right head peaks > 25º/s: ',num2str(lPeakN),'/',num2str(rPeakN),' || Left/Right velocity SD of head peaks: ',num2str(std(lHeadPeaks)),'/',num2str(std(rHeadPeaks)));
 else
     resultPR = 'PR score is not available for VORS - supression - testing';
 end
@@ -212,3 +275,5 @@ function analize(t,e,h,s)
 disp(resultPR);
 end
 
+
+

prScoreVVR.m

@@ -33,7 +33,8 @@
             left = false;
         end
         [preCheck,~] = size(eInt);
-        if preCheck > 3
+        [preCheckSec,~] = max(abs(hInt));
+        if preCheck > 3 && preCheckSec > 15 % Set a minimum eye width of 3 samples and minimum head peak velocity to 15
             if isOctave
                 [peaks,locs] = findpeaks(abs(eInt),'MinPeakHeight',140,'MinPeakDistance',30);
             else
@@ -52,7 +53,7 @@
                 %disp('No peak found') %Debug
             end
         else
-            disp(["Not enought data in iteration: ",num2str(n)])
+            disp(["Not enought data in segment: ",num2str(n)])
         end
         n = n + 1;
     end

splitTest.m

@@ -0,0 +1,35 @@
+function [dataEyeR,dataEyeL,dataHeadR,dataHeadL] = splitTest(e,h,headThreshold)
+sigPos = logical(h > 0);%Determine positive and negatve values
+cros = sigPos - circshift(sigPos,1); %get sign changes
+crosPos = find(cros);
+[a,~] = size(crosPos);
+n = 2;
+dataEyeR = [];
+dataEyeL = [];
+dataHeadR = [];
+dataHeadL = [];
+while n <= a
+    dataHeadInt = h(crosPos(n-1):crosPos(n));
+    dataEyeInt = e(crosPos(n-1):crosPos(n));
+    if cros(crosPos(n-1)) == 1
+        left = true;
+    else
+        left = false;
+    end
+    [preCheck,~] = max(abs(dataHeadInt));
+    if preCheck > headThreshold
+        if left
+            dataEyeL = vertcat(dataEyeL,dataEyeInt);
+            dataEyeL = vertcat(dataEyeL,(-1*dataEyeInt));
+            dataHeadL = vertcat(dataHeadL,dataHeadInt);
+            dataHeadL = vertcat(dataHeadL,(-1*dataHeadInt));
+        else
+            dataEyeR = vertcat(dataEyeR,dataEyeInt);
+            dataEyeR = vertcat(dataEyeR,(-1*dataEyeInt));
+            dataHeadR = vertcat(dataHeadR,dataHeadInt);
+            dataHeadR = vertcat(dataHeadR,(-1*dataHeadInt));
+        end
+    end
+    n = n + 1;
+end
+end