Using iFFt to convert acceleration PSD ((m/s²)²/Hz) to random acceleration time series (m/s²)
15 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (2 件)
William Rose
2024 年 6 月 27 日
編集済み: William Rose
2024 年 7 月 2 日
[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 件のコメント
William Rose
2024 年 6 月 29 日
編集済み: William Rose
2024 年 7 月 2 日
[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
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.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Spectral Measurements についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!