> > 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.
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 time-series displacement/velocity data from a PSD % using the method described here:
% 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)
%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:N-1) 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)
%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).