[chronojump] MIF. Starting of the test uses the RFD



commit f3d3de94fe6ee62e49abb8f5d22f486f6251c221
Author: Xavier Padullés <x padulles gmail com>
Date:   Mon Apr 10 15:43:26 2017 +0200

    MIF. Starting of the test uses the RFD

 r-scripts/maximumIsometricForce.R |  190 +++++++++++++++++++++++--------------
 1 files changed, 117 insertions(+), 73 deletions(-)
---
diff --git a/r-scripts/maximumIsometricForce.R b/r-scripts/maximumIsometricForce.R
index cf99a3f..bb1d8c9 100644
--- a/r-scripts/maximumIsometricForce.R
+++ b/r-scripts/maximumIsometricForce.R
@@ -27,6 +27,7 @@ prepareGraph <- function(os, pngFile, width, height)
                Cairo(width, height, file = pngFile, type="png", bg="white")
        else
                png(pngFile, width=width, height=height)
+               #pdf(file = "/tmp/maxIsomForce.pdf", width=width, height=height)
 }
 
 endGraph <- function()
@@ -53,8 +54,8 @@ assignOptions <- function(options) {
                    rfd200.raw          = as.numeric(options[15]),
                    rfd0_200.raw        = as.numeric(options[16]),
                    rfd0_200.fitted     = as.numeric(options[17]),
-                   rfd50fmax.raw       = as.numeric(options[18]),
-                   rfd50fmax.fitted    = as.numeric(options[19])
+                   rfd50pfmax.raw      = as.numeric(options[18]),
+                   rfd50pfmax.fitted   = as.numeric(options[19])
                    ))
 }
 
