Code covered by the BSD License

### Highlights from Lomb (Lomb-Scargle) Periodogram

4.66667
4.7 | 6 ratings Rate this file 68 Downloads (last 30 days) File Size: 2.11 KB File ID: #20004

# Lomb (Lomb-Scargle) Periodogram

### Dmitry Savransky (view profile)

21 May 2008 (Updated )

Computes the Lomb normalized periodogram (spectral power as a function of frequency).

File Information
Description

LOMB(T,H,OFAC,HIFAC) computes the Lomb normalized periodogram (spectral
power as a function of frequency) of a sequence of N data points H,
sampled at times T, which are not necessarily evenly spaced. T and H must
be vectors of equal size. The routine will calculate the spectral power
for an increasing sequence of frequencies (in reciprocal units of the
time array T) up to HIFAC times the average Nyquist frequency, with an
oversampling factor of OFAC (typically >= 4).

The returned values are arrays of frequencies considered (f), the
associated spectral power (P) and estimated significance of the power
values (prob). Note: the significance returned is the false alarm
probability of the null hypothesis, i.e. that the data is composed of
independent gaussian random variables. Low probability values indicate a
high degree of significance in the associated periodic signal.

Although this implementation is based on that described in Press,
Teukolsky, et al. Numerical Recipes In C, section 13.8, rather than using
trigonometric rercurrences, this takes advantage of MATALB's array
operators to calculate the exact spectral power as defined in equation
13.8.4 on page 577. This may cause memory issues for large data sets and
frequency ranges.

Example
[f,P,prob] = lomb(t,h,4,1);
plot(f,P)
[Pmax,jmax] = max(P)
disp(['Most significant period is',num2str(1/f(jmax)),' with FAP of ',num2str(prob(jmax))])

MATLAB release MATLAB 7.0.1 (R14SP1)
Tags for This File   Please login to tag files.
Comments and Ratings (12)
04 Jul 2014 Hunter

### Hunter (view profile)

Thanks for the script, it really helped a lot. I am a bit confused as what should be the x-axis in the plots: FAP or power spectrum?

23 Apr 2014 John

### John (view profile)

Thank you for creating this script. I have observed significantly different results compared to Matlab's "pwelch" when analyzing the same data (uniformly sampled). Can you please elaborate on the differences between the "lomb" result and Matlab's "pwelch" result? Thank you

To clarify, the differences observed are in spectrum amplitude. The frequencies associated with the peaks are the same.

Comment only
20 Feb 2014 Boris

### Boris (view profile)

To Sophia:
Try to transpose the input arrays.

Comment only
17 Nov 2013 sophia

### sophia (view profile)

In your code, tau, cterm and sterm, i have error massage of " Error using ==> mtimes
Inner matrix dimensions must agree."

how can that possible?

06 Apr 2012 Yoshiki

23 Mar 2010 Lei

### Lei (view profile)

Although the input of the data is the same, the signicant level is different so much when the time resolution is different (for expample, hourly, or daily resolution). I understand the signicant level depends on the number of total frequency. it seems not make sense to me. how to resolve this problem?

Comment only
03 Mar 2010 Sarah Halliday

### Sarah Halliday (view profile)

Thank you very much for this code it has been very useful so far.

I am attempting to analysis 2 1/2 years worth of data collected on a daily timestep. I just wanted to know exactly how I should format the input TIME vector and what the power output units would be??

Additionally, is there any guidance for selecting the best HIFAC and OFAC factors.

Thank you very much, any help would be greatly appreciated.

05 Nov 2009 Filipe

### Filipe (view profile)

I have one question regarding the units that the lomb.m outputs.
The power is in ms^2, s^2, ms^2/Hz, s^2/Hz????

Any help on this would be greatly appreaciated.

Best regards!

Thank you very much in advance!

25 Jun 2009 guj

### guj (view profile)

How you can reconstruct the original data

Comment only
31 Jul 2008 Benj Hodel

Have you thought about coding a faster implementation of this as detailed in W.H. Press and G.B. Rybicki, "Fast Algorithm for Spectral Analysis of Unevenly Sampled Data," Astrophysical Journal, 338, 277 (1989)? I think it would be much quicker and more useful. See http://adsbit.harvard.edu/cgi-bin/nph-iarticle_query?1989ApJ%2E%2E%2E338%2E%2E277P

23 May 2008 Dmitry Savransky

As explained in the description, this implementation calculates the exact normalized periodogram equation in a single step, using all matrix and array operations. lombscargle.m is an almost exact copy of the algorithm from Numerical Recipes in C, which uses trigonometric recurrences and many loops, and runs somewhat slower. Also, this program does not search out significant frequencies or try to reconstruct the waveform.

Comment only
22 May 2008 Scott Miller

How does this compare with lombscargle.m?

Comment only