How to apply FFT or LombScargle to a discontinous signal?
2 views (last 30 days)
Show older comments
Hi! I have a code that gives me a signal like the one in the figure where I have some empty gaps because the signal is not continous. I need to apply the FFT or Lomb Scargle or any other method in order to find the predominant frequencies of my motion. I have applied basic MATLAB commands for the fft but I get an empty plot. Can anyone help me? Thanks in advance
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
plot(elapsed_seconds(steps),m_app(steps));
axis([0 600 2 10.5])
xlabel('Times(s)')
ylabel('Apparent Magnitude')
title('Light Curve')
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
0 Comments
Accepted Answer
Mathieu NOE
on 24 Nov 2021
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
8 Comments
More Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!