@@ -92,40 +93,46 @@ getForceModel <- function(time, force, startTime, # startTime is the instant whe
 getDynamicsFromLoadCellFile <- function(inputFile, averageLength = 0.1, percentChange = 5)
 {
         originalTest = read.csv(inputFile, header = F, dec = ".", sep = ";", skip = 2)
-        colnames(originalTest) <- c("Time", "Force")
-        originalTest$Time = as.numeric(originalTest$Time)
+        colnames(originalTest) <- c("time", "force")
+        originalTest$time = as.numeric(originalTest$time)
+        
+        #Instantaneous RFD
+        rfd = getRFD(originalTest)
         
         #Finding the decrease of the foce to detect the end of the maximum voluntary force
-        trimmingSamples = getTrimmingSamples(originalTest, averageLength = averageLength, percentChange = 
percentChange)
+        trimmingSamples = getTrimmingSamples(originalTest, rfd, averageLength = averageLength, percentChange 
= percentChange)
         startSample = trimmingSamples$startSample
-        startTime = originalTest$Time[startSample]
+        startTime = originalTest$time[startSample]
         
         endSample = trimmingSamples$endSample
-        endTime = originalTest$Time[endSample]
+        endTime = originalTest$time[endSample]
         
         # Initial force. It is needed to perform an initial steady force to avoid jerks and great peaks in 
the force
-        initf = mean(originalTest$Force[1:(startSample - 10)]) #ATENTION. This value is different from f0.raw
-        fmax.raw = max(originalTest$Force)
+        initf = mean(originalTest$force[1:(startSample - 10)]) #ATENTION. This value is different from f0.raw
+        fmax.raw = max(originalTest$force)
         
         #Trimming the data before and after contraction
+        print(paste("StartSample :", startSample))
+        print(paste("EndSample :", endSample))
         test = originalTest[startSample:endSample,]
 
         f.smoothed = getMovingAverageForce(originalTest, averageLength = averageLength) #Running average 
with equal weight averageLength seconds
         fmax.smoothed = max(f.smoothed, na.rm =T)
         
-        model = getForceModel(test$Time, test$Force, startTime, fmax.smoothed, initf)
-        f.fitted = initf + model$fmax*(1-exp(-model$K*(originalTest$Time - startTime)))
+        model = getForceModel(test$time, test$force, startTime, fmax.smoothed, initf)
+        f.fitted = initf + model$fmax*(1-exp(-model$K*(originalTest$time - startTime)))
+ 
         
-        f0.raw = test$Force[1]                                                       #Force at t=0ms. 
ATENTION. This value is different than initf
-        f100.raw = test$Force[which.min(abs(test$Time - (startTime + 0.1)))]         #Force at t=100ms
+        f0.raw = test$force[1]                                                       #Force at t=0ms. 
ATENTION. This value is different than initf
+        f100.raw = test$force[which.min(abs(test$time - (startTime + 0.1)))]         #Force at t=100ms
         f100.fitted = model$fmax*(1 - exp( - model$K * 0.1)) + initf
-        f200.raw = test$Force[which.min(abs(test$Time - (startTime + 0.2)))]         #Force at t=200ms
+        f200.raw = test$force[which.min(abs(test$time - (startTime + 0.2)))]         #Force at t=200ms
         f200.fitted = model$fmax*(1 - exp( - model$K * 0.2)) + initf
-        f90.raw = test$Force[which.min(abs(test$Time - (startTime + 0.09)))]        #Force at t=90ms
-        f110.raw = test$Force[which.min(abs(test$Time - (startTime + 0.11)))]       #Force at t=110ms
-        f190.raw = test$Force[which.min(abs(test$Time - (startTime + 0.19)))]       #Force at t=190ms
-        f210.raw = test$Force[which.min(abs(test$Time - (startTime + 0.21)))]       #Force at t=210ms
-        f50fmax.raw = test$Force[which.min(abs(test$Force - fmax.raw/2))]
+        f90.raw = test$force[which.min(abs(test$time - (startTime + 0.09)))]         #Force at t=90ms
+        f110.raw = test$force[which.min(abs(test$time - (startTime + 0.11)))]        #Force at t=110ms
+        f190.raw = test$force[which.min(abs(test$time - (startTime + 0.19)))]        #Force at t=190ms
+        f210.raw = test$force[which.min(abs(test$time - (startTime + 0.21)))]        #Force at t=210ms
+        f50fmax.raw = test$force[which.min(abs(test$force - fmax.raw/2))]
         f50fmax.fitted = (model$fmax + initf)/2
         
         rfd0.fitted = model$fmax * model$K # RFD at t=0ms using the exponential model
@@ -137,14 +144,15 @@ getDynamicsFromLoadCellFile <- function(inputFile, averageLength = 0.1, percentC
         rfd0_200.fitted = (f200.fitted - initf) / 0.2
         rfd100_200.raw = ((f200.raw - f100.raw) * 10)
         
-        t50fmax.raw = test$Time[which.min(abs(test$Force - fmax.raw/2))]
-        t50fmax.fitted = originalTest$Time[which.min(abs(f.fitted - (model$fmax + initf)/ 2))]
-        rfd50fmax.raw = (test$Force[which.min(abs(test$Force - fmax.raw/2)) + 1] - 
test$Force[which.min(abs(test$Force - fmax.raw/2)) - 1]) * 100 / 2 #RFD at the moment that the force is 50% 
of the fmax
-        rfd50fmax.fitted = model$fmax * model$K * exp(-model$K*(t50fmax.fitted - startTime))
+        t50fmax.raw = test$time[which.min(abs(test$force - fmax.raw/2))]
+        t50fmax.fitted = originalTest$time[which.min(abs(f.fitted - (model$fmax + initf)/ 2))]
+        rfd50pfmax.raw = (test$force[which.min(abs(test$force - fmax.raw/2)) + 1] - 
test$force[which.min(abs(test$force - fmax.raw/2)) - 1]) * 100 / 2 #RFD at the moment that the force is 50% 
of the fmax
+        rfd50pfmax.fitted = model$fmax * model$K * exp(-model$K*(t50fmax.fitted - startTime))
         
-        return(list(nameOfFile = inputFile, time = originalTest[, "Time"], fmax.fitted = model$fmax, 
k.fitted = model$K,
+        return(list(nameOfFile = inputFile, time = originalTest[, "time"], fmax.fitted = model$fmax, 
k.fitted = model$K,
                     startTime = startTime, entTime = endTime,
                     startSample = startSample, endSample = endSample,
+                    totalSample = length(originalTest$time),
                     initf = initf,
                     f0.raw = f0.raw,
                     fmax.raw = fmax.raw, fmax.smoothed = fmax.smoothed,
@@ -161,9 +169,10 @@ getDynamicsFromLoadCellFile <- function(inputFile, averageLength = 0.1, percentC
                     rfd0_200.raw = rfd0_200.raw,
                     rfd0_200.fitted = rfd0_200.fitted,
                     rfd200.raw = rfd200.raw,
-                    t50fmax.raw = t50fmax.raw, rfd50fmax.raw = rfd50fmax.raw,
-                    t50fmax.fitted = t50fmax.fitted, rfd50fmax.fitted = rfd50fmax.fitted,
-                    f.raw = originalTest$Force, f.smoothed = f.smoothed, f.fitted = f.fitted,
+                    rfd = rfd,
+                    t50fmax.raw = t50fmax.raw, rfd50pfmax.raw = rfd50pfmax.raw,
+                    t50fmax.fitted = t50fmax.fitted, rfd50pfmax.fitted = rfd50pfmax.fitted,
+                    f.raw = originalTest$force, f.smoothed = f.smoothed, f.fitted = f.fitted,
                     endTime = endTime))
 }
 
@@ -172,21 +181,25 @@ drawDynamicsFromLoadCell <- function(
                                     dynamics, vlineT0=T, vline50fmax.raw=F, vline50fmax.fitted=F,
                                     hline50fmax.raw=F, hline50fmax.fitted=F,
                                      rfd0.fitted=F, rfd100.raw=F, rfd0_100.raw=F, rfd0_100.fitted = F, 
rfd200.raw=F, rfd0_200.raw=F,
-                                    rfd0_200.fitted = F, rfd50fmax.raw=F, rfd50fmax.fitted=F,
+                                    rfd0_200.fitted = F, rfd50pfmax.raw=F, rfd50pfmax.fitted=F,
                                      xlimits=NA)
 {
-#        pdf(file = pdfFilename, paper="a4r", width = 12, height = 12)
+        par(mar = c(5, 4, 6, 4))
         
         #Plotting raw data from startTime to endTime (Only the analysed data)
         yHeight = dynamics$fmax.raw * 1.1
         if (!is.na(xlimits[1])){
                 xWidth = xlimits[2] - xlimits[1]
                 plot(dynamics$time[dynamics$startSample:dynamics$endSample] , 
dynamics$f.raw[dynamics$startSample:dynamics$endSample],
-                     type="l", xlab="Time[s]", ylab="Force[N]", xlim = xlimits, ylim=c(0, yHeight),
+                     type="l", xlab="Time[s]", ylab="Force[N]",
+                     xlim = xlimits, ylim=c(0, yHeight),
                      main = dynamics$nameOfFile, yaxs= "i", xaxs = "i")
+                xmin = xlimits[1]
+                xmax = xlimits[2]
+                points(dynamics$time[dynamics$startSample:dynamics$endSample] , 
dynamics$f.raw[dynamics$startSample:dynamics$endSample])
         } else if (is.na(xlimits[1])){
-                xmin = dynamics$startTime*0.97 - dynamics$time[1]*0.03
-                xmax = dynamics$endTime*1.1 - dynamics$startTime*0.1
+                xmin = 0
+                xmax = min(c(dynamics$endTime*1.1 - dynamics$startTime*0.1, dynamics$t0 + 1))
                 xWidth = xmax - xmin
                 plot(dynamics$time[dynamics$startSample:dynamics$endSample] , 
dynamics$f.raw[dynamics$startSample:dynamics$endSample],
                      type="l", xlab="Time[s]", ylab="Force[N]",
@@ -196,18 +209,21 @@ drawDynamicsFromLoadCell <- function(
         }
         
         #Plotting not analysed data
-        lines(dynamics$time[1:dynamics$startSample] , dynamics$f.raw[1:dynamics$startSample], type = "c")
-        lines(dynamics$time[dynamics$endSample: length(dynamics$time)] , dynamics$f.raw[dynamics$endSample: 
length(dynamics$time)], type = "c")
+        lines(dynamics$time[1:dynamics$startSample] , dynamics$f.raw[1:dynamics$startSample], col = "grey") 
#Pre-analysis
+        lines(dynamics$time[dynamics$endSample: dynamics$totalSample] , dynamics$f.raw[dynamics$endSample: 
dynamics$totalSample], col = "grey") #Post-analysis
         
-        text( x = min(which(dynamics$f.raw == max(dynamics$f.raw))/100), y = dynamics$fmax.raw, labels = 
paste("Fmax = ", round(dynamics$fmax.raw, digits=2), " N", sep=""), pos = 3)
+        text( x = min(which(dynamics$f.raw == max(dynamics$f.raw))/100), y = dynamics$fmax.raw,
+              labels = paste("Fmax = ", round(dynamics$fmax.raw, digits=2), " N", sep=""), pos = 3)
         
         #Plotting fitted data
         lines(dynamics$time, dynamics$f.fitted, col="blue")
-        text(x = dynamics$time[length(dynamics$time)], y = dynamics$fmax.fitted - 30, labels = paste("Fmax 
=", round(dynamics$fmax.fitted, digits = 2), "N"), pos = 4, col="blue")
-        axis(4, at = dynamics$fmax.fitted + dynamics$initf, labels = round(dynamics$fmax.fitted + 
dynamics$initf, digits = 2), col = "blue")
+        text(x = dynamics$time[dynamics$totalSample], y = dynamics$fmax.fitted - 30,
+             labels = paste("Fmax =", round(dynamics$fmax.fitted, digits = 2), "N"), pos = 4, col="blue")
+        axis(2, at = dynamics$fmax.fitted + dynamics$initf, labels = round(dynamics$fmax.fitted + 
dynamics$initf, digits = 2),
+             line = 2, col = "blue")
         
         #Plottting smoothed data
-        lines(dynamics$time, dynamics$f.smoothed, col="grey")
+        #lines(dynamics$time, dynamics$f.smoothed, col="grey")
         
         if(vlineT0){
                 abline(v = dynamics$startTime, lty = 2)
@@ -245,7 +261,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercepstartTime.fitted, b = dynamics$rfd0.fitted, col="blue", lty = 2)
                 text(x = (yHeight - 5 - intercepstartTime.fitted)/dynamics$rfd0.fitted, y = yHeight - 5,
                      label=paste("RFD0 =", round(dynamics$rfd0.fitted, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75, col = "blue")
+                     srt=atan(windowSlope)*180/pi, pos = 2, col = "blue")
                 points(x = dynamics$startTime, y = dynamics$initf, col = "blue")
         }
         
@@ -256,7 +272,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercept100.raw, b = dynamics$rfd100.raw, lty = 2)
                 text(x = (yHeight - 5 - intercept100.raw)/dynamics$rfd100.raw, y = yHeight - 5,
                      label=paste("RFD100 =", round(dynamics$rfd100, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75)
+                     srt=atan(windowSlope)*180/pi, pos = 2)
                 points(x = dynamics$startTime + 0.1, y = dynamics$f100.raw)
         }
         
@@ -267,7 +283,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercepstartTime_100.raw, b = dynamics$rfd0_100.raw, lty = 2)
                 text(x = (yHeight - 5 - intercepstartTime_100.raw)/dynamics$rfd0_100.raw, y = yHeight - 5,
                      label=paste("RFD0_100 =", round(dynamics$rfd0_100.raw, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75)
+                     srt=atan(windowSlope)*180/pi, pos = 2)
                 points(x = c(dynamics$startTime, dynamics$startTime + 0.1), y = c(dynamics$f0.raw, 
dynamics$f100.raw))
         }
         
@@ -278,7 +294,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercepstartTime_100.fitted, b = dynamics$rfd0_100.fitted, lty = 2, col = "blue")
                 text(x = (yHeight - 5 - intercepstartTime_100.fitted)/dynamics$rfd0_100.fitted, y = yHeight 
- 5,
                      label=paste("RFD0_100 =", round(dynamics$rfd0_100.fitted, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75, col = "blue")
+                     srt=atan(windowSlope)*180/pi, pos = 2, col = "blue")
                 points(x = c(dynamics$startTime, dynamics$startTime + 0.1), y = c(dynamics$initf, 
dynamics$f100.fitted), col = "blue")
         }
         
@@ -289,7 +305,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercept200.raw, b = dynamics$rfd200.raw, lty = 2)
                 text(x = (yHeight - 5 - intercept200.raw)/dynamics$rfd200.raw, y = yHeight - 5,
                      label=paste("RFD200 =", round(dynamics$rfd200.raw, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75)
+                     srt=atan(windowSlope)*180/pi, pos = 2)
                 points(x = dynamics$startTime + 0.2, y = dynamics$f200.raw)
         }
         if(rfd0_200.raw)
@@ -299,7 +315,7 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercepstartTime_200.raw, b = dynamics$rfd0_200.raw, lty = 2)
                 text(x = (yHeight - 5 - intercepstartTime_200.raw)/dynamics$rfd0_200.raw, y = yHeight - 5,
                      label=paste("RFD0_200 =", round(dynamics$rfd0_200.raw, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75)
+                     srt=atan(windowSlope)*180/pi, pos = 2)
                 points(x = c(dynamics$startTime, dynamics$startTime + 0.2), y = 
c(dynamics$f0.raw,dynamics$f200.raw))
         }
         
@@ -310,30 +326,41 @@ drawDynamicsFromLoadCell <- function(
                 abline(a = intercepstartTime_200.fitted, b = dynamics$rfd0_200.fitted, lty = 2, col = "blue")
                 text(x = (yHeight - 5 - intercepstartTime_200.fitted)/dynamics$rfd0_200.fitted, y = yHeight 
- 5,
                      label=paste("RFD0_200 =", round(dynamics$rfd0_200.fitted, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75, col = "blue")
+                     srt=atan(windowSlope)*180/pi, pos = 2, col = "blue")
                 points(x = c(dynamics$startTime, dynamics$startTime + 0.2), y = 
c(dynamics$initf,dynamics$f200.fitted), col = "blue")
         }
         
-        if(rfd50fmax.raw){
-                windowSlope = dynamics$rfd50fmax.raw*(plotHeight/yHeight)/(plotWidth/xWidth)
-                intercept50fmax.raw = dynamics$f50fmax.raw - dynamics$rfd50fmax.raw * dynamics$t50fmax.raw
-                abline(a = intercept50fmax.raw, b = dynamics$rfd50fmax.raw, lty = 2)
-                text(x = (yHeight - 5 - intercept50fmax.raw)/dynamics$rfd50fmax.raw, y = yHeight - 5,
-                     label=paste("RFD50fmax =", round(dynamics$rfd50fmax.raw, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos = 2, cex = 0.75)
+        if(rfd50pfmax.raw){
+                windowSlope = dynamics$rfd50pfmax.raw*(plotHeight/yHeight)/(plotWidth/xWidth)
+                intercept50fmax.raw = dynamics$f50fmax.raw - dynamics$rfd50pfmax.raw * dynamics$t50fmax.raw
+                abline(a = intercept50fmax.raw, b = dynamics$rfd50pfmax.raw, lty = 2)
+                text(x = (yHeight - 5 - intercept50fmax.raw)/dynamics$rfd50pfmax.raw, y = yHeight - 5,
+                     label=paste("rfd50%fmax =", round(dynamics$rfd50pfmax.raw, digits=0), "N/s"),
+                     srt=atan(windowSlope)*180/pi, pos = 2)
                 points(x =dynamics$t50fmax.raw, y = dynamics$f50fmax.raw)
         }
         
-        if(rfd50fmax.fitted){
-                windowSlope = dynamics$rfd50fmax.fitted*(plotHeight/yHeight)/(plotWidth/xWidth)
-                intercept50fmax.fitted = dynamics$f50fmax.fitted - dynamics$rfd50fmax.fitted * 
dynamics$t50fmax.fitted
-                abline(a = intercept50fmax.fitted, b = dynamics$rfd50fmax.fitted, col ="blue", lty = 2)
-                text(x = (yHeight - 5 - intercept50fmax.fitted)/dynamics$rfd50fmax.fitted, y = yHeight - 5, 
col = "blue",
-                     label=paste("RFD50fmax =", round(dynamics$rfd50fmax.fitted, digits=0), "N/s"),
-                     srt=atan(windowSlope)*180/pi, pos=2, cex = 0.75)
+        if(rfd50pfmax.fitted){
+                windowSlope = dynamics$rfd50pfmax.fitted*(plotHeight/yHeight)/(plotWidth/xWidth)
+                intercept50fmax.fitted = dynamics$f50fmax.fitted - dynamics$rfd50pfmax.fitted * 
dynamics$t50fmax.fitted
+                abline(a = intercept50fmax.fitted, b = dynamics$rfd50pfmax.fitted, col ="blue", lty = 2)
+                text(x = (yHeight - 5 - intercept50fmax.fitted)/dynamics$rfd50pfmax.fitted, y = yHeight - 5, 
col = "blue",
+                     label=paste("rfd50%fmax =", round(dynamics$rfd50pfmax.fitted, digits=0), "N/s"),
+                     srt=atan(windowSlope)*180/pi, pos=2)
                 points(x = dynamics$t50fmax.fitted, y = dynamics$f50fmax.fitted, col = "blue")
         }
-        #dev.off()
+        
+        #Plotting instantaneous RFD
+        par(new = T)
+        print(dynamics$time)
+        print(dynamics$f.raw)
+        print(dynamics$rfd)
+        plot(dynamics$time[(dynamics$startSample):dynamics$endSample], 
dynamics$rfd[(dynamics$startSample):dynamics$endSample],
+             type = "l", col = "red", axes = F,
+             xlim = c(xmin, xmax), xlab = "", ylab = "", yaxs= "i", xaxs = "i")
+        lines(dynamics$time[1:(dynamics$startSample)], dynamics$rfd[1:dynamics$startSample], col = "orange")
+        lines(dynamics$time[dynamics$endSample:dynamics$totalSample], 
dynamics$rfd[dynamics$endSample:dynamics$totalSample], col = "orange")
+        axis(4, col = "red", line = 1)
 }
 
 getDynamicsFromLoadCellFolder <- function(folderName, resultFileName, export2Pdf)
@@ -344,7 +371,7 @@ getDynamicsFromLoadCellFolder <- function(folderName, resultFileName, export2Pdf
         colnames(results)=c("fileName", "fmax.fitted", "k.fitted", "fmax.raw", "startTime", "initf", 
"fmax.smoothed",
                             "rfd0.fitted", "rfd100.raw", "rfd0_100.raw", "rfd0_100.fitted",
                             "rfd200.raw", "rfd0_200.raw", "rfd0_200.fitted",
-                            "rfd50fmax.raw", "rfd50fmax.fitted")
+                            "rfd50pfmax.raw", "rfd50pfmax.fitted")
         
         results[,"fileName"] = originalFiles
         
@@ -366,8 +393,8 @@ getDynamicsFromLoadCellFolder <- function(folderName, resultFileName, export2Pdf
                 results[i, "rfd200.raw"] = dynamics$rfd200.raw
                 results[i, "rfd0_200.raw"] = dynamics$rfd0_200.raw
                 results[i, "rfd0_200.fitted"] = dynamics$rfd0_200.fitted
-                results[i, "rfd50fmax.raw"] = dynamics$rfd50fmax.rawfilter(test$Force, rep(1/19, 19), sides 
= 2)
-                results[i, "rfd50fmax.fitted"] = dynamics$rfd50fmax.fitted
+                results[i, "rfd50pfmax.raw"] = dynamics$rfd50pfmax.rawfilter(test$force, rep(1/19, 19), 
sides = 2)
+                results[i, "rfd50pfmax.fitted"] = dynamics$rfd50pfmax.fitted
         }
         write.table(results, file = resultFileName, sep = ";", dec = ",", col.names = NA)
         return(results)
@@ -376,31 +403,48 @@ getDynamicsFromLoadCellFolder <- function(folderName, resultFileName, export2Pdf
 
 #Finds the sample in which thpercentChangee force decrease a given percentage of the maximum force.
 #The maximum force is calculed from the moving average of averageLength seconds
-getTrimmingSamples <- function(test, movingAverageForce, averageLength = 0.1, percentChange = 5)
+getTrimmingSamples <- function(test, rfd, movingAverageForce, averageLength = 0.1, percentChange = 5)
 {
         movingAverageForce = getMovingAverageForce(test, averageLength = 0.1)
         maxAverageForce = max(movingAverageForce, na.rm = T)
         maxSample = min(which(movingAverageForce == maxAverageForce), na.rm = T)
+        maxRFD = max(rfd[2:(length(rfd) - 1)])
         
-        #Detecting an increase of percentChange% of the maximum force in the initial signal
-        startSample = 1
-        while((test$Force[startSample] + test$Force[startSample +2]) / 2 <= 
maxAverageForce*1.5*percentChange/100)
+        #Detecting an RFD 
+        startSample = 2
+        
+        while(rfd[startSample] <= maxRFD*0.2)
         {
                 startSample = startSample + 1
-        } 
+        }
         
         #Detecting a decrease of percentChange% in the maximum force
         endSample = min(which((movingAverageForce < maxAverageForce*(100 - percentChange) / 100 &
-                                       test$Time > test$Time[maxSample])), na.rm = T)
+                                       test$time > test$time[maxSample])), na.rm = T)
+        if(is.infinite(endSample))
+        {
+                endSample = length(test$time)
+        }
         
         return(list(startSample = startSample, endSample = endSample))
 }
 
+getRFD <- function(test)
+{
+        #Instantaneous RFD
+        rfd = rep(NA, length(test$time))
+        print(paste("RFD length :",length(test$time)))
+        for (n in 2:(length(test$time) - 1))
+        {
+                rfd[n] = (test$force[n + 1] - test$force[n - 1])/(test$time[n + 1] - test$time[n - 1])
+        }
+        return(rfd)
+}
 getMovingAverageForce <- function(test, averageLength = 0.1)
 {
-        sampleRate = (length(test$Time) - 1) / (test$Time[length(test$Time)] - test$Time[1])
+        sampleRate = (length(test$time) - 1) / (test$time[length(test$time)] - test$time[1])
         lengthSamples = round(averageLength * sampleRate, digits = 0)
-        movingAverageForce = filter(test$Force, rep(1/lengthSamples, lengthSamples), sides = 2)
+        movingAverageForce = filter(test$force, rep(1/lengthSamples, lengthSamples), sides = 2)
         return(movingAverageForce)
 }
 
@@ -408,11 +452,11 @@ prepareGraph(op$os, pngFile, op$graphWidth, op$graphHeight)
 dynamics = getDynamicsFromLoadCellFile(dataFile, op$averageLength, op$percentChange)
 drawDynamicsFromLoadCell(dynamics, op$vlineT0, op$vline50fmax.raw, op$vline50fmax.fitted, 
op$hline50fmax.raw, op$hline50fmax.fitted,
                          op$rfd0.fitted, op$rfd100.raw, op$rfd0_100.raw, op$rfd0_100.fitted, op$rfd200.raw, 
op$rfd0_200.raw, op$rfd0_200.fitted,
-                         op$rfd50fmax.raw, op$rfd50fmax.fitted)
+                         op$rfd50pfmax.raw, op$rfd50pfmax.fitted)
 endGraph()
 
 #dynamics = getDynamicsFromLoadCellFile("~/ownCloud/Xavier/Recerca/Yoyo-Tests/Galga/RowData/APl1", 
averageLength = 0.1, percentChange = 5, sep = ";", dec = ",")
 #drawDynamicsFromLoadCell(dynamics, vlineT0=F, vline50fmax.raw=F, vline50fmax.fitted=T, hline50fmax.raw=F, 
hline50fmax.fitted=T, 
 #                         rfd0.fitted=T, rfd100.raw=F, rfd0_100.raw=F, rfd0_100.fitted = F, rfd200.raw=F, 
rfd0_200.raw=F, rfd0_200.fitted = F,
-#                         rfd50fmax.raw=F, rfd50fmax.fitted=T)
+#                         rfd50pfmax.raw=F, rfd50pfmax.fitted=T)
 #getDynamicsFromLoadCellFolder("~/Documentos/RowData/", resultFileName = "~/Documentos/results.csv")


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]