hello
I used fillmissing to fill the Inf values (created by your log10 of zero in your function)
then the fft will show this spectrum
expanded code :
clc; clear all; close all;
A = load('DatiECICosmoSkymed.txt');
C = load('versori_tumbling.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)];
% phi_time(i_time) = acos(dot(u_sun_time,u_obs_time));
% phi_time(i_time) = convang(phi_time(i_time),'rad','deg');
m_app(i_time) = Brightness_Magnitude(u_sun_time,u_obs_time,u_normale,number_facets,rho,area,height);
end
steps = (1:40:size(elapsed_seconds,1));
figure(1)
plot(elapsed_seconds(steps),m_app(steps));
axis([0 600 2 10.5])
xlabel('Times(s)')
ylabel('Apparent Magnitude')
title('Light Curve')
%% FFT
time = elapsed_seconds(steps);
signal = m_app(steps);
ind = find(signal>1e6); % find Inf values
signal(ind) = NaN; % replace Inf values by NaN
signal2 = fillmissing(signal,'linear'); % replace NaN by linear interpolation
figure(2)
plot(time,signal,'*',time,signal2);
axis([0 600 2 10.5])
xlabel('Times(s)')
ylabel('Apparent Magnitude')
title('Light Curve')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 64; %
OVERLAP = 0.75;
Fs = 1./mean(diff(time));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum] = myfft_peak(signal2,Fs,NFFT,OVERLAP);
figure(3),semilogy(freq,spectrum);
df = freq(2)-freq(1); % frequency resolution
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude');
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
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
