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.
