Inverting Logic in an Algorithm

作者 Loren Shure, January 18, 2006

1 次查看（过去 30 天） | 0 个赞 | 7 个评论

A while ago, a colleague needed to use the sinc function to interpolate some channel data posted by the IEEE 802.15.3a working group. The files they supplied were large and the interpolation took a while. Seems like a job worthy inspecting the code and using the profiler (see references below). Here's what you'll find.

The sinc function used to contain this code, avoiding directly evaluating sinc(0):

y=ones(size(x));
i = find(x);
y(i)=sin(pi*x(i))./(pi*x(i));

The output is preloaded with ones and the function is only evaluated for nonzero input values. The problem with this is that the indexing in the third line is very expensive in terms of time.

Thinking about using the sinc function, it often gets used for signal processing waveform creation. So ask yourself, in a signal, how many of the time elements are 0 – and the answer is typically a very small number, usually either 0 or 1.

The following approach seems to run about 30 to 40% faster, based on a very unscientific study, and is what ships today in the Signal Processing Toolbox.

i=find(x==0);
x(i)= 1;  % From LS: don't need this if /0 warning is off
y = sin(pi*x)./(pi*x);
y(i) = 1;

In this case, we look for the zero values (most likely a small number of them), replace them with ones, and then evaluate the expression on the whole input array. After that's complete, we replace the output values corresponding to input values of zero with ones.

Note the comment on the second line. If warnings for zero-division are off, then removing that line gives even more speed improvement. Why? Because a copy of the x array needs to be made if the second line is there since some values in x are being modified; the local copy no longer matches the array passed into the function. MATLAB will only use extra memory for an input array if it gets modified inside the function.

Finally, an even shorter version, where we don't worry about the 0 division warning might look like this:

y = sin(pi*x)./(pi*x);
y(x==0) = 1;

In this version, we are performing the computation on the original array x, and then operating on the output, y using logical indexing to replace values where x = 0 with the value 1, taking advantage of MATLAB's scalar expansion.