function x = f_alpha_gaussian ( n, q_d, alpha )

%*****************************************************************************80
%
%% F_ALPHA_GAUSSIAN generates noise using a Gaussian distribuion.
%
%  Discussion:
%
%    This function generates a discrete colored noise vector X of size N
%    with a power spectrum distribution that behaves like 1 / f^ALPHA.
%    The white noise used in the generation is sampled from a Gaussian
%    distribution with zero mean and a variance of Q_D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2011
%
%  Author:
%
%    Miroslav Stoyanov
%
%  Reference:
%
%    Jeremy Kasdin,
%    Discrete Simulation of Colored Noise and Stochastic Processes
%    and 1/f^a Power Law Noise Generation,
%    Proceedings of the IEEE,
%    Volume 83, Number 5, 1995, pages 802-827.
%    
%    Miroslav Stoyanov, Max Gunzburger, John Burkardt,
%    Pink Noise, 1/F^Alpha Noise, and Their Effect on Solution
%    of Differential Equations,
%    submitted, International Journal for Uncertainty Quantification.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the discrete noise vector.
%
%    Input, real Q_D, the variance of the Gaussian distribution.  A standard
%    Gaussian distribution has a variance of 1.  The variance is the square
%    of the standard deviation.
%
%    Input, real ALPHA, specifies that the computed noise is to have a 
%    power spectrum that behaves like 1/f^alpha.  Normally 0 <= ALPHA <= 2.
%
%    Output, real X(N), the computed discrete noise vector.
%

%
%  Set the standard deviation of the white noise.
%
  stdev = sqrt ( abs ( q_d ) );
%
%  Generate the coefficients.
%
  hfa = zeros ( 2 * n, 1 );
  hfa(1) = 1.0; 
  for i = 2 : n
    hfa(i) = hfa(i-1) * ( 0.5 * alpha + ( i - 2 ) ) / ( i - 1 );
  end
  hfa(n+1:2*n) = 0.0;
%
%  Sample the white noise.
%
  wfa = stdev * randn ( n, 1 );
%
%  Pad the array with zeros in anticipation of the FFT process.
%
  z = zeros ( n, 1 );
  wfa = [ wfa; z ];
%
%  Perform the discrete Fourier transforms.
%
  fh = fft ( hfa );
  fw = fft ( wfa );
%
%  Multiply the two complex vectors.
%
  fh = fh(1:n+1);
  fw = fw(1:n+1);
      
  fw = fh .* fw;
%
%  Scale to match the conventions of the Numerical Recipes FFT code.
%
  fw(1)   = fw(1)   / 2.0;
  fw(end) = fw(end) / 2.0; 
%
%  Pad the array with zeros in anticipation of the FFT process.
%
  z = zeros ( n - 1, 1 );
  fw = [ fw; z ];
%
%  Take the inverse Fourier transform.
%
  x = ifft ( fw );
%
%  Only the first half of the inverse Fourier transform is useful.
%
  x = 2.0 * real ( x(1:n) );

  return
end