Both functions caculate the Lomb normalized periodogram (aka Lomb-Scargle, Gauss-Vanicek or Least-Squares spectrum) of a vector x with coordinates in t, which is essentially a generalization of the DFT for unevenly sampled data.
The codes are transcriptions from Fortran of the subroutines found in Section 13.8 (pp. 569-577) of "Numerical recipes in Fortran 77: the art of scientific computing", 2nd ed., vol. 1, Cambridge University Press, NY, USA, 2001 by WH Press, SA Teukolsky, WT Vetterling and BP Flannery,
However, Matlab's characteristics have been taken into account in order to make it fast for Matlab.
FASTLOMB is much faster than LOMB (especially when the length of the input increases) but even LOMB is faster than any other implementation I found in FileExchange. Also they both do not suffer from memory problems (I tested them both for inputs of 100,000 samples).
I'd also like to acknowledge file ID: 20004 (for some reason I can't get two file IDs in the acknowledgements)
Christos Saragiotis (2021). Lomb normalized periodogram (https://www.mathworks.com/matlabcentral/fileexchange/22215-lomb-normalized-periodogram), MATLAB Central File Exchange. Retrieved .
Thanks a lot for sharing the functions, this implementation is so fast I can now plot spectrograms in real-time!
Erratum: my previous post indicated the default MACC=4. That was the value used by Press et al. However, Saragiotis used MACC=10 in his code (which I then used as the baseline).
There are potential problems with the technique used to eliminate duplicates as consideration has not been given that the duplicate time values might be associated with different values of "x".
Hence comparisons of accuracy or timing between these m-files and alternative routines needs to be careful of different processing of duplicate input data points.
On another matter, it would also be advantageous to open the parameter MACC to the user to choose a value. As it stands, errors from the fast algorithm could be of order 2% (with limited testing) using the default MACC=4 [following Press et al.]. However, the discrepancy can be reduced by a factor of circa 100 (to an error of order 0.02%) by increasing MACC to 100; while this does reduce the speed of the "fast" algorithm, it remained quicker than the "slow" algorithm (in limited testing), and therefore also remained quicker than the two older routines from Shoelson and Savransky respectively.
For me, on an old PC (Intel Core2 Duo E8500 @ 3.16/3.17 GHz, 8 GB RAM, Win7 64bit) both of these algorithms are so far (with limited testing) significantly faster than the alternatives of Brett Shoelson (File ID: #993) and Dmitry Savransky (File ID: #20004).
FOR EXAMPLE, with a data set of 4403 data points (with no time duplicates):
Shoelson ~ 69 s
Shoelson (with processing for duplicates) ~ 84 s
Savransky ~ 93 s
Saragiotis (slow) ~ 10.5 s
Saragiotis (fast) ~ 0.2 s
Saragiotis (fast, with MACC=100) ~ 2.3 s
hello Christos (or anybody), can you specify what are the output units? is it a power spectrum?
or a power spectrum density? thanks in advance, Edo
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
if length(x)~=nt, warning('fastlomb:Duplicates','Double entries have been eliminated.'); end
so that it can be disabled issuing warning('off','fastlomb:Duplicates')
Daniel, thank you for your commenting.
Regarding the plotting part, I am not sure I understand what you mean, but you can always set fig to 0 and then take the output of the m-files and plot them as you wish.
Sorry it took me so long to answer, I thought I was getting an automatic email every time someone commented but it seems this is not the case
Excellent function. However, in fastlomb.m, it seems like on line 212,
ck = 1 + mod(fix(t(j)*fac),2*nfreq);
ck = 1 + mod(t(j)*fac,2*nfreq);
If not changed, it seems like the else statement on line 228 is never reached. Thanks!
Nice implelemtation of the Lomb method. The only reason I am not giving this a five is that the built-in plotting functionality is far to over-specified. Opening up a new figure and then drawing raw lines in it is not particularly flexible or helpful for integrating into larger applications.
Inspired by: lombscargle.m
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!