How can I generate a sin wave with power law noise?

I am interested in simulating a signal with different types of power law noise, as entitled in wikipedia, "Color of Noise". (The audio examples are pretty cool).

Basically, white noise, which many are familiar with, has a power spectral distribution (psd) that changes as $$f^0$$. Pink noise has a psd that changes as $$f^{-1}$$, brown noise's psd changes as $$f^{-2}$$, and blue noise has a psd that changes as $$f^1$$. What I want is to be able to add a noise that I specify as $$\beta$$, where $$f^\beta$$.

It seems like this should be possible with the function WhiteNoiseProcess[]Help page here. Also, here are two relevant links:

A wolfram link where they show different distributions.

A SE question about pink noise that has a comment:

It seems straightforward to implement using WhiteNoiseProcess though I don't have confidence that I understand it well enough to post as an answer.

So here is the signal: $$s(t)=\sin[2\pi f_o t + \phi(t)]$$

where $$f_o$$ is nominal frequency, $$t$$ is time, and $$\phi(t)$$ is the phase deviation. I believe we want to add the noise to $$\phi(t)$$.

Below, I take a stab at what the code should be:

s0 = Table[{t, Cos[2*Pi*4*t]}, {t, 0, 1, .001}]; (* signal w/o noise *)
stDev = .008;
phi = 2*Pi*stDev*Accumulate[RandomVariate[NormalDistribution[0, 1], Length[s0]  ]  ];
s1 = Table[{t, Cos[2*Pi*4*t + phi[[(t + .001)/.001]]  ]}, {t, 0, 1, .001}];(* signal w/noise *)
ListLinePlot[{s0, s1}](* Plot the two signals *) Great! Blue shows the s0 without noise, and s1 shown in orange has some noise. But what color is this noise? Some code to take the Fourier transform, so that we can see what $$\beta$$ is:

ft = Transpose[{Table[f, {f, 0, Length[s0] - 1}], Abs[Fourier[phi ]] }][[1 ;; 500]]; (* Fourier of the Phi data *)
orangeLine = Table[{n, 10/n}, {n, 1, 500}];
ListLogLogPlot[{ft, orangeLine  }, PlotRange -> All, Joined -> True] It looks sorta like $$\beta=-1$$, so it might be pink noise.

TLDR

Is there an easy way to generate the different colors of noise, and put them in a $$\sin$$ wave?

