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

11 views (last 30 days)
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; ...
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);
end
signal = m_app;
%% LOMB SCARGLE PERIODOGRAM
time = elapsed_seconds;
[pxx,f] = plomb(signal,time);
[pk,f0] = findpeaks(pxx,f,'MinPeakHeight',10);
figure
plot(f,pxx,'k',f0,pk,'o','LineWidth',1.1)
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;
else
F_obs_faccetta(i_face) = 0;
end
end
m_app = -26.7 - 2.5*log10(sum(F_obs_faccetta)/455);
end

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:
plot(time(1:end-1),diff(time),'-r.');
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):
data2save=[time,signal];
save('brightness.txt','data2save','-ascii');
Future scripts can now load the brightness data without having to compute anything.
clear;
data=load('brightness.txt');
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
figure;
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.
  2 Comments
Loren99
Loren99 on 10 Dec 2021
Edited: Loren99 on 10 Dec 2021
@William Rose first of all, thank you for your answer! I have added the missing part of the function that you have supposed to be m_app=sum(facets). However if i try a motion with a very lower frequency (0.032 Hz) with respect to the sampling rate (0.333 Hz), it doesn't seem to me that i get the interference of two sine waves, but only one sinusoidale wave as you can see from the image. In this case, however, I continue to get a signal split by the lomb scargle and I cannot explain it to myself. I attach the lowest frequency file, named spinning02. You can try it in the code in place of spinning06. Thanks in advance!

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.
  2 Comments
William Rose
William Rose on 10 Dec 2021
@Loren99, I get an error when I try to run your code, because function Brightness_Magnitude() does not assign a vaue to the output variable, mapp. There should be a line at the end of Brightness_Magnitude() that looks like
mapp=42; %or
mapp=sum(F_obs_faccetta); %or something else
Please fix it.
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.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!