Calculating number of peaks within window - HELP translating MATLAB script to R script

[sap1403]

Dear all,

I am a R newbie and looking for some help to translate a freely available MATLAB script into R.

The MATLAB script is designed to analyse a 27650x1 vector which consists of a bunch of peaks and troughs. Each element within the vector corresponds to a time of 0.04s (so first element is at 0s, 2nd is at 0.04s, 3rd is 0.08 etc). The MATLAB script is supposed to find the number of peaks in a ten second window, then divides it by 10, then stores that value in a separate vector, then the window should shift along by one second and the process is repeated. Original Link

The proposed MATLAB script, that I want to translate, looks like this:

% fake data
data = rand(27650, 1);

% Set parameters
sampleRate = 0.04; %time interval (sec)
windowSize = 10; %time (seconds), moving window size
windowIncr = 1.0; %time (seconds), window increment
n = length(data); %number of data points
start = 0; %staring point (time, sec)

% Calculate end time
stopTime = n * sampleRate - start;

% create vector of time stamps (this isn't needed but you might want it for plots, etc.)
time = start : sampleRate : stopTime;

% Window index size
windowIdxSz = windowSize / sampleRate; % !! This must be an integer!
if mod(windowIdxSz,1) ~= 0
error('Window size divided by sample rate must be an integer')
end

% initialize window index
windIdx = 1 : windowIdxSz;

% Number of indicies to shift window
windShift = windowIncr/sampleRate; % !! This must be an integer
if mod(windShift,1) ~= 0
error('Window increment divided by sample rate must be an integer')
end

% determine number of loops needed
nLoops = floor((n-start)/windowIdxSz); %floor() results in skipping leftovers outside of final window

% Loop through your data
nPeaks = nan(nLoops, 1);
for i = 1:nLoops
pks = findpeaks(data(windIdx)); % find peaks within your windowed data
nPeaks(i) = length(pks)/10; % store the number of peaks, divided by 10
windIdx = windIdx + windShift; % shift window
end

My attempt to translate this MATLAB script into an R script looks currently like this:

set parameters

sampleRate = 0.04
windowSize = 10
windowIncr = 1.0
n <- length(data)
start = 0

calculate end time

stopTime <- n*sampleRate-start

create vector of time stamps

time <-start:sampleRate:stopTime

window index size

windowIdxSz <- windowSize / sampleRate
if (mod(windowIdxSz,1) != 0)
{error('Window siize divided by samplerate must be an integer')
}

initialize window index

windIdx <- 1:windowIdxSz

Number of indices to shift window

windshift <- windowIncr / sampleRate
if (mod(windshift,1) !=0)
{error('Window increment divided by sample rate must be an integer')
}

determine number of loops needed

nLoops <- floor((n-start)/windowIdxSz)

loop through data

nPeaks <- matrix(nLoops,1)
for (i in 1:nLoops)
{
pks <- findpeaks(data[windIdx], threshold = 1.5)
nPeaks[i] <- length(pks)/windowSize
windIdx <- windIdx + windshift
}

My questions now:

1. Is my translation thus far correct?
2. What is the appropriate expression in R for the "nan" in MATLAB at the beginning of the loop when nPeaks is defined? I simply went with "matrix" now but is this doing the same thing?
3. How do you translate this MATLAB term from the loop "pks = findpeaks(data(windIdx))" into R. I think I understand what the line is suppose to do: "calculate the peaks in a defined window which through each iteration loop will shift"... but how can you do this kind of referencing in R? I tried to simply use [] instead of () but did not help at all.

I would be so grateful if anyone could help me with this!!!
Best wishes
Stef

[AlexisW]

You need to go step-by-step and check what each function is doing. For example for your question 2, you can test with nLoops = 3:

matrix(3, 1)
#>      [,1]
#> [1,]    3

This is probably not what you want. You can run ?matrix to see the documentation, in that case you might mean:

matrix(NA, nrow = 3, ncol = 1)
#>      [,1]
#> [1,]   NA
#> [2,]   NA
#> [3,]   NA

Note however that in MATLAB, everything is a matrix, but in R, everything is a vector. So since your data seems to be mostly 1-dimensional, I don't think you need any matrix anywhere.

I think that will also solve your question 3: [a] is appropriate to subset a vector, whereas if you're trying to subset a matrix you would need [a,b] (specifying both the rows and columns).

Also make sure that you visually check your results at each step, you can use plot(x,y) and abline(v=x) for that.

Feel free to ask pointed questions about a specific step that is not giving the result you want: you will get much better (more) answers with precise questions than just asking us to read and understand all of your untested code.