src/noise/FalphaNoise.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 /// Framework include files
0015 #include <DDDigi/noise/FalphaNoise.h>
0016
0017 /// C/C++ include files
0018 #include <cmath>
0019 #include <iostream>
0020 #include <stdexcept>
0021
0022 using namespace dd4hep::detail;
0023
0024 #ifdef  __GSL_FALPHA_NOISE
0025 #ifndef __CNOISE_C
0026 #define __CNOISE_C
0027
0028 /*
0029 * This code is distributed under GPLv3
0030 *
0031 * This code generates correlated or colored noise with a 1/f^alpha power
0032 * spectrum distribution.
0033 *
0034 * It uses the algorithm by:
0035 *
0036 * Jeremy Kasdin
0037 * Discrete Simulation of Colored Noise and Stochastic Processes and $ 1/f^\alpha $ Power Law Noise Generation
0038 * Proceedings of the IEEE
0039 * Volume 83, Number 5, 1995, pages 802-827.
0040 *
0041 * This code uses GSL fast Fourier transform gsl_fft_complex_forward(...)
0042 * and the GCC rand() functions
0043 *
0044 * Code Author: Miroslav Stoyanov, Jan 2011
0045 *
0046 * Copyright (C) 2011  Miroslav Stoyanov
0047 *
0048 * This program is free software: you can redistribute it and/or modify
0049 * it under the terms of the GNU General Public License as published by
0050 * the Free Software Foundation, either version 3 of the License, or
0051 * (at your option) any later version.
0052 *
0053 * This program is distributed in the hope that it will be useful,
0054 * but WITHOUT ANY WARRANTY; without even the implied warranty of
0055 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0056 * GNU General Public License for more details.
0057 *
0058 * Since the GNU General Public License is longer than this entire code,
0059 * a copy of it can be obtained separately at <http://www.gnu.org/licenses/>
0060 *
0061 * updated June 2011: fixed a small bug causing "nan" to be returned sometimes
0062 *
0063 */
0064
0065 //#include <cnoise.h>
0066 #ifndef __CNOISE_H
0067 #define __CNOISE_H
0068
0069 /*
0070 * This code isdistrubuted under GPLv3
0071 *
0072 * This code generates corelated or colored noise with 1/f^alpha power spectrum distribution.
0073 * It uses the algorithm by:
0074 *
0075 * Jeremy Kasdin
0076 * Discrete Simulation of Colored Noise and Stochastic Processes and $ 1/f^\alpha $ Power Law Noise Generation
0077 * Proceedings of the IEEE
0078 * Volume 83, Number 5, 1995, pages 802-827.
0079 *
0080 * This code uses GSL fast Fourier transform gsl_fft_complex_forward(...) and the GCC rand() functions
0081 *
0082 * Code Author: Miroslav Stoyanov, Jan 2011
0083 *
0084 * Copyright (C) 2011  Miroslav Stoyanov
0085 *
0086 * This program is free software: you can redistribute it and/or modify
0087 * it under the terms of the GNU General Public License as published by
0088 * the Free Software Foundation, either version 3 of the License, or
0089 * (at your option) any later version.
0090 *
0091 * This program is distributed in the hope that it will be useful,
0092 * but WITHOUT ANY WARRANTY; without even the implied warranty of
0093 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0094 * GNU General Public License for more details.
0095 *
0096 * Since the GNU General Public License is longer than this entire code,
0097 * a copy of it can be obtained separately at <http://www.gnu.org/licenses/>
0098 *
0099 */
0100
0101 #include <stdlib.h>
0102 #include <stdio.h>
0103 #include <math.h>
0104 #include <gsl/gsl_errno.h>
0105 #include <gsl/gsl_fft_complex.h>
0106
0107 double cnoise_uniform01 ( );
0108
0109 void cnoise_two_gaussian_truncated ( double *a, double *b, double std, double range );
0110
0111 void cnoise_generate_colored_noise ( double *x, int N, double alpha, double std );
0112
0113 // generates a vector x of size N with a 1/f^alpha frequency distribution
0114 // std is the standard deviation of the underlying gaussian distribution
0115 // Variables: double *x - Input: allocated vector of size N
0116 //                        Output: a realization of 1/f^alpha noise vector of size N, using an undelying Gaussian (0,std)
0117 //             int N    - Input: the size of the vector
0118 //            double alpha - Input: the decay rate for the power spectrum (i.e. 1/f^alpha
0119 //            double std - Input: the standard deviation of the underlying Gaussian distribution
0120 // NOTE: you should call srand( seed ) before you call dcnoise_generate_colored_noise
0121
0122 void cnoise_generate_colored_noise_uniform( double *x, int N, double alpha, double range );
0123 // same as cnoise_generate_colored_noise_uniform, except the white noise vector comes from
0124 // uniform distribution on (-range,+range)
0125
0126 void cnoise_generate_colored_noise_truncated( double *x, int N, double alpha, double std, double range );
0127 // same as cnoise_generate_colored_noise_uniform, except the white noise vector comes from
0128 // truncated Gaussian distribution with mean 0, standard deviation std and truncated to (-range,range)
0129
0130
0131
0132 #endif
0133
0134
0135 #define _CNOISE_PI 3.14159265358979323846 // pi
0136
0137 /******************************************************************************/
0138
0139 double cnoise_uniform01 ( )
0140
0141 /******************************************************************************/
0142 {
0143     return ( ( (double) rand() + 1.0 ) / ( (double) RAND_MAX + 1.0 ) );
0144 };
0145 /******************************************************************************/
0146
0147 void cnoise_two_gaussian_truncated ( double *a, double *b, double std, double range )
0148
0149 /******************************************************************************/
0150 {
0151     double gauss_u = cnoise_uniform01();
0152     double gauss_v = cnoise_uniform01();
0153     *a = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v );
0154     *b = std * sqrt( - 2 * log( gauss_u ) ) * sin( 2 * _CNOISE_PI * gauss_v );
0155     while ( (*a < -range) || (*a > range) ){
0156         gauss_u = cnoise_uniform01(); gauss_v = cnoise_uniform01();
0157         *a = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v );
0158     };
0159     while ( (*b < -range) || (*b > range) ){
0160         gauss_u = cnoise_uniform01(); gauss_v = cnoise_uniform01();
0161         *b = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v );
0162     };
0163 };
0164 /******************************************************************************/
0165
0166 void cnoise_generate_colored_noise ( double *x, int N, double alpha, double std )
0167
0168 /******************************************************************************/
0169 {
0170     int i;
0171     double tmp;
0172     double gauss_u, gauss_v;
0173     double * fh = (double*)malloc( 4*N*sizeof(double) );
0174     double * fw = (double*)malloc( 4*N*sizeof(double) );
0175
0176     fh[0] = 1.0; fh[1] = 0.0;
0177     for( i=1; i<N; i++ ){
0178         fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i );
0179         fh[2*i+1] = 0.0;
0180     };
0181     for( i=0; i<N; i+=2 ){
0182         gauss_u = cnoise_uniform01();
0183     gauss_v = cnoise_uniform01();
0184         fw[2*i]   = std * sqrt( - 2 * log( gauss_u ) ) * cos( 2 * _CNOISE_PI * gauss_v );
0185     fw[2*i+1] = 0.0;
0186         fw[2*i+2] = std * sqrt( - 2 * log( gauss_u ) ) * sin( 2 * _CNOISE_PI * gauss_v );
0187     fw[2*i+3] = 0.0;
0188     };
0189     for( i=2*N; i<4*N; i++ ){
0190     fh[i] = fw[i] = 0.0;
0191   };
0192
0193     gsl_fft_complex_wavetable* wavetable = gsl_fft_complex_wavetable_alloc(2*N);
0194     gsl_fft_complex_workspace* workspace = gsl_fft_complex_workspace_alloc(2*N);
0195
0196     gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace);
0197     gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace);
0198
0199     for( i=0; i<N+1; i++ ){
0200         tmp = fw[2*i];
0201         fw[2*i]   = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1];
0202         fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i];
0203     };
0204
0205     fw[0] /= 2.0; fw[1] /= 2.0;
0206     fw[2*N] /= 2.0; fw[2*N+1] /= 2.0;
0207
0208     for( i=N+1; i<2*N; i++ ){
0209         fw[2*i] = 0.0; fw[2*i+1] = 0.0;
0210     };
0211
0212     gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace);
0213
0214     for( i=0; i<N; i++ ){
0215         x[i] = 2.0*fw[2*i];
0216     };
0217
0218     free( fh );
0219     free( fw );
0220
0221     gsl_fft_complex_wavetable_free (wavetable);
0222     gsl_fft_complex_workspace_free (workspace);
0223 };
0224 /******************************************************************************/
0225
0226 void cnoise_generate_colored_noise_uniform( double *x, int N, double alpha,
0227   double range )
0228
0229 /******************************************************************************/
0230 {
0231     gsl_fft_complex_wavetable * wavetable;
0232     gsl_fft_complex_workspace * workspace;
0233
0234     int i;
0235     double tmp;
0236     double * fh = (double*)malloc( 4*N*sizeof(double) );
0237     double * fw = (double*)malloc( 4*N*sizeof(double) );
0238
0239     fh[0] = 1.0; fh[1] = 0.0;
0240     fw[0] = 2.0*range*cnoise_uniform01()-range; fw[1] = 0.0;
0241     for( i=1; i<N; i++ ){
0242         fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i );
0243         fh[2*i+1] = 0.0;
0244         fw[2*i] = 2.0*range*cnoise_uniform01()-range; fw[2*i+1] = 0.0;
0245     };
0246     for( i=2*N; i<4*N; i++ ){ fh[i] = 0.0; fw[i] = 0.0; };
0247
0248     wavetable = gsl_fft_complex_wavetable_alloc(2*N);
0249     workspace = gsl_fft_complex_workspace_alloc(2*N);
0250
0251     gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace);
0252     gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace);
0253
0254     for( i=0; i<N+1; i++ ){
0255         tmp = fw[2*i];
0256         fw[2*i]   = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1];
0257         fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i];
0258     };
0259
0260     fw[0] /= 2.0; fw[1] /= 2.0;
0261     fw[2*N] /= 2.0; fw[2*N+1] /= 2.0;
0262
0263     for( i=N+1; i<2*N; i++ ){
0264         fw[2*i] = 0.0; fw[2*i+1] = 0.0;
0265     };
0266
0267     gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace);
0268
0269     for( i=0; i<N; i++ ){
0270         x[i] = 2.0*fw[2*i];
0271     };
0272
0273     free( fh );
0274     free( fw );
0275
0276     gsl_fft_complex_wavetable_free (wavetable);
0277     gsl_fft_complex_workspace_free (workspace);
0278 };
0279 /******************************************************************************/
0280
0281 void cnoise_generate_colored_noise_truncated( double *x, int N, double alpha, double std, double range )
0282
0283 /******************************************************************************/
0284 {
0285     gsl_fft_complex_wavetable * wavetable;
0286     gsl_fft_complex_workspace * workspace;
0287
0288     int i;
0289     double tmp;
0290     double * fh = (double*)malloc( 4*N*sizeof(double) );
0291     double * fw = (double*)malloc( 4*N*sizeof(double) );
0292
0293     fh[0] = 1.0; fh[1] = 0.0;
0294     for( i=1; i<N; i++ ){
0295         fh[2*i] = fh[2*(i-1)] * ( 0.5 * alpha + (double)(i-1) ) / ((double) i );
0296         fh[2*i+1] = 0.0;
0297     };
0298     for( i=0; i<N; i+=2 ){
0299         cnoise_two_gaussian_truncated( &fw[2*i], &fw[2*i+2], std, range );
0300     };
0301     for( i=2*N; i<4*N; i++ ){ fh[i] = 0.0; fw[i] = 0.0; };
0302
0303     wavetable = gsl_fft_complex_wavetable_alloc(2*N);
0304     workspace = gsl_fft_complex_workspace_alloc(2*N);
0305
0306     gsl_fft_complex_forward (fh, 1, 2*N, wavetable, workspace);
0307     gsl_fft_complex_forward (fw, 1, 2*N, wavetable, workspace);
0308
0309     for( i=0; i<N+1; i++ ){
0310         tmp = fw[2*i];
0311         fw[2*i]   = tmp*fh[2*i] - fw[2*i+1]*fh[2*i+1];
0312         fw[2*i+1] = tmp*fh[2*i+1] + fw[2*i+1]*fh[2*i];
0313     };
0314
0315     fw[0] /= 2.0; fw[1] /= 2.0;
0316     fw[2*N] /= 2.0; fw[2*N+1] /= 2.0;
0317
0318     for( i=N+1; i<2*N; i++ ){
0319         fw[2*i] = 0.0; fw[2*i+1] = 0.0;
0320     };
0321
0322     gsl_fft_complex_inverse( fw, 1, 2*N, wavetable, workspace);
0323
0324     for( i=0; i<N; i++ ){
0325         x[i] = 2.0*fw[2*i];
0326     };
0327
0328     free ( fh );
0329     free ( fw );
0330
0331     gsl_fft_complex_wavetable_free (wavetable);
0332     gsl_fft_complex_workspace_free (workspace);
0333 };
0334
0335 #endif
0336 #endif
0337
0338 /// Initializing constructor. Leaves the enerator full functional
0339 FalphaNoise::FalphaNoise(size_t poles, double alpha, double variance)   {
0340   init(poles, alpha, variance);
0341 }
0342
0343 /// Initializing constructor with 5 poles. Leaves the enerator full functional
0344 FalphaNoise::FalphaNoise(double alpha, double variance)   {
0345   init(5, alpha, variance);
0346 }
0347
0348 /// Initialize the 1/f**alpha random generator. If already called reconfigures.
0349 void FalphaNoise::init(size_t poles, double alpha, double variance)    {
0350   m_poles    = poles;
0351   m_alpha    = alpha;
0352   m_variance = variance;
0353   if ( alpha < 0e0 || alpha > 2e0 )   {
0354     throw std::runtime_error("FalphaNoise: Invalid value for alpha: must be inside the interval [0,2]");
0355   }
0356   else if ( variance < 1e-10 )   {
0357     throw std::runtime_error("FalphaNoise: Invalid variance. Must be bigger than NULL!");
0358   }
0359   m_values.clear();
0360   m_multipliers.clear();
0361   m_distribution = std::normal_distribution<double>(0.0, m_variance);
0362   m_values.resize(m_poles+1,0e0);
0363   m_multipliers.reserve(m_poles);
0364
0365   double val = 1e0;
0366   for ( size_t i=0; i<m_poles; ++i)   {
0367     val *= (double(i) - alpha/2e0) / double(i+1);
0368     m_multipliers.emplace_back(val);
0369   }
0370   // Fill the history with random values
0371   std::default_random_engine generator;
0372   for ( size_t i=0; i < 5*m_poles; i++)
0373     (*this)(generator);
0374 }
0375
0376 /// Approximative compute proper variance and apply computed value
0377 void FalphaNoise::normalize(size_t shots)    {
0378   std::default_random_engine engine;
0379   normalize(engine, shots);
0380 }
0381
0382 /// Internal usage to normalize using user defined random engine
0383 void FalphaNoise::normalizeVariance(const random_engine_wrapper& generator, size_t shots)  {
0384   double mean = 0e0, mean2 = 0e0, var;
0385   auto tmp    = m_multipliers;
0386   double val  = 1e0;
0387   /// We have to correct for the case of white noise: alpha=0
0388   for ( size_t i=0; i<m_poles; ++i)   {
0389     val *= (double(i) - 0.5) / double(i+1);
0390     tmp[i] = val;
0391   }
0392   for ( size_t i=0; i < shots; i++)  {
0393     var = (*this)(generator);
0394     mean  += var;
0395     mean2 += var*var;
0396   }
0397   mean  /= double(shots);
0398   var = std::sqrt(mean2/double(shots) - mean*mean);
0399   m_variance *= m_variance/var;
0400   m_distribution = std::normal_distribution<double>(0.0, m_variance);
0401 }
0402
0403 /// Retrieve the next random number of the sequence
0404 double FalphaNoise::compute(double rndm_value)   {
0405 #ifdef  __GSL_FALPHA_NOISE
0406   if ( arr == nullptr )  {
0407     arr = new double[1000];
0408   }
0409   if ( ptr < 0 )  {
0410     cout << "Alpha: " << m_alpha << " sigma: " << m_variance << endl;
0411     cnoise_generate_colored_noise (arr, 1000, m_alpha, m_variance);
0412     ptr = 999;
0413   }
0414   --ptr;
0415   return arr[ptr+1];
0416 #else
0417   for ( int i=m_poles-1; i >= 0; --i )   {
0418     rndm_value -= m_multipliers[i] * m_values[i];
0419     m_values[i+1] = m_values[i]; // highest index always drops off...
0420   }
0421   m_values[0] = rndm_value;
0422   return rndm_value;
0423 #endif
0424 }