fract_flick.cc

obs网络试验平台

/****************************************************************
**                     fract_flick.cc                          **
*****************************************************************
**
**
**  Filename:   fract_flick.cc
**  Purpose:    Create a fractional flicker noise random proccess
**              This is a random process whose power spectral density
**              is of the form S(f) = 1/f**gamma, where gamma is
**              a (fractional) exponent between 0 and 2 (0 < gamma < 2).
**              Note that we do not allow the cases gamma = 0 (white
**              noise) and gamma = 2 (random walk) here, as these are
**              dealt with more simply as special cases.
**
**              The user specifies the number of frequency decades for the
**              filter, the number of poles per decade, and supplies
**              pointers to arrays of appropriate size to hold filter
**              coefficients and states.  The highest frequency zero is
**              automatically placed at a frequency corresponding to 1/10
**              the sampling rate.  The lowest frequency zero and pole
**              are then determined by the number of frequency decades the
**              user supplies.  The generated noise will have the above PSD
**              over time intervals that do not exceed the order of
**              1/(lowest frequency).
**
**              The noise process is generated using a bank of filters
**              as described by Corsini and Saletti (IEEE Trans Instrum
**              and Measuremetn, vol. 37, no. 4, Dec. 1988), which is
**              a generalization of the Barnes/Jarvis/Greenhall technique
**              used for the case gamma = 1.
**
**
** Author:      G. M. Garner
**
******************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define TWO_PI   6.28318530718 /* 2*PI    */

void init_flick (int n, double h, double *a, double *b, double gamma)
/*
 * Initialize the arrays of filter coefficients
 *
 * n = number of frequency decades for the filter
 * h = number of poles per decade (may be a fraction;
 *     ceil(n*h) = N = number of filter stages)
 * a = pointer to array of poles (array size = N)
 *     (phi array of Barnes and Greenhall)
 * b = pointer to array of zeros (array size = N)
 *     (theta array of Barnes and Greenhall)
 * gamma = noise exponent in power spectral density (1/f**gamma)
 */
{
     int N;       /* filter order */
     int i;
     double f, ae, be;
     double *apt, *bpt;
     double W, R;

     N = ceil (n*h);
     *a = 0.13;  /* From Barnes and Greenhall paper */
     R = pow (10.0, 1/(2.0*h));
     W = ((1.0 - (*a))/sqrt(*a));
     *b = 0.0;

     for (i = 2, apt = a+1, bpt = b+1; i <= N; i++, apt++, bpt++)
     {
        W = W / pow (R, 2.0 - gamma);
        (*bpt) = 1.0 + ((0.5*W) * (W - sqrt((W * W) + 4.0)));
        W = W / pow (R, gamma);
        (*apt) = 1.0 + ((0.5*W) * (W - sqrt((W * W) + 4.0)));
     }
}


void flicker (int N, double *y, double *y_old, double *a, double *b)
/*
 * Inegrate fractional flicker noise generator over 1 time step
 *
 * N = number of filter stages (= n*h, see init_flick above)
 * y = pointer to array of current state vector plus white noise input in
 *     element zero (array size = N+1)
 * y_old = pointer to array of previous time step state vector plus white
 *         noise input in element zero (array size = N+1)
 * a = pointer to array of poles (array size = N)
 *     (phi array of Barnes and Greenhall)
 * b = pointer to array of zeros (array size = N)
 *     (theta array of Barnes and Greenhall)
 */
{

     int i;
     double *ypt, *y_oldpt, *apt, *bpt;

     *(y+1) = ( (*a) - (*b) ) * (*(y_old+1)) + (*(y));
     for (i = 2, ypt =  y+2, y_oldpt = y_old+2, apt = a+1, bpt = b+1; i <= N;
                                   i++, ypt++, y_oldpt++, apt++, bpt++)
     {
        *ypt = (*apt) * (*y_oldpt) + ( (*(ypt-1)) -
                                            (*bpt) * (*(y_oldpt - 1)) );
     }

     for (i = 0, ypt =  y, y_oldpt = y_old, apt = a, bpt = b; i <= N;
          i++, ypt++, y_oldpt++, apt++, bpt++)
        *y_oldpt = (*ypt);
}