synthesize PSD, scaling and units of ifft(X)

I want to recover a time signals from a given power spectral density, assuming a normal distribution of the original signal:

PSD;                                 % [(m/s)^2/Hz] given spectrum
T = 60;                              % [s] length of original signal
dt = 0.005;                          % [s] time step of original signal
N = T/dt;                            % [-] number of samples
NFFT = 2^nextpow2(N);                % [-] number of bins for FFT
fs = 1/dt;                           % [Hz] sampling frequency
ASD = sqrt(PSD);                     % [(m/s)/sqrt(Hz)] get amplitude spectrum 
omega = 2*pi*rand(NFFT/2,1);         % [rad] generate phase vector
Z = ASD.*exp(1i*omega);              % create complex amplitude vector
Z = [0;Z;flipud(conj(Z))];           % extend to satisfy symmetry
Y = real(ifft(Z));                   % inverse FFT
[PSDY,f] = pwelch(Y,[],[],NFFT,fs);  % generate PSD from Y to compare

The results show a power spectrum several orders of magnitude lower than the original, but the shape matches very good. I guess there is something wrong with the units or there might be a scaling factor missing. I'm not sure about the units of the time signal after ifft, since the amplitude has [(m/s)/sqrt(Hz)].

Answer 1

[Rick Rosson]

   ASD = sqrt(PSD*fs/N);

[Felipe]

There's still a large offset between the spectra. But I found the following scaling that worked:

ASD = sqrt(PSD*fs/2);
...
Y = real(ifft(Z)*sqrt(NFFT+1));

Since I don't fully understand the scaling factors, I'm still interested in any comments.

[Maximilian Weber]

I've been bugged by this for quite some time. Would also appreciate input on this.