git.gag.com Git - debian/openrocket/blob - core/src/net/sf/openrocket/util/PinkNoise.java

   1 package net.sf.openrocket.util;
   2 import java.util.Random;
   3
   4
   5 /**
   6  * A class that provides a source of pink noise with a power spectrum density
   7  * proportional to 1/f^alpha.  The values are computed by applying an IIR filter to
   8  * generated Gaussian random numbers.  The number of poles used in the filter may be
   9  * specified.  Values as low as 3 produce good results, but using a larger number of
  10  * poles allows lower frequencies to be amplified.  Below the cutoff frequency the
  11  * power spectrum density if constant.
  12  * <p>
  13  * The IIR filter use by this class is presented by N. Jeremy Kasdin, Proceedings of
  14  * the IEEE, Vol. 83, No. 5, May 1995, p. 822.
  15  *
  16  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  17  */
  18 public class PinkNoise {
  19         private final int poles;
  20         private final double[] multipliers;
  21
  22         private final double[] values;
  23         private final Random rnd;
  24
  25
  26         /**
  27          * Generate pink noise with alpha=1.0 using a five-pole IIR.
  28          */
  29         public PinkNoise() {
  30                 this(1.0, 5, new Random());
  31         }
  32
  33
  34         /**
  35          * Generate a specific pink noise using a five-pole IIR.
  36          *
  37          * @param alpha         the exponent of the pink noise, 1/f^alpha.
  38          */
  39         public PinkNoise(double alpha) {
  40                 this(alpha, 5, new Random());
  41         }
  42
  43
  44         /**
  45          * Generate pink noise specifying alpha and the number of poles.  The larger the
  46          * number of poles, the lower are the lowest frequency components that are amplified.
  47          *
  48          * @param alpha         the exponent of the pink noise, 1/f^alpha.
  49          * @param poles         the number of poles to use.
  50          */
  51         public PinkNoise(double alpha, int poles) {
  52                 this(alpha, poles, new Random());
  53         }
  54
  55
  56         /**
  57          * Generate pink noise specifying alpha, the number of poles and the randomness source.
  58          *
  59          * @param alpha    the exponent of the pink noise, 1/f^alpha.
  60          * @param poles    the number of poles to use.
  61          * @param random   the randomness source.
  62          */
  63         public PinkNoise(double alpha, int poles, Random random) {
  64                 this.rnd = random;
  65                 this.poles = poles;
  66                 this.multipliers = new double[poles];
  67                 this.values = new double[poles];
  68
  69                 double a = 1;
  70                 for (int i=0; i < poles; i++) {
  71                         a = (i - alpha/2) * a / (i+1);
  72                         multipliers[i] = a;
  73                 }
  74
  75                 // Fill the history with random values
  76                 for (int i=0; i < 5*poles; i++)
  77                         this.nextValue();
  78         }
  79
  80
  81
  82         public double nextValue() {
  83                 double x = rnd.nextGaussian();
  84 //              double x = rnd.nextDouble()-0.5;
  85
  86                 for (int i=0; i < poles; i++) {
  87                         x -= multipliers[i] * values[i];
  88                 }
  89                 System.arraycopy(values, 0, values, 1, values.length-1);
  90                 values[0] = x;
  91
  92                 return x;
  93         }
  94
  95
  96         public static void main(String[] arg) {
  97
  98                 PinkNoise source;
  99
 100                 source = new PinkNoise(1.0, 100);
 101                 double std = 0;
 102                 for (int i=0; i < 1000000; i++) {
 103
 104                 }
 105
 106
 107 //              int n = 5000000;
 108 //              double avgavg=0;
 109 //              double avgstd = 0;
 110 //              double[] val = new double[n];
 111 //
 112 //              for (int j=0; j < 10; j++) {
 113 //                      double avg=0, std=0;
 114 //                      source = new PinkNoise(5.0/3.0, 2);
 115 //
 116 //                      for (int i=0; i < n; i++) {
 117 //                              val[i] = source.nextValue();
 118 //                              avg += val[i];
 119 //                      }
 120 //                      avg /= n;
 121 //                      for (int i=0; i < n; i++) {
 122 //                              std += (val[i]-avg)*(val[i]-avg);
 123 //                      }
 124 //                      std /= n;
 125 //                      std = Math.sqrt(std);
 126 //
 127 //                      System.out.println("avg:"+avg+" stddev:"+std);
 128 //                      avgavg += avg;
 129 //                      avgstd += std;
 130 //              }
 131 //              avgavg /= 10;
 132 //              avgstd /= 10;
 133 //              System.out.println("Average avg:"+avgavg+" std:"+avgstd);
 134 //
 135                 // Two poles:
 136
 137         }
 138
 139
 140 }