Using iFFt to convert acceleration PSD ((m/s²)²/Hz) to random acceleration time series (m/s²)

45 次查看（过去 30 天）
Tiasa Ghosh
Tiasa Ghosh 2024-6-27
编辑： William Rose 2024-7-2
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)

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

回答（2 个）

William Rose
William Rose 2024-6-27
编辑：William Rose 2024-7-2
@Tiasa Ghosh,
[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 个评论
Tiasa Ghosh
Tiasa Ghosh 2024-6-27
Thank you for your prompt response. @William Rose
The 'psd_disp_3150repeats' is available as .mat variable in the attached MATLAB file 'psd_values.mat'. To avoid confusion, I have attached the variable as .mat again to this reply.
I expect the 'max values around 130-140m/s²' because I am comparing the results from another software where I input the same PSD to get a time series and it gives a series with those appx values.
I agree, RMS would be a more appropriate measure. In this case, I compared the RMS of the resulting time series from MATLAB with the reference. The RMS of the time signal from the reference is 36 while that of the signal from my code is 0.0029.
I would still like to keep the sampling rate of the time series as 2000Hz and rather use the power from the PSD till 500Hz.
The attached document 'spectrumAnalysisNotes' include Apendices only till 8. Could you please resend the document with two further Appendices?
Thank you again.
William Rose
William Rose 2024-6-29
编辑：William Rose 2024-7-2
@Tiasa Ghosh,
[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
Star Strider 2024-6-27
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.

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by