Generate a time series from power spectral density
34 views (last 30 days)
Show older comments
I have a force PSD of 2X2X63 acting at 2 points of 63 frequencies defined. For each frequency, there is a 2x2 matrix representing the auto and cross variances of the force. I want to construct a time signal of 2xN size where N is the size of the time signal. How do I do it?
4 Comments
Star Strider
on 24 Apr 2024
It is not possible to invert a power spectral density signal. First, the units are in dB/Hz, and second, all the phase information is missing.
Accepted Answer
Star Strider
on 24 Apr 2024
Edited: Star Strider
on 25 Apr 2024
In that event, just create a sine curve of the appropriate frequency or frequencies as determined by the PSD result. This will be a ‘synthesis’ rather than a true inversion, however it is likely as close as you can get with only the PSD and frequency vectors.
Illustrating that —
psd_freq = 0:500;
gp = @(x,d) exp(-(x-d).^2/200);
psd_mag = sum(cell2mat(arrayfun(@(d)gp(psd_freq,d), [10 50 120 250 450], 'Unif',0).')) + 0.01;
figure
plot(psd_freq, psd_mag)
grid
title('Original PSD')
figure
plot(psd_freq, mag2db(psd_mag))
grid
title('Original PSD (dB)')
Fs = max(psd_freq)*2; % Calculated Sampling Frequency (Hz)
L = 5.25; % Length Of Signal (s)
t = linspace(0, Fs*L, Fs*L+1)/Fs; % Time Vector
s = psd_mag(:).*sin(2*pi*t.*psd_freq(:)); % Synthesized Signal Matrix
s = sum(s); % Create One Vector, Add Noise
figure
plot(t, s)
grid
xlim([min(t) max(t)])
xlabel('Time')
ylabel('Amplitude')
figure
pwelch(s,[],[],[],Fs)
Calculatng the PSD of the synthesized signal is not a perfect reproduction of the original PSD, however is reasonably close to it. That implies that the reconstructed time-domain signal is probably representative of the original time-domain signal that created the original PSD. It is not possible to determine how closely the reconstructed time-domain signal reproduces the original time-domain signal without having the original to compare it with.
I encourage you to experiment with this idea.
EDIT — (25 Apr 2024 at 06:10)
Simplified code, expanded discussion.
.
3 Comments
More Answers (0)
See Also
Categories
Find more on Parametric Spectral Estimation 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!