plomb ifft sequence??? please help me!!!!

Plomb was conducted, but modification was not performed to fill the gap here.

After that, if you do an ift, the position of the gap and the position of the peak must match, but as shown in the last picture, the position is misaligned.

Where is the problem?

clc;
clear;
close all;

sampling = 1000;

t = linspace(0, 50, sampling);  
a = [8, 10, 7, 6, 11];
b = [1, 2, 3, 4, 5];   
omega = 2 * pi * [0.2, 0.4, 0.6, 0.8, 1.0]; 
Theta = [pi/6, pi/4, pi/3, pi/2, pi/8];  

y = zeros(size(t));
for k = 1:5
    A_k = sqrt(a(k)^2 + b(k)^2); 
    phi_k = atan(b(k) / a(k));    
    y = y + A_k * cos(omega(k) * t - Theta(k) - phi_k);  
end     

noise = 10 * randn(size(y)); 
y_noisy = y + noise;      
y_gap = y_noisy;

start_idx_1 = round(0.1 * sampling);
end_idx_1 = round(0.2 * sampling);  
start_idx_2 = round(0.4 * sampling);
end_idx_2 = round(0.55 * sampling);  

y_gap(start_idx_1:end_idx_1) = NaN;  
y_gap(start_idx_2:end_idx_2) = NaN;

y_zero_filled = y_gap;
y_zero_filled(isnan(y_zero_filled)) = 0;

t_valid = t(~isnan(y_gap)); 
y_valid = y_gap(~isnan(y_gap)); 
y_interpolated = interp1(t_valid, y_valid, t, 'linear'); 

% 원래 신호 플롯
figure(1);
plot(t, y_interpolated, 'k');
title('a) Artificial series of 5 superposed sine waves');
xlabel('t [day]');
ylabel('y');
grid on;

%% Lomb-Scargle 분석을 위한 변수 준비
N = length(y_interpolated);
tave = (t(1) + t(end)) / 2;

sigma = var(y_interpolated, 1);
f_s = 1 / (t(2) - t(1)); 
f = (0:N-1) * (f_s / N); 

omega = 2 * pi * f; 
theta = 0.5 * atan(sum(sin(2 * omega .* t)) / sum(cos(2 * omega .* t)));


R = sum(y_interpolated .* cos(omega.' * t - theta), 2).';
I = sum(y_interpolated .* sin(omega.' * t - theta), 2).';
C = sum(cos(omega.' * t - theta).^2, 2).';
S = sum(sin(omega.' * t - theta).^2, 2).';

P = 1 / (2 * sigma) * (R.^2 ./ C + I.^2 ./ S);

figure(2);
plot(f, P);

threshold = mean(P)*20;
P_modu = P;
% P_modu(P_modu<threshold) = 0;

figure(3)
plot(f,P_modu)

A_FT = sqrt(N / 2 * sigma * P_modu);
phase = atan2(I, R) + omega * tave + theta;

R_FT = A_FT .* cos(phase); 
I_FT = A_FT .* sin(phase);  


F_positive = R_FT + 1i * I_FT;
F_negative = flip(conj(F_positive));  


F = [complex(0, 0), F_positive F_negative];  
% F = F(1:length(F)/2);

F_if = ifft(F, 'symmetric');  
t_reconstructed = linspace(0, 50, length(F_if));

figure;
plot(t, y_interpolated, 'k', 'DisplayName', 'Interpolated Signal'); 
hold on;
plot(t_reconstructed, real(F_if), 'r--', 'DisplayName', 'Reconstructed Signal'); 
title('Original and Reconstructed Signals');
xlabel('Time (day)');
ylabel('Amplitude');
legend;
grid on;

Answer 1

[Star Strider]

The plomb function returns the power spectral density of the input time-domaiin signal. In order to do in inversee Fourier transform, you need to have the phase, and the PSD proceess destroys that information. To the best of my knowledge, it is not possible to reconstruct the phase information with any accuracy. Calculating the inverse Fourier transform of the PSD will not reconstruct the orignal signal.

I am not certain what you want to do, however consider using the nufft function (introduced in R2020a) instead.

Answer 2

[Vilnis Liepins]

I suggest to try Extended DFT available on fileexchange https://www.mathworks.com/matlabcentral/fileexchange/11020-extended-dft

Added a few lines to your code:

clc;
clear;
close all;
rng(1,"twister"); % to generate the same random sequence
sampling = 1000;
t = linspace(0, 50, sampling);  
a = [8, 10, 7, 6, 11];
b = [1, 2, 3, 4, 5];   
omega = 2 * pi * [0.2, 0.4, 0.6, 0.8, 1.0]; 
Theta = [pi/6, pi/4, pi/3, pi/2, pi/8];  
y = zeros(size(t));
for k = 1:5
    A_k = sqrt(a(k)^2 + b(k)^2); 
    phi_k = atan(b(k) / a(k));    
    y = y + A_k * cos(omega(k) * t - Theta(k) - phi_k);  
end     
noise = 10 * randn(size(y)); 
y_noisy = y + noise;      
y_gap = y_noisy;
start_idx_1 = round(0.1 * sampling);
end_idx_1 = round(0.2 * sampling);  
start_idx_2 = round(0.4 * sampling);
end_idx_2 = round(0.55 * sampling);  
y_gap(start_idx_1:end_idx_1) = NaN;  
y_gap(start_idx_2:end_idx_2) = NaN;
y_zero_filled = y_gap;
y_zero_filled(isnan(y_zero_filled)) = 0;
t_valid = t(~isnan(y_gap)); 
y_valid = y_gap(~isnan(y_gap)); 
y_interpolated = interp1(t_valid, y_valid, t, 'linear'); 
% 원래 신호 플롯
figure(1);
plot(t, y_interpolated, 'k');
title('a) Artificial series of 5 superposed sine waves');
xlabel('t [day]');
ylabel('y');
grid on;

N = length(y_interpolated);
f_s = 1 / (t(2) - t(1)); 

% my code:
F=edft(y_gap); % Function edft.m is available on fileexchange. 

Unrecognized function or variable 'edft'.

Y=real(ifft(F));
hold on
plot(t,Y,'b');
hold off

As you can be seen in the figure, the IFFT applied to Extended DFT output returns Y where gaps (NaNs) in y_gap are filled with interpolated values while available samples are not changed.