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

42 次查看（过去 30 天）
종영
종영 2024-10-28，10:06
回答： Vilnis Liepins 2024-10-28，14:15
Ran in:
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;

请先登录，再回答此问题。

回答（2 个）

Star Strider
Star Strider 2024-10-28，11:18
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.
Vilnis Liepins
Vilnis Liepins 2024-10-28，14:15
Ran in:
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.

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by