1/f noise generation

Yu Xiaodong wrote:
> Hi All,
> 
> Is there anyone know how to generate 1/f noise for simulating the phase
> noise. Someone told me though passing a white noise to a FIR the noise can
> be generated. Any clue about the FIR?
> 
> Thanks alot
> 
> R. Y.

R.Y.:

Here's some C++ code that I wrote that will create phase noise samples 
from a user-specified phase noise profile. The code uses the GNU 
Scientific Library (freely available).

Good luck,
OUP

//--------------------------------------------------------------------------------------------------
//
//  phase_noise
//
//  Generate random deviates whose spectrum matches a phase noise 
profile. This code uses the
//  GNU Scientific Library (GSL) v1.1.1.
//
//  Revisions
//  ---------
//  25-Jun-02  Initial version written.
//
//--------------------------------------------------------------------------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_vector.h>

#include "phase_noise.h"

//--------------------------------------------------------------------------------------------------
//
//  main
//
//--------------------------------------------------------------------------------------------------

int main (int argc, char *argv[])
{
   int i, numSamples;
   double fs;
   double xi, yi;

   // number of samples
   numSamples = atoi(argv[1]);

   // sample frequency in Hz
   fs = atof(argv[2]);

   InitPhaseProfile();

   double *phi = new double [numSamples];

   GenPhaseNoiseDeviates(phi, numSamples, fs);

   for (i = 0; i < numSamples; i++)
   {
     fprintf(stdout,"%d %19.16e\n", i, phi[i]);
   }

}

//--------------------------------------------------------------------------------------------------
//
//  GenPhaseNoiseDeviates
//
//  Inputs
//  ------
//  numSamples  number of samples to generate
//  fs          sample frequency (a.k.a. PRF)
//
//  Input/Output
//  ------------
//  phiRandom   array to hold random deviates
//
//  Revisions
//  ---------
//  25-Jun-02  Initial version written.
//
//--------------------------------------------------------------------------------------------------
int GenPhaseNoiseDeviates(double phiRandom[], int numSamples, double fs)
{
   int i;

   extern int numPointsProfile;

   //-------------------------------
   // Generate white Gaussian noise
   //-------------------------------

   gsl_rng_env_setup();

   const gsl_rng_type *T = gsl_rng_default;
   gsl_rng *r = gsl_rng_alloc(T);

   double *wn = new double [numSamples];

   for (i = 0; i < numSamples; i++)
   {
     wn[i] = gsl_ran_gaussian(r, 1.0);  // zero mean, unity sigma
     //cout << wn[i] << "\n";
   }

   // Store the white noise samples in a GSL "packed" array. A "packed" 
array holds complex data
   // in a double arrray. The order is real, imag, real, imag, and so on.

   double *packedWn = new double [2*numSamples];

   for (i = 0; i < numSamples; i++)
   {
     REAL_PACKED(packedWn,i) = wn[i];
     IMAG_PACKED(packedWn,i) = 0.0;
   }

   //------------------------------------------------------
   // forward Fourier transform of the white noise samples
   //------------------------------------------------------

   gsl_fft_complex_workspace *workspace = 
gsl_fft_complex_workspace_alloc(numSamples);
   gsl_fft_complex_wavetable *wavetable = 
gsl_fft_complex_wavetable_alloc(numSamples);  // we'll use this for the 
inverse too

   gsl_fft_complex_forward(packedWn, 1, numSamples, wavetable, workspace);

   //--------------------------------------------
   // Construct the aliased phase noise response
   //--------------------------------------------

   int nFold = int( fInputPhaseProfile[numPointsProfile-1] / fs ) + 1;

   if (nFold > nFoldMax) nFold = nFoldMax;

   int nHalf = numSamples/2;

   double *hMag = new double [nHalf];
   double *f    = new double [nHalf];
   double *mSSB = new double [nHalf];

   double fMin = 1.0;
   double fMax = fs / 2.0;
   double fStep = (fMax - fMin)/(nHalf-1);

   int k;
   double summation;

   for (i = 0; i < nHalf; i++)
   {
     summation = 0.0;
     f[i] = fMin + i*fStep;
     for (k = -nFold; k <= nFold; k++)
     {
       summation = summation + GetPhaseNoiseValue(f[i] + k*fs);
     }
     hMag[i] = summation;
     mSSB[i] = GetPhaseNoiseValue(f[i]);
   }

//   for (i = 0; i < nHalf; i++)
//   {
//     fprintf(stdout,"%f %19.16e %19.16e\n", f[i], mSSB[i], hMag[i]);
//   }

   // Construct a double-sided phase noise profile. This has even 
symmetry. Note that the
   // GSL stores the positive frequencies first, then the negative 
frequencies in reverse
   // order. We'll construct our DSB profile the same way.

   double *mDSB = new double [numSamples];

   for (i = 0; i < nHalf; i++)
   {
     mDSB[i]       = hMag[i];
     mDSB[i+nHalf] = hMag[nHalf - 1 - i];
   }

   double *packedDSB = new double [2*numSamples];

   for (i = 0; i < numSamples; i++)
   {
     REAL_PACKED(packedDSB,i) = mDSB[i];
     IMAG_PACKED(packedDSB,i) = 0.0;
     //cout << mDSB[i] << "\n";
   }

   //----------------------------
   // Convolve the two sequences
   //----------------------------

   double *packedY = new double [2*numSamples];

   gsl_complex h, x, z;
//   double hx, hy, wx, wy;

   for (i = 0; i < numSamples; i++)
   {
//     wx = REAL_PACKED(packedWn,i);
//     wy = IMAG_PACKED(packedWn,i);

//     hx = REAL_PACKED(packedDSB,i);
//     hy = IMAG_PACKED(packedDSB,i);

//     GSL_SET_COMPLEX(&x, wx, wy);
//     GSL_SET_COMPLEX(&h, hx, hy);

     GSL_SET_COMPLEX(&x, REAL_PACKED(packedWn,i), IMAG_PACKED(packedWn,i));
     GSL_SET_COMPLEX(&h, REAL_PACKED(packedDSB,i), 
IMAG_PACKED(packedDSB,i));

     z = gsl_complex_mul(h,x);

     REAL_PACKED(packedY,i) = GSL_REAL(z);
     IMAG_PACKED(packedY,i) = GSL_IMAG(z);

//     packedY[i] = packedDSB[i] * packedWn[i];  // this only works 
because imaginary part is zero!
   }

   //---------------------------
   // Inverse Fourier transform
   //---------------------------

   gsl_fft_complex_inverse(packedY, 1, numSamples, wavetable, workspace);

   for (i = 0; i < numSamples; i++)
   {
     phiRandom[i] = REAL_PACKED(packedY,i);
   }

   // free some stuff

   gsl_fft_complex_workspace_free(workspace);
   gsl_fft_complex_wavetable_free(wavetable);

   gsl_spline_free (spline);
   gsl_interp_accel_free(acc);

   return 0;
}