• Maybe write an expression in the frequency domain giving the distribution you want and use the inverse transform to put that back into the time domain? – Bill May 9 at 20:09
• For the power spectrum you need to look at Abs[Fourier[phi]]^2, which gives you $\beta=-2$ in your case. The spectral power is the squared modulus of the Fourier coefficient. – Roman May 10 at 17:38
• As far as I know there is no simple (reductionistic) model that generates $1/f$ noise in particular, and $f^{-\beta}$ noise in general (except for $\beta=-2$ as you've discovered: Wiener process etc.). You can fake it in certain ways (as @Bill suggests), but never fully. Have a look at the classic Bak/Tang/Wiesenfeld paper: doi.org/10.1103/PhysRevLett.59.381 or the Wikipedia article en.wikipedia.org/wiki/Pink_noise – Roman May 10 at 17:44

Following the comment by @Bill above, and with suggestions from @Roman I have implemented a simple noise generator:

noiseColor[sigma_, beta_, numberOfPoints_] :=
Module[{data, sig = sigma, b = beta, numPoin = numberOfPoints},
data = RandomFunction[WhiteNoiseProcess[sig*Sqrt], {1, numPoin*2}];
data = Table[data[[2, 1, 1]][[n]]*(1/n^-(b/2)), {n, 1, numPoin - 1}]; (*random list *)
data = Prepend[data, 0];  (* Removes DC offset *)
data = Re[InverseFourier[data]];  (*Brings to time domain *)
Return[data]
];

Note that we Prepend a zero to remove the DC bias, and have used a $$\frac{b}{2}$$ power to account for the square in the spectral density. (Thanks @Roman for these comments).

There are two problem with this module that will need updating when I gain an understanding of how to fix them. Firstly, the standard deviation sigma only works for white noise. Secondly, it doesn't seem to work for Blue noise or Violet noise.

Ok, first test is the easy one. It show that when we request white noise with a standard deviation of 1, we get a mean that is $$\approx$$0 and a standard deviation that is $$\approx$$1 :

white = noiseColor[1, 0, 10001];
Mean[white]
StandardDeviation[white]

Now, confirming that it works:

formatting = {ImageSize -> Medium, PlotStyle -> {{Gray}, {Red}}, Joined -> True, Frame -> True, LabelStyle -> Black};
plot1 = ListLinePlot[white, formatting, FrameLabel -> {"time", "amplitude"}];
pdf = (Abs[Fourier[white]])^2; (*pdf = probability density function*)
redLine = Table[x^-0, {x, 1, 10000}];
plot2 = ListLogLogPlot[{Drop[pdf, 1], redLine}, PlotRange -> {{100, 5000}, All}, formatting, FrameLabel -> {"frequency", "amplitude"}];
GraphicsRow[{plot1, plot2}]

Brown and pink work pretty well too, as shown next. In these figures below, the red line shows a $$f^{-\beta}$$ relationship, and the power spectrum follows that line very well. This proves that it is working.

pink = noiseColor[1, -1, 10001];
formatting = {ImageSize -> Medium, PlotStyle -> {{Pink}, {Red}}, Joined -> True, Frame -> True, LabelStyle -> Black};
plot1 = ListLinePlot[pink, formatting, FrameLabel -> {"time", "amplitude"}];
pdf = (Abs[Fourier[pink]])^2;
redLine = Table[x^-1, {x, 1, 10000}];
plot2 = ListLogLogPlot[{Drop[pdf, 1], redLine}, PlotRange -> {{100, 5000}, All}, formatting, FrameLabel -> {"frequency", "amplitude"}];
GraphicsRow[{plot1, plot2}]

brown = noiseColor[1, -2, 10001];
formatting = {ImageSize -> Medium, PlotStyle -> {{Brown}, {Red}}, Joined -> True, Frame -> True, LabelStyle -> Black};
plot1 = ListLinePlot[brown, formatting, FrameLabel -> {"time", "amplitude"}];
pdf = (Abs[Fourier[brown]])^2;
redLine = Table[x^-2, {x, 1, 10000}];
plot2 = ListLogLogPlot[{Drop[pdf, 1], redLine}, PlotRange -> {{100, 5000}, All}, formatting, FrameLabel -> {"frequency", "amplitude"}];
GraphicsRow[{plot1, plot2}]

Moving on to Blue and Violet noise. In these figures below, the red line shows a $$f^-\beta$$ relationship, and the power spectrum for some reason does not follow that line very well. What is wrong?

blue = noiseColor[1, 1, 10001];
formatting = {ImageSize -> Medium, PlotStyle -> {{Blue}, {Red}}, Joined -> True, Frame -> True, LabelStyle -> Black};
plot1 = ListLinePlot[blue, formatting, FrameLabel -> {"time", "amplitude"}];
pdf = (Abs[Fourier[blue]])^2;
redLine = Table[x^1, {x, 1, 10000}];
plot2 = ListLogLogPlot[{Drop[pdf, 1], redLine}, PlotRange -> {{100, 5000}, All}, formatting, FrameLabel -> {"frequency", "amplitude"}];
GraphicsRow[{plot1, plot2}]

violet = noiseColor[1, 2, 10001];
formatting = {ImageSize -> Medium, PlotStyle -> {{Purple}, {Red}}, Joined -> True, Frame -> True, LabelStyle -> Black};
plot1 = ListLinePlot[violet, formatting, FrameLabel -> {"time", "amplitude"}];
pdf = (Abs[Fourier[violet]])^2;
redLine = Table[x^2, {x, 1, 10000}];
plot2 = ListLogLogPlot[{Drop[pdf, 1], redLine}, PlotRange -> {{100, 5000}, All}, formatting, FrameLabel -> {"frequency", "amplitude"}];
GraphicsRow[{plot1, plot2}]

Summary

So this answer presents a power law noise generator that can do the following things correctly:

1. Make zero-mean white noise with correct standard deviation, with $$f^0$$ power spectrum.
2. Make zero-mean pink noise, with $$f^{-1}$$ power spectrum.
3. Make zero-mean brown noise, and $$f^{-2}$$ power spectrum.

And does the following which is wrong:

1. Standard deviation for pink, brown, blue, and violet is incorrect
2. Make zero-mean blue noise, without $$f^{1}$$ power spectrum.
3. Make zero-mean violet noise, without $$f^{2}$$ power spectrum.
• The data you get from RandomFunction[] is in the time domain. I think @Bill's comment is to take the Fourier of the data and then do the coloring and eventually get it back to the time domain? – Anjan Kumar May 10 at 19:48
• Same comment as above: a spectral power dependence of $1/f^{\beta}$ means that the amplitude scales as $1/f^{\beta/2}$ and you need to use 1/n^(b/2) in your code to scale the amplitudes, not 1/n^b. – Roman May 10 at 20:04
• @Roman I think you are right about this. Any suggestions for a reference where I can read about the amplitude scaling by $1/f^{\beta/2}$ in depth? – axsvl77 May 10 at 21:48
• @axsvl77 the power spectral density is explained here: en.wikipedia.org/wiki/Spectral_density under "Power spectral density". You can see that there's a squared absolute value involved when relating the amplitude spectral density (the result of a Fourier transform) to the power spectral density. – Roman May 11 at 6:36
• You have an off-by-one error in your amplitudes. From the manual of InverseFourier: "Note that the zero frequency term must appear at position 1 in the input list." You have the frequency-1 term at position 1 though. – Roman May 11 at 10:53

This is a bit of a hack but approximates a power spectrum of $$f^{\beta}$$ by generating uniform complex Gaussian random noise amplitudes and weighing them by $$f^{\beta/2}$$, then inverse-Fourier-transforming to the time domain. I think there are much better algorithms out there though. This is a well-researched topic.

noise[n_Integer /; n >= 1, β_?NumericQ] :=
Re[InverseFourier[(RandomVariate[NormalDistribution[], {n, 2}].{1,I})*
Prepend[Range[n-1]^(β/2), 0]]]

Notice that I've forced the zero-frequency component to zero amplitude because otherwise we'd get an error for $$\beta<0$$.

Try it out: $$1/f$$ noise:

ListLinePlot[noise[10^3, -1]] You can add such generated noise vectors to a sine wave or whatever your signal is.

• Trying to get this to work. Can you tell me why ListLogLogPlot[Drop[Abs[Fourier[noise[10^4, 2]]]^2, 1], PlotRange -> All] Doesn't give a $f^2$ relationship? It seems flat. – axsvl77 May 13 at 21:39
• Something to do with taking the real part of the noise in the definition. If you remove the Re from the noise definition it works, right? Anyway this answer was meant as a quick hack, not as a good algorithm. – Roman May 13 at 21:46
• Works fine with ListLogLogPlot[Drop[Abs[Fourier[noise[10^4, -2]]]^2, 1], PlotRange -> All] though – axsvl77 May 14 at 1:10
• @axsvl77 the problem seems to be with $\beta>0$. – Roman May 14 at 11:34
• Could be. For mine, it works find before it gets in the time domain, but can't come back. – axsvl77 May 14 at 11:52