MATLAB Central Newsreader - help vectorizing a FOR loop

help vectorizing a FOR loop

GG — Thu, 04 Sep 2008 14:51:01 -0400

Hi Experts,

I'm hoping someone can help me vectorize the following
code. The tricky part is that bin_width is a function of
the loop-index i.

for i = 1:1:10000
smoothline(i) = mean(data(i-bin_width(i)/2:i+bin_width
(i)/2)
end

This code is just a sliding-window average to smooth-out a
noisy data array. The smoothed line is created at each
data point by looking at bin_width(i)/2 datapoints ahead
and bin_width(i)/2 datapoints behind, and taking the mean
of all this data. Then it steps to the next datapoint and
repeats with a new bin_width(i) value.

Thanks in advance, -GG

Re: help vectorizing a FOR loop

Walter Roberson — Thu, 04 Sep 2008 15:05:44 -0400

GG wrote:

> I'm hoping someone can help me vectorize the following
> code. The tricky part is that bin_width is a function of
> the loop-index i.

> for i = 1:1:10000
> smoothline(i) = mean(data(i-bin_width(i)/2:i+bin_width
> (i)/2)
> end

> This code is just a sliding-window average to smooth-out a
> noisy data array. The smoothed line is created at each
> data point by looking at bin_width(i)/2 datapoints ahead
> and bin_width(i)/2 datapoints behind, and taking the mean
> of all this data. Then it steps to the next datapoint and
> repeats with a new bin_width(i) value.

cumsum() the data. Suppose we call that csdata.
Then the mean "around" point i is
(csdata(i+bin_width(i)/2) - csdata(i-bin_width(i)/2)) / (bin_width(i)+1)

Then provided that bin_width is vectorized, the transformed
expression should be vectorizable:

csdata = cumsum(data);
maxi = 10000;
bw = bin_width(1:maxi);
idxhigh = (1:maxi) + bw/2;
idxlow = (1:maxi) - bw/2;
smoothline = (csdata(idxhigh) - csdata(idxlow)) ./ (bw+1);

--
Q = quotation(rand);
if isempty(Q); error('Quotation server filesystem problems')
else sprintf('%s',Q), end

Re: help vectorizing a FOR loop

GG — Fri, 05 Sep 2008 03:49:02 -0400

Hi Walter,

Thanks SOOOOOOO much for your reply. Your code speeds up
my function a factor of 153 times! That's amazing.

Just two small corrections I found in reducing it to
practice that perhaps will save others a little time....

csdata = cumsum(data);
maxi = 10000;
bw = bin_width(1:maxi);
idxhigh = (1:maxi) + bw/2;

% idxlow = (1:maxi) - bw/2;
idxlow = (1:maxi) - bw/2 -1; % this -1 is needed keep ...
% starting datapt of sliding window in the average ...
% calculation

% smoothline = (csdata(idxhigh) - csdata(idxlow)) ./
(bw+1);
smoothline = (csdata(idxhigh) - csdata(idxlow)) ./ (bw'+1);
% I needed to replace bw with bw' to get dimensions right

Thanks a million! -GG

Re: help vectorizing a FOR loop

GG — Fri, 05 Sep 2008 15:18:02 -0400

I have another loop maybe someone can help me with as well. I'm not too optimistic that it can be vectorized but I thought I'd post it as I wouldn't have expected the previous loop above to have been vectorized and I was wrong.

It is a variation of what's above... instead of a centered-sliding window average, this loop is a backward-looking sliding window average. The window size for averaging increases with frequency. The tricky thing is this data has a lot of noise (some noise shoots up quite high, which I call spurs). I identify the spurs (if they're above some threshold level) and don't use them when I calculate the mean_line as I step through the data. The data is in array 'd'.

dx=1:20000;
bin_width = round(0.2*sqrt(dx)) - 2; % makes bin_width
                                      % increase w/freq
bin_width(1:2)=0; % prevents negative bin_width
first_data = 1:1:length(bin_width); % initialize array
first_data = first_data - bin_width +1; % index of 1st
                          % datapt in window for avg'g
mean_line(4) = 1e-5; %give starting pt

for i = 5:1:10000
    threshold(i) = 2 * mean_line(i-1);
    if (d(i) > threshold(i)) % identifies spurs
        spurs = [spurs, i];
    end

% Next, calc the avg by summing all data pts in
% window (including spurs), then subtract the spur
% data, then divide this by the number of non-spur
% data points (which were used to compute this avg).

running_spurs = spurs(find(spurs >= first_data(i)));
% finds spurs in window for avg'g

mean_line(i) = (sum(d(first_data(i):i)) - ...
  sum(d(running_spurs))) /(bin_width(i)- ...
  sum(running_spurs>=first_data(i))); %% computes avg

end

Thanks in advance, -GG