1 package net.sf.openrocket.util;
2 import java.util.Random;
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.
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.
16 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
18 public class PinkNoise {
19 private final int poles;
20 private final double[] multipliers;
22 private final double[] values;
23 private final Random rnd;
27 * Generate pink noise with alpha=1.0 using a five-pole IIR.
30 this(1.0, 5, new Random());
35 * Generate a specific pink noise using a five-pole IIR.
37 * @param alpha the exponent of the pink noise, 1/f^alpha.
39 public PinkNoise(double alpha) {
40 this(alpha, 5, new Random());
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.
48 * @param alpha the exponent of the pink noise, 1/f^alpha.
49 * @param poles the number of poles to use.
51 public PinkNoise(double alpha, int poles) {
52 this(alpha, poles, new Random());
57 * Generate pink noise specifying alpha, the number of poles and the randomness source.
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.
63 public PinkNoise(double alpha, int poles, Random random) {
66 this.multipliers = new double[poles];
67 this.values = new double[poles];
70 for (int i=0; i < poles; i++) {
71 a = (i - alpha/2) * a / (i+1);
75 // Fill the history with random values
76 for (int i=0; i < 5*poles; i++)
82 public double nextValue() {
83 double x = rnd.nextGaussian();
84 // double x = rnd.nextDouble()-0.5;
86 for (int i=0; i < poles; i++) {
87 x -= multipliers[i] * values[i];
89 System.arraycopy(values, 0, values, 1, values.length-1);
96 public static void main(String[] arg) {
100 source = new PinkNoise(1.0, 100);
102 for (int i=0; i < 1000000; i++) {
109 // double avgstd = 0;
110 // double[] val = new double[n];
112 // for (int j=0; j < 10; j++) {
113 // double avg=0, std=0;
114 // source = new PinkNoise(5.0/3.0, 2);
116 // for (int i=0; i < n; i++) {
117 // val[i] = source.nextValue();
121 // for (int i=0; i < n; i++) {
122 // std += (val[i]-avg)*(val[i]-avg);
125 // std = Math.sqrt(std);
127 // System.out.println("avg:"+avg+" stddev:"+std);
133 // System.out.println("Average avg:"+avgavg+" std:"+avgstd);