File Exchange

image thumbnail

Lomb normalized periodogram

version 1.2.0.0 (6.57 KB) by Christos Saragiotis
Both functions calculate the Lomb-Scargle periodogram (aka Gauss-Vanicek/Least-squares spectrum)

1 Download

Updated 04 Dec 2008

View Version History

View License

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)

Cite As

Christos Saragiotis (2021). Lomb normalized periodogram (https://www.mathworks.com/matlabcentral/fileexchange/22215-lomb-normalized-periodogram), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (12)

Stefan Spulber

Thanks a lot for sharing the functions, this implementation is so fast I can now plot spectrograms in real-time!

Uli Wortmann

David Verrelli

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).

David Verrelli

ACCURACY COMMENTS
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.

David Verrelli

SPEED COMMENTS
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

D

Edoardo Alberti

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

Felipe G. Nievinski

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')

Christos Saragiotis

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

Benjamin Bayes

Benjamin Bayes

Excellent function. However, in fastlomb.m, it seems like on line 212,

ck = 1 + mod(fix(t(j)*fac),2*nfreq);

should be

ck = 1 + mod(t(j)*fac,2*nfreq);

If not changed, it seems like the else statement on line 228 is never reached. Thanks!

Daniel Armyr

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.

MATLAB Release Compatibility
Created with R2008a
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired by: lombscargle.m

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!