(* Input CSV file volume-time curve *)

inputfile = "inputfile.csv";

data=Import[inputfile][[All,2;;3]];

(* Example data for JOVE paper *)

data = {{0.`,8.72`},{3.31`,8.48`},{6.61`,8.11`},{9.92`,7.97`},{13.23`,7.79`},{16.53`,7.41`},{19.84`,6.7`},{23.15`,5.41`},{26.45`,4.62`},{29.76`,4.11`},{33.07`,3.54`},{36.38`,3.02`},{39.68`,2.3`},{42.99`,1.93`},{46.3`,1.72`},{49.6`,1.47`},{52.91`,1.52`},{56.22`,1.78`},{59.52`,3.02`},{62.83`,3.99`},{66.14`,4.41`},{69.44`,4.77`},{72.75`,5.47`},{76.06`,6.2`},{79.36`,6.81`},{82.67`,7.28`},{85.98`,7.48`},{89.28`,7.79`},{92.59`,8.01`},{95.9`,8.15`},{99.21`,8.35`},{102.51`,8.56`}};

Dimensions[data]

(* SavitzkyGolay Filter *)

ClearAll[SavitzkyGolayHatMatrix];
SavitzkyGolayHatMatrix[windowHalfWidth_Integer,polynomialOrder_Integer,derivativeOrder_Integer,opts___]/;0≤derivativeOrder≤polynomialOrder≤2 windowHalfWidth:=Block[{wprec=WorkingPrecision/.{opts}~Join~Options[SavitzkyGolayHatMatrix,WorkingPrecision],vandermondeMatrix},vandermondeMatrix=Table[If[i==j==0,1,i^j],{i,-windowHalfWidth,+windowHalfWidth},{j,0,polynomialOrder}]~N~wprec;
(-1)^derivativeOrder derivativeOrder! Drop[vandermondeMatrix,None,-derivativeOrder].Drop[PseudoInverse[vandermondeMatrix],derivativeOrder]];
Options[SavitzkyGolayHatMatrix]={WorkingPrecision->Infinity};
Attributes[SavitzkyGolayHatMatrix]={NHoldAll};

ClearAll[SavitzkyGolayFilter];
SavitzkyGolayFilter[data_?VectorQ,windowHalfWidth_Integer,polynomialOrder_Integer,derivativeOrder_Integer,opts___]/;0≤derivativeOrder≤polynomialOrder≤2 windowHalfWidth<Length[data]:=Block[{wprec=WorkingPrecision/.{opts}~Join~Options[SavitzkyGolayFilter,WorkingPrecision],dprec=Precision[data],ndat,hatMatrix},ndat=If[wprec<dprec,N[data,wprec],data];
hatMatrix=SavitzkyGolayHatMatrix[windowHalfWidth,polynomialOrder,derivativeOrder,WorkingPrecision->Min[wprec,dprec]];
Flatten[{ListConvolve[#,Take[ndat,2 windowHalfWidth+1]]&/@Reverse@Take[hatMatrix,-windowHalfWidth],ListConvolve[hatMatrix[[windowHalfWidth+1]],ndat],ListConvolve[#,Take[ndat,-(2 windowHalfWidth+1)]]&/@Reverse@Take[hatMatrix,windowHalfWidth]}]];
Options[SavitzkyGolayFilter]={WorkingPrecision->Infinity};
Attributes[SavitzkyGolayFilter]={NHoldRest};

(* Normalize volume to 100% *)

normfactor=Mean[Sort[data[[All,2]],Greater][[1;;3]]];
datanorm=Table[100*data[[i,2]]/normfactor,{i,1,Length[data]}];
timestep = data[[2,1]]-data[[1,1]];

(* Manually shift the data to begin systole *)

datan=RotateRight[datanorm,-0];                
ListPlot[datan,PlotRange->All]

(* Calculation of the time-derivative dV/dt using Savitzky-Golay filter*)

ddata=SavitzkyGolayFilter[datan,2 (*window half-width*),2(*polynomial order*),1 (*derivative order*)];

(* Peak values *)
peaks=FindPeaks[ddata,0,0,-∞];

ListPlot[ddata,Joined->True,Epilog->{Orange,PointSize[0.03],Point[peaks]}]

(* Manual selection of E-peak and A-peak *)

Epeak =2;
Apeak=3;

(* E' in % / ms *)

EP= peaks[[Epeak,2]] /timestep ;

(* A' in % / ms *)

AP = peaks[[Apeak,2]]/timestep ;

(* FINAL RESULT: E'/A' ratio *)

EA=EP/AP

(* Export data *)

datalist ={
{"E' (% / ms)",EP},
{"A' (% / ms)",AP},
{"E/A ratio",EA},
{"  ","  "}
};

Export[inputfile<>"_parameters.xlsx",datalist];