Lomb Scargle Periodogram gives me an unexpected peak in the final plot

3 views (last 30 days)
Loren99 on 9 Dec 2021
Commented: William Rose on 10 Dec 2021
Hi everyone! I have applied plomb to obtain the lomb scargle periodogram of my function. However in addition to the most significant peak, I get another important and lower peak than the first. I can't explain the existence of this second peak and how to remove it. Here my code:
clc; clear all; close all;
A = load('DatiECICosmoSkymednuovaf.txt');
C = load('spinning06.txt');
global Req
Req = 6378;
elapsed_seconds = A(:,1);
num_times = size(elapsed_seconds,1);
m_app = zeros(num_times,1);
%% Geometric characteristics of the spacecraft
height = 600000; % m
number_facets = 6; %cubo + solar arrays
u_normale_body = [1,0,0; ...
-1,0,0; ...
0,1,0; ...
0,-1,0; ...
0,0,1; ...
u_normale = [u_normale_body] ;
rho_body = [0.5,0.5,0.5,0.5,0.5,0.5];
rho = [rho_body] ;
area_facets = [4,4,4,4,4,4]; %m^2
area = [area_facets];
%% Light curve generation and Phase angle
for i_time = 1:num_times
u_sun_time = [C(i_time,2),C(i_time,3),C(i_time,4)];
u_obs_time = [C(i_time,5),C(i_time,6),C(i_time,7)];
m_app(i_time) = Brightness_Magnitude(u_sun_time,u_obs_time,u_normale,number_facets,rho,area,height);
signal = m_app;
time = elapsed_seconds;
[pxx,f] = plomb(signal,time);
[pk,f0] = findpeaks(pxx,f,'MinPeakHeight',10);
xlabel('f (Hz)')
ylabel('Power Spectral Density')
function m_app = Brightness_Magnitude(u_sun,u_obs,u_normale,number_facets,rho,area_facets,height)
F_obs_faccetta = zeros(number_facets,1);
for i_face = 1:number_facets
prodottoscalare1 = dot(u_sun,u_normale(i_face,:));
prodottoscalare2 = dot(u_obs,u_normale(i_face,:));
if prodottoscalare1 > 0 && prodottoscalare2 > 0
F_obs_faccetta(i_face) = (455*rho(i_face)*area_facets(i_face)*prodottoscalare1*prodottoscalare2)/height^2;
F_obs_faccetta(i_face) = 0;
m_app = -26.7 - 2.5*log10(sum(F_obs_faccetta)/455);

Answers (2)

William Rose
William Rose on 10 Dec 2021
Lomb-Scargle is not the way I would estimate the spectrum from this data. Lomb-Scargle is great when there is missing data - which is often the case in spacecraft appplications such as yours, because of observing windows, clouds, etc. But your data is sampled at 3 second intervals, to high precision. You can see that this is true by plotting the sampling interval (i.e. difference between successive times) versus time:
xlabel('Time (s)'); ylabel('\Delta t (s)'); grid on;
Since your data is sampled at a constant rate, you might as well use a standard spectral esitmator such as periodogram() or an autoregressive model (more later).
First, I saved the data in a file, so I would not have to compute it each time, by doing the following after running your script (after I fixed your script by adding m_app=sum(...); to the function):
Future scripts can now load the brightness data without having to compute anything.
time=data(:,1); signal=data(:,2);
We are ready to roll now. It is always a good idea to look at the signal in the time domain, if you have questions about the spectrum. Sometimes something jumps out at you. That is definitely true here.
figure; plot(time,signal,'-r.');
xlabel('Time (s)'); ylabel('Brightness'); grid on;
makes plot below.
This looks like two sine waves of very close frequencies added together , resulting in "beating": they slowly progress between constrictuve and destructiv interference. The math is simple. Look up beat frequency and angle sum formula and amplitude modulaiton. On top of that, there is an interesting sampling thing going on: the sine wave has a period of about 8.8 seconds (which I esitmate by counting 34 peaks in the first 300 seconds). The sampling interval is 3.00 seconds. So you have almost exactly 3 samples per cycle, but not exactly. SO the sampling is slowly changing phase relative to the signal. This is another kind of beating. Technically, it's not aliasing, since the sampling rate (0.333 Hz) is more than twice the main frequency (34/300=0.113 Hz). Nonetheless, the inteaction of the sampling rate with the signal produces weird results. I always tell students to sample at 5 times the highest frequency, or more, to avoid this kind of thing.
I think the split peak in the spectral estimate is an accurate representation of a system with two nearby sinusoids which are creating a beat frequency.
To see if I'm correct, let's make two sine waves with nearby frequencies and add them together. The L-S spectrum shows peaks at about 0.1126 Hx and 0.1148 Hz, so I'll use those frequencies. I'll make a standard sine wave at each frequency, sampled at 2 Hz. THis means about 18 samples per cycle which is enough to be sure I'm not missing any details due to too-low a sampling rate.
t=0:.5:600; %time vector
f1=0.1126; x1=sin(2*pi*f1*t); %first sine wave
f2=0.1148; x2=sin(2*pi*f2*t); %second sine wave
subplot(3,1,1); plot(t,x1,'-k'); %plot first
subplot(3,1,2); plot(t,x2,'-k'); %plot second
subplot(3,1,3); plot(t,x1+x2,'-k'); %plot the sum
This signal has a slow beat frequency, like your signal. Notice the similarity of the time scales of this plot and the time-dmain plot yof your signal.
Now sample the sum of sine waves at 3 second intervals, as in your system:
t3=0:3:600; %time vector: sample every 3 seconds
x1=sin(2*pi*f1*t3); x2=sin(2*pi*f2*t3); %sine waves sampled every 3 sec
figure; plot(t3,x1+x2,'-k'); %plot the sum
This plot has features very similar to your data!
We have made data that looks a lot like yours, by an amazingly simple process: make unit-amplitude sine waves at nearby frequencies, and sample the sum of them at about 1/3 of the average frequency. I didn't even try to adjust the phase angles or the amplitudes of the sine waves. This convinces me that the split peak in the L-S spectral estimate is real: the split peak tells us that the signal is the sum of two sinusoids with almost, but not quite, identical frequencies.

Sign in to comment.

William Rose
William Rose on 9 Dec 2021
I would interpolate your signal at 3 second intervals.* Then I would estimate the spectrum of the interpolated signal with periodogram(), with different window lengths, to see how the spectral estimate differs from the Lomb-Scargle estimate. Or you could use a parametric spectral estimator. For the parametric estimator, you must estimate or guess the model order.
Parametric spectral estimators sometimes yield a split peak where there really is no split. Perhaps the Lomb Scargle estimator does also. In the case of parametric estimators, the cure is to reduce the model order, in many cases. See here for example.
*I choose 3 second sampling because that will produce the same maximm frequency in the spectrum as what you show in your plot.
William Rose
William Rose on 10 Dec 2021
My guess that m_app=sum(facets) was a lucky good guess. The spectrum that results when I put that into the code is very similar to the one you posted, except the magnitude is aprroximately 10^-17. Some of the little bumps are not exactly the same, but it is very close, and it has a prominent split spectral peak.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!