Using iFFt to convert acceleration PSD ((m/s²)²/Hz) to random acceleration time series (m/s²)
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)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) 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) ifft, psd, time series, acceleration psd, psd to time series, fft MATLAB Answers — New Questions