Christos Saragiotis

Dear William,

1. Regarding the "slowness issue":

This function is in fact FASTER than the submissions by Telasula and Barrows (btw Barrowes' submission is not something you can run, unless you correct it).

For example for an input of 200,000 samples, in my machine I get:
      Barrows: 5.5 sec, Telasula: 3.1 sec, this function: 1.8 sec

This function is indeed slower than Barrowes' and Talasula's ones when you have inputs with very few samples (say 20 samples or so) but in that case the time of execution is of the order of msec. I think that 3 msec vs. 0.3 msec is not a waste of anyone's time.

2. Regarding the "accuracy issue":

Abramowitz's and Stegun's tables are accurate to the 10th decimal. Interpolation for other values was not my idea it was Abramowitz's and Stegun's idea. The accuracy of the interpolation was commented in the description of the file above. I don't see how your "findings" on the accuracy contribute anything new. By the way, what is the accuracy of Barrowes' or Telasula's submissions?

In my opinion it is very unfortunate that you gave a bad rating based
 - on something which doesn't hold ("slowness") and
 - on a limitation ("accuracy") which is however stated in both the description and the help of the function.

Very good and an excellent example of terse programming.

I have two remarks though:
1. when selecting the 'both' option some peaks or troughs may be lost. E.g consider the vector [0, 5 -5, 5, 0]. Clearly -5 is a trough but localpeaks will not pick it.
2. if we consider again the vector [ 0 5 5 5 0], localpeaks will pick the first and last 5 as a peak but not the second one. I'm not sure if this was intentional but the way I see it it should either pick all of them or none.

Does the job.
You might want to replace line 198:

    if length(x)~=nt, disp(sprintf('WARNING %s: Double entries have been eliminated',mfilename)); end

for:

    if length(x)~=nt, warning('fastlomb:Duplicates','Double entries have been eliminated.'); end

so that it can be disabled issuing warning('off','fastlomb:Duplicates')

This certainly is worth more than only the 1 star claimed by William Sailor, but it can be improved.

The help and error checking are completely satisfactory in my opinion, and an H1 line is provided.

The idea of interpolating the tables can work, but better tables, as well as a better choice of interpolant are necessary. (This code uses pchip, i.e., interp1, with the 'cubic' option, so only a C1 interpolant.) Also, it is important to note that the interpolation error will grow towards the upper end if equally spaced tables are employed.

First, I'll look at the accuracy issue.

n = 10000;
T = sort(rand(n,1)*10);
FresnelCObj = @(t) cos(pi*t.^2/2);

FCquad = zeros(size(T));
for i = 1:n
  FCquad(i) = quadgk(FresnelCObj,0,T(i),'abstol',1e-16);
end
FCpred = fresnel(T,'c',0);

max(FCquad - FCpred)
ans =
      0.000354015139247155

min(FCquad - FCpred)
ans =
     -0.000445656718834175

So in fact, the maximal error of interpolation is significant, and apparently as large as 4e-4 in places. I have verified that quadgk is indeed producing the correct results to within a precision that is close to that requested.

A plot of the interpolation error shows part of the problem, i.e., the use of the interp1 'cubic' option, which is only a C1 interpolant. The errors produced are classically what I would expect to see. One problem is that pchip is simply not designed to interpolate oscillatory curves as we have here.

plot(T,FCquad - FCpred,'.','markersize',4)
grid on
title 'Fresnel prediction errors on test data'

The author can gain an extra digit or so of accuracy in my tests simply by changing to the 'spline' option with interp1. Of course, that would tend to slow down the code slightly too, but that problem too can be resolved with careful coding. Additional accuracy can also be gained by improving the tables used for interpolation, and more accuracy yet obtained by a careful choice of how to build the interpolation tables.

This brings up the speed issue. The fresnel code is indeed reasonably fast, but the issue is whether you call it with one isolated point or a million. Since fresnel is vectorized, it is extremely efficient when called on a large set of points.

For one point on my CPU, I get a time (as predicted by timeit by Steve Eddins) of:

timeit(@() fresnel(pi,'c',0))
ans =
            0.000694679339

In this test, T is a vector with 1e6 elements in it.

timeit(@() fresnel(T,'c',0))
ans =
            0.266534082114

.26653/.000686
ans =
          388.527696793003

Thus the time required per point on the large sample is now only roughly 0.00000027 seconds per point. Well vectorized code pays dividends.

What should I rate this? Since the accuracy can be improved with little effort on the part of the author, I'll give it 4 stars.

This program is a factor of ten slower than the submissions by either Barrowes or Telasula. I am not sure that reading in the tables of Abramowitz and Stegun and interpolating is a good idea, either. The results from this code differs in typically the 5th decimal place vs the other two codes, which agree well with each other.