How to modify a for loop script to conduct calculation on a sliding window

I'm working with some tables I've made that are a continuous time-record of hour-by-hour observations at a given site. What I want to do is take this long term record, calculate weekly rates of change (i.e. slope), and do so with a sliding window that conducts this calculation every day. The goal of the latter is to create a continuous record of changes in rates.

The first problem I have solved, admittedly inelegantly, with the following code (I apologize for the sloppiness; I'm new to both MATLAB and programming and any tips on how to streamline this code would be appreciated). DLTT is the name of my rounded and continuous record of observations, stored in a table datatype. I should note also that my data contains NaN's, which is why I filter the data before running polyfit in the script:

r = height(DLTT);
coefficients = 0;
count = 0

for i = 1 : r
    if rem(i,168) == 0
        t = DLTT(i-167:i,:);
        validdata = ~isnan(t.TiltMag);
        t = t(validdata,:);
        x = datenum(t.Date_Time);
        x = x - x(1);
        y = t.TiltMag;
        coefficients = polyfit(x,y,1);
        count = count + 1
        slope(count) = coefficients(1);
        t = []; 
end

I really want to know how to make this script, which does the job of calculating weekly rates wonderfully, able to do the same calculation on a sliding window. Since this is hour by hour data, this would mean after index 168 it should repeat the calculation every 24 rows, but I'm not sure how to tell it to do that.

Any help would be appreciated, thank you.

Answer 1

[Image Analyst]

Why not just use movmean() or conv()?

[Nicholas Barber]

I have looked at both, but I'm not sure in the case of movmean it would get what I want. I want slope on a weekly base, moving day by day. Movmean could give me the mean every 24 hours, but not the weekly rate of change, right? In regards to conv(), I am not exactly sure how it would apply either. If either of those does what I want, then my further point of confusion is where to put that code to make it operate properly.

[Image Analyst]

They average within a sliding window of fixed width. If you need the number of elements in the window to vary (like because your time points are not evenly spaced), then they would not work. If the time between elements is a constant, then a fixed number of elements in the sliding window should work fine, and thus conv() and movmean() would also.

[Nicholas Barber]

Fortunately the time between each row is exactly one hour, so I should be able to use these functions. However, I'm not looking for averages. I'm trying to compute the slope, or the rate of change every week moving from day to day. So starting on Day 7, I want to calculate the rate of change from Day 1 to Day 7. Then on Day 8 I want to calculate the rate of change from Day 2 to Day 8 and so on. Does that make sense? I admit it's kind of a strange request, but for the data I'm working with it would be great to have a continuous estimate of weekly change as we move from day to day.

[Chris Turnes]

To make a slight correction to Image Analyst's answer, you don't need the time points to be evenly spaced if you use the 'SamplePoints' Name-Value pair to specify your timestamps. If you say, for example, movmean(x,hours(24*7),'SamplePoints',t), then the output value y(i) corresponding to timestamp t(i) would be the average of all points x(j) whose corresponding timestamps t(j) are within the 24*7 hour window centered at t(i).

[Chris Turnes]

It sounds like what you are looking for is a linear fit of the data in each window, and you care about the slope parameter. Since your data is evenly spaced at hour intervals, you can consider that for any 24*7 = 168 hour, you can fit a least-squares linear model of your data with:

     a = V\x;

where V(i,:) = [i-84 1];

So what you're trying to do here is essentially compute s = V\x for a sliding window over the x values.

One approach here is to create the V matrix, factor it, and use it to compute the least-squares fit in every window:

     % Take the Q-R decomposition
     [Q,R] = qr([(-84:83)', ones(168,1)], 0);
     % Convolve the columns of Q with the observations.
     y = conv(x, flip(Q(:,1)), 'valid');
     y(:,2) = conv(x, flip(Q(:,2)), 'valid');
     % Compute the slope and intercepts.
     s = y / R;

I've used the 'valid' option here, so you're only getting the windows where there's a full 168 hours worth of data to fit (you'd have to handle the endpoints in a more nuanced way), but you'll see it's equivalent to fitting the slope and intercept for each window:

     >> x = randn(300,1);
     >> % Create the linear fit matrix
     >> V = [(-84:83)', ones(168, 1)];
     >> [Q,R] = qr(V,0);
     >> % Use it in a sliding window method
     >> y = conv(x, flip(Q(:,1)), 'valid');
     >> y(:,2) = conv(x, flip(Q(:,2)), 'valid');
     >> % Use the result to compute the slope and intercept
     >> s = y / R;
     >> s(1,:)
     ans =
         0.0030   -0.0739
     >> % Compare against the first valid window.
     >> V\x(1:168)
     ans =
         0.0030
        -0.0739