//--------------------------------------------------------------------------------------------------
//
//  InitPhaseProfile
//
//  Initialize stuff for the phase profile interpolation.
//
//  Revisions
//  ---------
//  25-Jun-02  Initial version written.
//
//--------------------------------------------------------------------------------------------------
int InitPhaseProfile()
{
   int i;
   extern int numPointsProfile;

   numPointsProfile = sizeof(fInputPhaseProfile)/sizeof(double);

   double *fLog = new double [numPointsProfile];

   extern double fInputPhaseProfile[], mInputPhaseProfile[];
   extern gsl_interp_accel *acc;
   extern gsl_spline *spline;

   for (i = 0; i < numPointsProfile; i++)
   {
     fLog[i] = log10(fInputPhaseProfile[i]);
   }

   //   printf ("#m=0,S=2\n");
   //   for (i = 0; i < numPointsProfile; i++)
   //   {
   //     printf ("%g %g\n", fInputPhaseProfile[i], mInputPhaseProfile[i]);
   //   }

   acc    = gsl_interp_accel_alloc ();
   spline = gsl_spline_alloc (gsl_interp_linear, numPointsProfile);

   gsl_spline_init (spline, fLog, mInputPhaseProfile, numPointsProfile);

   return 0;
}

//--------------------------------------------------------------------------------------------------
//
//  GetPhaseNoiseValue
//
//  Given the desired frequency, this function uses linear interpolation 
in log-space to
//  calculate the phase noise power.
//
//  Input
//  -----
//  f      frequency (Hz)
//
//  Output
//  ------
//  sf     phase noise power in absolute units
//
//  Revisions
//  ---------
//  25-Jun-02  Initial version written.
//
//--------------------------------------------------------------------------------------------------
double GetPhaseNoiseValue(double f)
{
   double sf;
   const double noiseFloor = -150.0;

   extern int numPointsProfile;
   extern double fInputPhaseProfile[];
   extern gsl_interp_accel *acc;
   extern gsl_spline *spline;

   f = abs(f);

   if (f > fInputPhaseProfile[numPointsProfile-1])
   {
     sf = noiseFloor;   // dBc/Hz
   }
   else
   {
     // interpolate to get phase noise power at desired frequency
     sf = gsl_spline_eval (spline, log10(f), acc);
   }

   sf = pow(10.0, 0.1*sf);

   return(sf);
}

