matlabsignal-processingspectral-clustering

implementation of Lomb-Scargle periodogram

from matlab official site , Lomb-Scargle periodogram is defined as

http://www.mathworks.com/help/signal/ref/plomb.html#lomb

suppose we have some random signal let say

 x=rand(1,1000);

average of this signal can be easily implemented as

average=mean(x); variance can be implemented as

>> average_vector=repmat(average,1,1000);
>> diff=x-average_vector;
>> variance= sum(diff.*diff)/(length(x)-1);

how should i continue? i mean how to choose frequencies ? calculation of time offset is not problem,let us suppose that we have time vector

t=0:0.1:99.9;

so that total length of time vector be 1000, in generally for DFT, frequencies bins are represented as a multiplier of 2*pi/N, where N is length of signal, but what about this case? thanks in advance

Solution

  • As can be seen from the provided link to MATLAB documentation, the algorithm does not depend on a specific sampling times tk selection. Note that for equally spaced sampling times (as you have selected), the same link indicates:

    The offset depends only on the measurement times and vanishes when the times are equally spaced.

    So, as you put it "calculation of time offset is not a problem".

    Similar to the DFT which can be obtained from the Discrete-Time Fourier Transform (DTFT) by selecting a discrete set of frequencies, we can also choose f[n] = n * sampling_rate/N (where sampling_rate = 10 for your selection of tk). If we disregard the value of PLS(f[n]) for n=mN where m is any integer (since it's value is ill-formed, at least in the formula posted in the link), then:

    \sum_{k=1}^N \cos^2(2\pi f_n t_k) = \sum_{k=1}^N \sin^2(2\pi f_n t_k) = N/

    Thus for real-valued data samples:

    P_{LS}(f_n) = \frac{1}{N\sigma^2} |Y(n)|^2

    where Y can be expressed in terms of the diff vector you defined as:

    Y = fft(diff);

    That said, as indicated on wikipedia, the Lomb–Scargle method is primarilly intended for use with unequally spaced data.