
Re: random timeseries displacement data from PSD and back again
Posted:
Jun 24, 2013 4:12 PM


> > Jim > > If PSD is the traditional power spectral density, then you have not prepared the Fourier transform correctly before ifft. > To see how to do this correctly, search "phase randomisation" (note spelling) on this newsgroup.
Hi Tideman,
Thanks for the response. I found the examples as suggested, and was able to duplicate fairly well, but still have a couple of questions. First, current code posted below:
%% Description: % This generates random timeseries displacement/velocity data from a PSD % using the method described here:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/42194
% 1. Take the sqrt of your PSD. This gives the amplitude spectrum. % 2. Generate uniform noise for the length of the spectrum using rand % and multiply by 2pi. This is a set of random phases. % 3. Make a set of complex noise using Y=amp.*exp(i*phi) % 4. Replicate complex conj Y=[Y;flipud(conj(Y))]; % 5. Give it zero mean Y=[0;Y] % 6. Inverse transform y=ifft(Y) % 7. Scale to get the desired the std devn
function [t,x,xdot, Fs] = CreateTimeSeries(N, FLAG_plot)
%% Create original PSD
N = 4000 Fs = N F_Max = 2000 gRMS = 2.5
freq = linspace(0,F_Max,N)'; PSD = ones(size(freq)); cumulative = cumtrapz(freq,PSD);
%Scale PSD to desired gRMS PSD = (gRMS^2/max(cumulative)).*PSD;
if(FLAG_plot) figure(11) subplot(2,1,1) semilogx(freq,PSD,'LineWidth',2) hold on end
%% Create ASD ASD = sqrt(PSD);
%% Create random phases rand_phase = 2*pi*rand(size(ASD));
%% Make complex noise Y =ASD.*exp(i*rand_phase);
%% Replicate Complex Conj Y=[Y;flipud(conj(Y))];
%% Give it 0 Mean: Y=[0;Y];
%% IFFT Fs = N dt = 1/N
%Length doubled when we replicated complex conjugate above t = 0:dt:2; y=ifft(Y);
%Parceval's theorem works as follows: % Sum(0:N1) x(i)^2 = (1/N)SUM(F(i)^2) = Average g^2/Hz % SO... sqrt(avg g^2/Hz)*sqrt(Hz) gets us the area = gRMS gRMS_time_series = sqrt(sum(abs(y).^2))*sqrt(F_Max)
if(FLAG_plot) figure(11) subplot(2,1,2) plot(t,y) xlabel('time (s)') ylabel('displacement (m)') title(['gRMS_{time series} = ',num2str(gRMS_time_series)]) end
disp(['Average g^2/Hz = ',num2str(sum(abs(y).^2))])
%% CONVERT BACK TO FREQUENCY DOMAIN FOR GOOD FEELINGS x = y;
N = length(x); xdft = fft(x); xdft = xdft(1:N/2+1); psdx = (1/(Fs*N)).*abs(xdft).^2; psdx(2:end1) = 2*psdx(2:end1); freq_new = 0:(2*F_Max)/length(x):F_Max;
psdx=psdx*(length(psdx))^2;
gRMS_recreate = sqrt(max(cumtrapz(freq_new,psdx)))
if(FLAG_plot) figure(11) subplot(2,1,1) semilogx(freq_new,psdx,'r'); grid on; title(['gRMS_{original} = ',num2str(gRMS_original),', gRMS_{recreate} = ',num2str(gRMS_recreate)]) legend('PSD_{orig}','PSD_{recreate}') xlabel('Frequency (Hz)'); ylabel('g^2/Hz'); end
%Where I'm confused now is the relative scaling when going between the ifft and fft. When I converted from the time domain back to the frequency domain, I'm just taking half the data, multiplying by 2, and that gets me to the correct g^2/Hz amplitude. You'll note that I multiply and divide the psdx value by N^2, which doesn't make sense. I thought going between we'd always scale by sqrt(N).
Thoughts?