>Yu Xiaodong wrote:
>> Hi All,
>> 
>> Is there anyone know how to generate 1/f noise for simulating the
phase
>> noise. Someone told me though passing a white noise to a FIR the noise
can
>> be generated. Any clue about the FIR?
>> 
>> Thanks alot
>> 
>> R. Y.
>
>R.Y.:
>
>Here's some C++ code that I wrote that will create phase noise samples 
>from a user-specified phase noise profile. The code uses the GNU 
>Scientific Library (freely available).
>
>Good luck,
>OUP
>
>//--------------------------------------------------------------------------------------------------
>//
>//  phase_noise
>//
>//  Generate random deviates whose spectrum matches a phase noise 
>profile. This code uses the
>//  GNU Scientific Library (GSL) v1.1.1.
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>
>#include <stdio.h>
>#include <stdlib.h>
>#include <iostream.h>
>#include <math.h>
>
>#include <gsl/gsl_complex.h>
>#include <gsl/gsl_complex_math.h>
>#include <gsl/gsl_errno.h>
>#include <gsl/gsl_fft_complex.h>
>#include <gsl/gsl_fft_halfcomplex.h>
>#include <gsl/gsl_fft_real.h>
>#include <gsl/gsl_randist.h>
>#include <gsl/gsl_rng.h>
>#include <gsl/gsl_spline.h>
>#include <gsl/gsl_vector.h>
>
>#include "phase_noise.h"
>
>//--------------------------------------------------------------------------------------------------
>//
>//  main
>//
>//--------------------------------------------------------------------------------------------------
>
>int main (int argc, char *argv[])
>{
>   int i, numSamples;
>   double fs;
>   double xi, yi;
>
>   // number of samples
>   numSamples = atoi(argv[1]);
>
>   // sample frequency in Hz
>   fs = atof(argv[2]);
>
>   InitPhaseProfile();
>
>   double *phi = new double [numSamples];
>
>   GenPhaseNoiseDeviates(phi, numSamples, fs);
>
>   for (i = 0; i < numSamples; i++)
>   {
>     fprintf(stdout,"%d %19.16e\n", i, phi[i]);
>   }
>
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  GenPhaseNoiseDeviates
>//
>//  Inputs
>//  ------
>//  numSamples  number of samples to generate
>//  fs          sample frequency (a.k.a. PRF)
>//
>//  Input/Output
>//  ------------
>//  phiRandom   array to hold random deviates
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>int GenPhaseNoiseDeviates(double phiRandom[], int numSamples, double fs)
>{
>   int i;
>
>   extern int numPointsProfile;
>
>   //-------------------------------
>   // Generate white Gaussian noise
>   //-------------------------------
>
>   gsl_rng_env_setup();
>
>   const gsl_rng_type *T = gsl_rng_default;
>   gsl_rng *r = gsl_rng_alloc(T);
>
>   double *wn = new double [numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     wn[i] = gsl_ran_gaussian(r, 1.0);  // zero mean, unity sigma
>     //cout << wn[i] << "\n";
>   }
>
>   // Store the white noise samples in a GSL "packed" array. A "packed" 
>array holds complex data
>   // in a double arrray. The order is real, imag, real, imag, and so
on.
>
>   double *packedWn = new double [2*numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     REAL_PACKED(packedWn,i) = wn[i];
>     IMAG_PACKED(packedWn,i) = 0.0;
>   }
>
>   //------------------------------------------------------
>   // forward Fourier transform of the white noise samples
>   //------------------------------------------------------
>
>   gsl_fft_complex_workspace *workspace = 
>gsl_fft_complex_workspace_alloc(numSamples);
>   gsl_fft_complex_wavetable *wavetable = 
>gsl_fft_complex_wavetable_alloc(numSamples);  // we'll use this for the 
>inverse too
>
>   gsl_fft_complex_forward(packedWn, 1, numSamples, wavetable,
workspace);
>
>   //--------------------------------------------
>   // Construct the aliased phase noise response
>   //--------------------------------------------
>
>   int nFold = int( fInputPhaseProfile[numPointsProfile-1] / fs ) + 1;
>
>   if (nFold > nFoldMax) nFold = nFoldMax;
>
>   int nHalf = numSamples/2;
>
>   double *hMag = new double [nHalf];
>   double *f    = new double [nHalf];
>   double *mSSB = new double [nHalf];
>
>   double fMin = 1.0;
>   double fMax = fs / 2.0;
>   double fStep = (fMax - fMin)/(nHalf-1);
>
>   int k;
>   double summation;
>
>   for (i = 0; i < nHalf; i++)
>   {
>     summation = 0.0;
>     f[i] = fMin + i*fStep;
>     for (k = -nFold; k <= nFold; k++)
>     {
>       summation = summation + GetPhaseNoiseValue(f[i] + k*fs);
>     }
>     hMag[i] = summation;
>     mSSB[i] = GetPhaseNoiseValue(f[i]);
>   }
>
>//   for (i = 0; i < nHalf; i++)
>//   {
>//     fprintf(stdout,"%f %19.16e %19.16e\n", f[i], mSSB[i], hMag[i]);
>//   }
>
>   // Construct a double-sided phase noise profile. This has even 
>symmetry. Note that the
>   // GSL stores the positive frequencies first, then the negative 
>frequencies in reverse
>   // order. We'll construct our DSB profile the same way.
>
>   double *mDSB = new double [numSamples];
>
>   for (i = 0; i < nHalf; i++)
>   {
>     mDSB[i]       = hMag[i];
>     mDSB[i+nHalf] = hMag[nHalf - 1 - i];
>   }
>
>   double *packedDSB = new double [2*numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     REAL_PACKED(packedDSB,i) = mDSB[i];
>     IMAG_PACKED(packedDSB,i) = 0.0;
>     //cout << mDSB[i] << "\n";
>   }
>
>   //----------------------------
>   // Convolve the two sequences
>   //----------------------------
>
>   double *packedY = new double [2*numSamples];
>
>   gsl_complex h, x, z;
>//   double hx, hy, wx, wy;
>
>   for (i = 0; i < numSamples; i++)
>   {
>//     wx = REAL_PACKED(packedWn,i);
>//     wy = IMAG_PACKED(packedWn,i);
>
>//     hx = REAL_PACKED(packedDSB,i);
>//     hy = IMAG_PACKED(packedDSB,i);
>
>//     GSL_SET_COMPLEX(&x, wx, wy);
>//     GSL_SET_COMPLEX(&h, hx, hy);
>
>     GSL_SET_COMPLEX(&x, REAL_PACKED(packedWn,i),
IMAG_PACKED(packedWn,i));
>     GSL_SET_COMPLEX(&h, REAL_PACKED(packedDSB,i), 
>IMAG_PACKED(packedDSB,i));
>
>     z = gsl_complex_mul(h,x);
>
>     REAL_PACKED(packedY,i) = GSL_REAL(z);
>     IMAG_PACKED(packedY,i) = GSL_IMAG(z);
>
>//     packedY[i] = packedDSB[i] * packedWn[i];  // this only works 
>because imaginary part is zero!
>   }
>
>   //---------------------------
>   // Inverse Fourier transform
>   //---------------------------
>
>   gsl_fft_complex_inverse(packedY, 1, numSamples, wavetable,
workspace);
>
>   for (i = 0; i < numSamples; i++)
>   {
>     phiRandom[i] = REAL_PACKED(packedY,i);
>   }
>
>   // free some stuff
>
>   gsl_fft_complex_workspace_free(workspace);
>   gsl_fft_complex_wavetable_free(wavetable);
>
>   gsl_spline_free (spline);
>   gsl_interp_accel_free(acc);
>
>   return 0;
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  InitPhaseProfile
>//
>//  Initialize stuff for the phase profile interpolation.
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>int InitPhaseProfile()
>{
>   int i;
>   extern int numPointsProfile;
>
>   numPointsProfile = sizeof(fInputPhaseProfile)/sizeof(double);
>
>   double *fLog = new double [numPointsProfile];
>
>   extern double fInputPhaseProfile[], mInputPhaseProfile[];
>   extern gsl_interp_accel *acc;
>   extern gsl_spline *spline;
>
>   for (i = 0; i < numPointsProfile; i++)
>   {
>     fLog[i] = log10(fInputPhaseProfile[i]);
>   }
>
>   //   printf ("#m=0,S=2\n");
>   //   for (i = 0; i < numPointsProfile; i++)
>   //   {
>   //     printf ("%g %g\n", fInputPhaseProfile[i],
mInputPhaseProfile[i]);
>   //   }
>
>   acc    = gsl_interp_accel_alloc ();
>   spline = gsl_spline_alloc (gsl_interp_linear, numPointsProfile);
>
>   gsl_spline_init (spline, fLog, mInputPhaseProfile, numPointsProfile);
>
>   return 0;
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  GetPhaseNoiseValue
>//
>//  Given the desired frequency, this function uses linear interpolation

>in log-space to
>//  calculate the phase noise power.
>//
>//  Input
>//  -----
>//  f      frequency (Hz)
>//
>//  Output
>//  ------
>//  sf     phase noise power in absolute units
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>double GetPhaseNoiseValue(double f)
>{
>   double sf;
>   const double noiseFloor = -150.0;
>
>   extern int numPointsProfile;
>   extern double fInputPhaseProfile[];
>   extern gsl_interp_accel *acc;
>   extern gsl_spline *spline;
>
>   f = abs(f);
>
>   if (f > fInputPhaseProfile[numPointsProfile-1])
>   {
>     sf = noiseFloor;   // dBc/Hz
>   }
>   else
>   {
>     // interpolate to get phase noise power at desired frequency
>     sf = gsl_spline_eval (spline, log10(f), acc);
>   }
>
>   sf = pow(10.0, 0.1*sf);
>
>   return(sf);
>}
>
>
>

Hi,
You can read this article for generate 1/f noise
http://www.istia.univ-angers.fr/~chapeau/papers/jautotel.pdf (in french)
http://www.istia.univ-angers.fr/~chapeau/papers/psiplrc2.pdf (in english)

ok




>Yu Xiaodong wrote:
>> Hi All,
>> 
>> Is there anyone know how to generate 1/f noise for simulating the
phase
>> noise. Someone told me though passing a white noise to a FIR the noise
can
>> be generated. Any clue about the FIR?
>> 
>> Thanks alot
>> 
>> R. Y.
>
>R.Y.:
>
>Here's some C++ code that I wrote that will create phase noise samples 
>from a user-specified phase noise profile. The code uses the GNU 
>Scientific Library (freely available).
>
>Good luck,
>OUP
>
>//--------------------------------------------------------------------------------------------------
>//
>//  phase_noise
>//
>//  Generate random deviates whose spectrum matches a phase noise 
>profile. This code uses the
>//  GNU Scientific Library (GSL) v1.1.1.
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>
>#include <stdio.h>
>#include <stdlib.h>
>#include <iostream.h>
>#include <math.h>
>
>#include <gsl/gsl_complex.h>
>#include <gsl/gsl_complex_math.h>
>#include <gsl/gsl_errno.h>
>#include <gsl/gsl_fft_complex.h>
>#include <gsl/gsl_fft_halfcomplex.h>
>#include <gsl/gsl_fft_real.h>
>#include <gsl/gsl_randist.h>
>#include <gsl/gsl_rng.h>
>#include <gsl/gsl_spline.h>
>#include <gsl/gsl_vector.h>
>
>#include "phase_noise.h"
>
>//--------------------------------------------------------------------------------------------------
>//
>//  main
>//
>//--------------------------------------------------------------------------------------------------
>
>int main (int argc, char *argv[])
>{
>   int i, numSamples;
>   double fs;
>   double xi, yi;
>
>   // number of samples
>   numSamples = atoi(argv[1]);
>
>   // sample frequency in Hz
>   fs = atof(argv[2]);
>
>   InitPhaseProfile();
>
>   double *phi = new double [numSamples];
>
>   GenPhaseNoiseDeviates(phi, numSamples, fs);
>
>   for (i = 0; i < numSamples; i++)
>   {
>     fprintf(stdout,"%d %19.16e\n", i, phi[i]);
>   }
>
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  GenPhaseNoiseDeviates
>//
>//  Inputs
>//  ------
>//  numSamples  number of samples to generate
>//  fs          sample frequency (a.k.a. PRF)
>//
>//  Input/Output
>//  ------------
>//  phiRandom   array to hold random deviates
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>int GenPhaseNoiseDeviates(double phiRandom[], int numSamples, double fs)
>{
>   int i;
>
>   extern int numPointsProfile;
>
>   //-------------------------------
>   // Generate white Gaussian noise
>   //-------------------------------
>
>   gsl_rng_env_setup();
>
>   const gsl_rng_type *T = gsl_rng_default;
>   gsl_rng *r = gsl_rng_alloc(T);
>
>   double *wn = new double [numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     wn[i] = gsl_ran_gaussian(r, 1.0);  // zero mean, unity sigma
>     //cout << wn[i] << "\n";
>   }
>
>   // Store the white noise samples in a GSL "packed" array. A "packed" 
>array holds complex data
>   // in a double arrray. The order is real, imag, real, imag, and so
on.
>
>   double *packedWn = new double [2*numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     REAL_PACKED(packedWn,i) = wn[i];
>     IMAG_PACKED(packedWn,i) = 0.0;
>   }
>
>   //------------------------------------------------------
>   // forward Fourier transform of the white noise samples
>   //------------------------------------------------------
>
>   gsl_fft_complex_workspace *workspace = 
>gsl_fft_complex_workspace_alloc(numSamples);
>   gsl_fft_complex_wavetable *wavetable = 
>gsl_fft_complex_wavetable_alloc(numSamples);  // we'll use this for the 
>inverse too
>
>   gsl_fft_complex_forward(packedWn, 1, numSamples, wavetable,
workspace);
>
>   //--------------------------------------------
>   // Construct the aliased phase noise response
>   //--------------------------------------------
>
>   int nFold = int( fInputPhaseProfile[numPointsProfile-1] / fs ) + 1;
>
>   if (nFold > nFoldMax) nFold = nFoldMax;
>
>   int nHalf = numSamples/2;
>
>   double *hMag = new double [nHalf];
>   double *f    = new double [nHalf];
>   double *mSSB = new double [nHalf];
>
>   double fMin = 1.0;
>   double fMax = fs / 2.0;
>   double fStep = (fMax - fMin)/(nHalf-1);
>
>   int k;
>   double summation;
>
>   for (i = 0; i < nHalf; i++)
>   {
>     summation = 0.0;
>     f[i] = fMin + i*fStep;
>     for (k = -nFold; k <= nFold; k++)
>     {
>       summation = summation + GetPhaseNoiseValue(f[i] + k*fs);
>     }
>     hMag[i] = summation;
>     mSSB[i] = GetPhaseNoiseValue(f[i]);
>   }
>
>//   for (i = 0; i < nHalf; i++)
>//   {
>//     fprintf(stdout,"%f %19.16e %19.16e\n", f[i], mSSB[i], hMag[i]);
>//   }
>
>   // Construct a double-sided phase noise profile. This has even 
>symmetry. Note that the
>   // GSL stores the positive frequencies first, then the negative 
>frequencies in reverse
>   // order. We'll construct our DSB profile the same way.
>
>   double *mDSB = new double [numSamples];
>
>   for (i = 0; i < nHalf; i++)
>   {
>     mDSB[i]       = hMag[i];
>     mDSB[i+nHalf] = hMag[nHalf - 1 - i];
>   }
>
>   double *packedDSB = new double [2*numSamples];
>
>   for (i = 0; i < numSamples; i++)
>   {
>     REAL_PACKED(packedDSB,i) = mDSB[i];
>     IMAG_PACKED(packedDSB,i) = 0.0;
>     //cout << mDSB[i] << "\n";
>   }
>
>   //----------------------------
>   // Convolve the two sequences
>   //----------------------------
>
>   double *packedY = new double [2*numSamples];
>
>   gsl_complex h, x, z;
>//   double hx, hy, wx, wy;
>
>   for (i = 0; i < numSamples; i++)
>   {
>//     wx = REAL_PACKED(packedWn,i);
>//     wy = IMAG_PACKED(packedWn,i);
>
>//     hx = REAL_PACKED(packedDSB,i);
>//     hy = IMAG_PACKED(packedDSB,i);
>
>//     GSL_SET_COMPLEX(&x, wx, wy);
>//     GSL_SET_COMPLEX(&h, hx, hy);
>
>     GSL_SET_COMPLEX(&x, REAL_PACKED(packedWn,i),
IMAG_PACKED(packedWn,i));
>     GSL_SET_COMPLEX(&h, REAL_PACKED(packedDSB,i), 
>IMAG_PACKED(packedDSB,i));
>
>     z = gsl_complex_mul(h,x);
>
>     REAL_PACKED(packedY,i) = GSL_REAL(z);
>     IMAG_PACKED(packedY,i) = GSL_IMAG(z);
>
>//     packedY[i] = packedDSB[i] * packedWn[i];  // this only works 
>because imaginary part is zero!
>   }
>
>   //---------------------------
>   // Inverse Fourier transform
>   //---------------------------
>
>   gsl_fft_complex_inverse(packedY, 1, numSamples, wavetable,
workspace);
>
>   for (i = 0; i < numSamples; i++)
>   {
>     phiRandom[i] = REAL_PACKED(packedY,i);
>   }
>
>   // free some stuff
>
>   gsl_fft_complex_workspace_free(workspace);
>   gsl_fft_complex_wavetable_free(wavetable);
>
>   gsl_spline_free (spline);
>   gsl_interp_accel_free(acc);
>
>   return 0;
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  InitPhaseProfile
>//
>//  Initialize stuff for the phase profile interpolation.
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>int InitPhaseProfile()
>{
>   int i;
>   extern int numPointsProfile;
>
>   numPointsProfile = sizeof(fInputPhaseProfile)/sizeof(double);
>
>   double *fLog = new double [numPointsProfile];
>
>   extern double fInputPhaseProfile[], mInputPhaseProfile[];
>   extern gsl_interp_accel *acc;
>   extern gsl_spline *spline;
>
>   for (i = 0; i < numPointsProfile; i++)
>   {
>     fLog[i] = log10(fInputPhaseProfile[i]);
>   }
>
>   //   printf ("#m=0,S=2\n");
>   //   for (i = 0; i < numPointsProfile; i++)
>   //   {
>   //     printf ("%g %g\n", fInputPhaseProfile[i],
mInputPhaseProfile[i]);
>   //   }
>
>   acc    = gsl_interp_accel_alloc ();
>   spline = gsl_spline_alloc (gsl_interp_linear, numPointsProfile);
>
>   gsl_spline_init (spline, fLog, mInputPhaseProfile, numPointsProfile);
>
>   return 0;
>}
>
>//--------------------------------------------------------------------------------------------------
>//
>//  GetPhaseNoiseValue
>//
>//  Given the desired frequency, this function uses linear interpolation

>in log-space to
>//  calculate the phase noise power.
>//
>//  Input
>//  -----
>//  f      frequency (Hz)
>//
>//  Output
>//  ------
>//  sf     phase noise power in absolute units
>//
>//  Revisions
>//  ---------
>//  25-Jun-02  Initial version written.
>//
>//--------------------------------------------------------------------------------------------------
>double GetPhaseNoiseValue(double f)
>{
>   double sf;
>   const double noiseFloor = -150.0;
>
>   extern int numPointsProfile;
>   extern double fInputPhaseProfile[];
>   extern gsl_interp_accel *acc;
>   extern gsl_spline *spline;
>
>   f = abs(f);
>
>   if (f > fInputPhaseProfile[numPointsProfile-1])
>   {
>     sf = noiseFloor;   // dBc/Hz
>   }
>   else
>   {
>     // interpolate to get phase noise power at desired frequency
>     sf = gsl_spline_eval (spline, log10(f), acc);
>   }
>
>   sf = pow(10.0, 0.1*sf);
>
>   return(sf);
>}
>
>
>

Flat white noise fed into a flat FM modulator will produce 1/F^2
shaped phase noise.

Most phase noise is usually 1/F^2.

1/F noise (flicker noise)  in transistors etc usually creates phase
noise that is 1/F^3 very close to the carrier.

Why do you want to create 1/F shaped phase noise?

Mark