>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);
>}
>
>
>