Using iFFt to convert acceleration PSD ((m/s²)²/Hz) to random acceleration time series (m/s²)
15 views (last 30 days)
Show older comments
Hello everybody,
I am trying to create a random acceleration time series from a given acceleration PSD ((m/s²)²/Hz). I have browsed some of the answers to similar questions in MATLAB community and put together few lines of code.. I do get a random acceleration time series but the acceleration amplitudes are really really small ( I would expect amplitudes at around max 130-140 m/s² and I get amplitudes around max 0.01-0.015 m/s² ). I have used randome phases for the time series. The time series is expected to have a sample rate of 2000Hz.
Please help me to understand where am I going wrong.
I have uploaded 2 MATLAB Files: psd_freq.mat contains the array of frequencies in the PSD (10-1000Hz), psd_values.mat contains the corresponding values of acceleration PSD with the units ((m/s²)²/Hz).
Thank you so much!
% This script generates a time series from a given PSD.
% Random phase angles are used to generate the time signal.
% The given PSD has freq from 10-1000Hz and the time series should also contain
% freq from 10-1000Hz.
% Desired time series duration is 12.8s.
%% Step 1
t_fin = 12.8; % duration of time signal (s)
Fr = 1/t_fin; % freq. resolution of PSD
%% Step 2
F0 = 10; % start freq. of time signal (Hz)
Fmax = 1000; % end freq. of time signal (Hz)
Fs = 2*Fmax; % sampling freq. of time signal
time = [0:1/Fs:t_fin]'; % time array
%% Step 3
% interpolating the given PSD
f = 0:Fr:Fmax;
PSD = [interp1(psd_freq, psd_disp_3150repeats, f)];
PSD(isnan(PSD))=0; % replace the non-existing frequencies in PSD with 0
%% Step 4
% calculating amplitudes from PSD
P1ampl = zeros(size(PSD));
for j=1:length(P1ampl)
if PSD(j)==0
P1ampl(j) = sqrt(Fr*PSD(j));
else
P1ampl(j) = sqrt(2*Fr*PSD(j));
end
end
P1ampl=P1ampl';
%% Step 5
% generating random phases for the time series
P1phase = [rand(1, numel(P1ampl))*(2*pi)]';
%% Step 6
P2 = [P1ampl.*exp(P1phase*1i)]'; % combining amplitudes and phases
P2=[P2,fliplr(conj(P2(2:end)))]; % flipping the array to make a fullsided array
%% Step 7
s=[ifft(P2)]'; % inverse FFT
plot(time,s)
0 Comments
Answers (2)
William Rose
on 27 Jun 2024
Edited: William Rose
on 2 Jul 2024
[edit: I uploaded an updated version of SpectrumAnalysisNotes. The new version has more formulas and examples, n Appendix 11, for creating a random signal with a specified standard deviation and range of frequencies.]
I understand that you wish to create a time domain signal with a desired spectral shape and an expected time domain approximate min and max values.
We cannot run the code above because we do not have "psd_disp_3150repeats".
You say "I would expect amplitudes at around max 130-140 m/s² and I get amplitudes around max 0.01-0.015 m/s²". I am not sure what you mean by "around max". The RMS value of the time domain signal is a more reproducible measure of amplitude or (when squared) power.
I cannot tell why you expect that the code above will generate a time domain signal whose "around max" values are 130-140 m/s^2.
You define F0=10 in your script, but you never use F0 after defining it.
You say you want the signal to have power from 10-1000 Hz. Your sampling rate is 2000 Hz. I strongly recommend a higher sampling rate. I realize the Nyquist criterion says that you only need sample at twice the highest frequency in the signal. But this means you are representing a 1000 Hz sine wave with just two points. Two points is not enough to give a reliable reconstruction of a sine wave. I would sample at 4 kHz to 5 kHz.
Check out appendices 9 and 10 in the attached document. These appendices discuss the expected power spectrum of a random signal, and different ways of scaling the power spectrum.
More later.
2 Comments
William Rose
on 29 Jun 2024
Edited: William Rose
on 2 Jul 2024
[edit: Replaced SpectrumAnalysisNotes.pdf with new version. Appendix 11 now has more formulas and examples for creating a random signal with specified standard deviation and frequency range.]
Here is the PDF with appendices 9 and 10 - and also appendix 11, which I just added, for you. Appendix 11 describes how to make a random signal with specified spectral properties. Parseval's theorem provides a way to determine the RMS signal amplitude from the power spectrum.
I realize that I have still not explained why the RMS acceleration in your time-domain signal is much smaller than you expect. I suspect it may be due to incorrect conversion from power spectral density to power spectrum to amplitude spectrum. But I am not sure.
Star Strider
on 27 Jun 2024
You cannot invert a power spectral density result, or any result using overlap-add to calculate the spectrum. The reason is simply that too much information is lost, specifically the phase information, and the spectrum resulting from any overlap-add approach is no longer the sort of one-to-one mapping tthat a Fourier transform provides, and that the inversion requires. Additionally, the power spectral density units are in power_unit/frequency_unit, and it is simply not possible to back-calculate that to a power spectrum, much less to a simple Fourier transform.
The Fourier transform and the power spectrum or power spectral density calculations and results should be kept separate.
0 Comments
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!