title: Synthesize PSD in MatLab
body: <p>I want to recover a time signals from a given power spectral density, assuming a normal distribution of the original signal:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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)].</p>