This notebook contains solutions to exercises in Chapter 4: Noise
Copyright 2015 Allen Downey
# Get thinkdsp.py
import os
if not os.path.exists('thinkdsp.py'):
!wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/thinkdsp.py
import numpy as np
import matplotlib.pyplot as plt
from thinkdsp import decorate
``A Soft Murmur'' is a web site that plays a mixture of natural noise sources, including rain, waves, wind, etc. At http://asoftmurmur.com/about/ you can find their list of recordings, most of which are at http://freesound.org.
Download a few of these files and compute the spectrum of each signal. Does the power spectrum look like white noise, pink noise, or Brownian noise? How does the spectrum vary over time?
if not os.path.exists('132736__ciccarelli__ocean-waves.wav'):
!wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/132736__ciccarelli__ocean-waves.wav
from thinkdsp import read_wave
wave = read_wave('132736__ciccarelli__ocean-waves.wav')
wave.make_audio()
I chose a recording of ocean waves. I selected a short segment:
segment = wave.segment(start=1.5, duration=1.0)
segment.make_audio()
And here's its spectrum:
spectrum = segment.make_spectrum()
spectrum.plot_power()
decorate(xlabel='Frequency (Hz)',
ylabel='Power')
Amplitude drops off with frequency, so this might be red or pink noise. We can check by looking at the power spectrum on a log-log scale.
spectrum.plot_power()
loglog = dict(xscale='log', yscale='log')
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
This structure, with increasing and then decreasing amplitude, seems to be common in natural noise sources.
Above $f = 10^3$, it might be dropping off linearly, but we can't really tell.
To see how the spectrum changes over time, I'll select another segment:
segment2 = wave.segment(start=2.5, duration=1.0)
segment2.make_audio()
And plot the two spectrums:
spectrum2 = segment2.make_spectrum()
spectrum.plot_power(alpha=0.5)
spectrum2.plot_power(alpha=0.5)
decorate(xlabel='Frequency (Hz)',
ylabel='Power')
Here they are again, plotting power on a log-log scale.
spectrum.plot_power(alpha=0.5)
spectrum2.plot_power(alpha=0.5)
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
So the structure seems to be consistent over time.
We can also look at a spectrogram:
segment.make_spectrogram(512).plot(high=5000)
decorate(xlabel='Time(s)', ylabel='Frequency (Hz)')
Within this segment, the overall amplitude drops off, but the mixture of frequencies seems consistent.
Exercise: In a noise signal, the mixture of frequencies changes over time. In the long run, we expect the power at all frequencies to be equal, but in any sample, the power at each frequency is random.
To estimate the long-term average power at each frequency, we can break a long signal into segments, compute the power spectrum for each segment, and then compute the average across the segments. You can read more about this algorithm at http://en.wikipedia.org/wiki/Bartlett%27s_method.
Implement Bartlett's method and use it to estimate the power
spectrum for a noise wave. Hint: look at the implementation
of make_spectrogram
.
from thinkdsp import Spectrum
def bartlett_method(wave, seg_length=512, win_flag=True):
"""Estimates the power spectrum of a noise wave.
wave: Wave
seg_length: segment length
"""
# make a spectrogram and extract the spectrums
spectro = wave.make_spectrogram(seg_length, win_flag)
spectrums = spectro.spec_map.values()
# extract the power array from each spectrum
psds = [spectrum.power for spectrum in spectrums]
# compute the root mean power (which is like an amplitude)
hs = np.sqrt(sum(psds) / len(psds))
fs = next(iter(spectrums)).fs
# make a Spectrum with the mean amplitudes
spectrum = Spectrum(hs, fs, wave.framerate)
return spectrum
bartlett_method
makes a spectrogram and extracts spec_map
, which maps from times to Spectrum objects. It computes the PSD for each spectrum, adds them up, and puts the results into a Spectrum object.
psd = bartlett_method(segment)
psd2 = bartlett_method(segment2)
psd.plot_power()
psd2.plot_power()
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
Now we can see the relationship between power and frequency more clearly. It is not a simple linear relationship, but it is consistent across different segments, even in details like the notches near 5000 Hz, 6000 Hz, and above 10,000 Hz.
if not os.path.exists('BTC_USD_2013-10-01_2020-03-26-CoinDesk.csv'):
!wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/BTC_USD_2013-10-01_2020-03-26-CoinDesk.csv
import pandas as pd
df = pd.read_csv('BTC_USD_2013-10-01_2020-03-26-CoinDesk.csv',
parse_dates=[0])
df
Currency | Date | Closing Price (USD) | 24h Open (USD) | 24h High (USD) | 24h Low (USD) | |
---|---|---|---|---|---|---|
0 | BTC | 2013-10-01 | 123.654990 | 124.304660 | 124.751660 | 122.563490 |
1 | BTC | 2013-10-02 | 125.455000 | 123.654990 | 125.758500 | 123.633830 |
2 | BTC | 2013-10-03 | 108.584830 | 125.455000 | 125.665660 | 83.328330 |
3 | BTC | 2013-10-04 | 118.674660 | 108.584830 | 118.675000 | 107.058160 |
4 | BTC | 2013-10-05 | 121.338660 | 118.674660 | 121.936330 | 118.005660 |
... | ... | ... | ... | ... | ... | ... |
2354 | BTC | 2020-03-22 | 5884.340133 | 6187.042146 | 6431.873162 | 5802.553402 |
2355 | BTC | 2020-03-23 | 6455.454688 | 5829.352511 | 6620.858253 | 5694.198299 |
2356 | BTC | 2020-03-24 | 6784.318011 | 6455.450650 | 6863.602196 | 6406.037439 |
2357 | BTC | 2020-03-25 | 6706.985089 | 6784.325204 | 6981.720386 | 6488.111885 |
2358 | BTC | 2020-03-26 | 6721.495392 | 6697.948320 | 6796.053701 | 6537.856462 |
2359 rows × 6 columns
ys = df['Closing Price (USD)']
ts = df.index
from thinkdsp import Wave
wave = Wave(ys, ts, framerate=1)
wave.plot()
decorate(xlabel='Time (days)')
spectrum = wave.make_spectrum()
spectrum.plot_power()
decorate(xlabel='Frequency (1/days)',
ylabel='Power',
**loglog)
spectrum.estimate_slope()[0]
-1.733254093675893
Red noise should have a slope of -2. The slope of this PSD is close to 1.7, so it's hard to say if we should consider it red noise or if we should say it's a kind of pink noise.
A Geiger counter is a device that detects radiation. When an ionizing particle strikes the detector, it outputs a surge of current. The total output at a point in time can be modeled as uncorrelated Poisson (UP) noise, where each sample is a random quantity from a Poisson distribution, which corresponds to the number of particles detected during an interval.
Write a class called UncorrelatedPoissonNoise
that inherits from _Noise
and provides evaluate
. It should use np.random.poisson
to generate random values from a Poisson distribution. The parameter of this function, lam
, is the average number of particles during each interval. You can use the attribute amp
to specify lam
. For example, if the framerate is 10 kHz and amp
is 0.001, we expect about 10 “clicks” per second.
Generate about a second of UP noise and listen to it. For low values of amp
, like 0.001, it should sound like a Geiger counter. For higher values it should sound like white noise. Compute and plot the power spectrum to see whether it looks like white noise.
from thinkdsp import Noise
class UncorrelatedPoissonNoise(Noise):
"""Represents uncorrelated Poisson noise."""
def evaluate(self, ts):
"""Evaluates the signal at the given times.
ts: float array of times
returns: float wave array
"""
ys = np.random.poisson(self.amp, len(ts))
return ys
Here's what it sounds like at low levels of "radiation".
amp = 0.001
framerate = 10000
duration = 1
signal = UncorrelatedPoissonNoise(amp=amp)
wave = signal.make_wave(duration=duration, framerate=framerate)
wave.make_audio()
To check that things worked, we compare the expected number of particles and the actual number:
expected = amp * framerate * duration
actual = sum(wave.ys)
print(expected, actual)
10.0 11
Here's what the wave looks like:
wave.plot()
And here's its power spectrum on a log-log scale.
spectrum = wave.make_spectrum()
spectrum.plot_power()
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
Looks like white noise, and the slope is close to 0.
spectrum.estimate_slope().slope
-0.0011109546801804139
With a higher arrival rate, it sounds more like white noise:
amp = 1
framerate = 10000
duration = 1
signal = UncorrelatedPoissonNoise(amp=amp)
wave = signal.make_wave(duration=duration, framerate=framerate)
wave.make_audio()
It looks more like a signal:
wave.plot()
And the spectrum converges on Gaussian noise.
import matplotlib.pyplot as plt
def normal_prob_plot(sample, fit_color='0.8', **options):
"""Makes a normal probability plot with a fitted line.
sample: sequence of numbers
fit_color: color string for the fitted line
options: passed along to Plot
"""
n = len(sample)
xs = np.random.normal(0, 1, n)
xs.sort()
ys = np.sort(sample)
mean, std = np.mean(sample), np.std(sample)
fit_ys = mean + std * xs
plt.plot(xs, fit_ys, color='gray', alpha=0.5, label='model')
plt.plot(xs, ys, **options)
spectrum = wave.make_spectrum()
spectrum.hs[0] = 0
normal_prob_plot(spectrum.real, label='real')
decorate(xlabel='Normal sample',
ylabel='Power')
normal_prob_plot(spectrum.imag, label='imag', color='C1')
decorate(xlabel='Normal sample')
The algorithm in this chapter for generating pink noise is conceptually simple but computationally expensive. There are more efficient alternatives, like the Voss-McCartney algorithm. Research this method, implement it, compute the spectrum of the result, and confirm that it has the desired relationship between power and frequency.
Solution: The fundamental idea of this algorithm is to add up several sequences of random numbers that get updates at different sampling rates. The first source should get updated at every time step; the second source every other time step, the third source ever fourth step, and so on.
In the original algorithm, the updates are evenly spaced. In an alternative proposed at http://www.firstpr.com.au/dsp/pink-noise/, they are randomly spaced.
My implementation starts with an array with one row per timestep and one column for each of the white noise sources. Initially, the first row and the first column are random and the rest of the array is Nan.
nrows = 100
ncols = 5
array = np.empty((nrows, ncols))
array.fill(np.nan)
array[0, :] = np.random.random(ncols)
array[:, 0] = np.random.random(nrows)
array[0:6]
array([[0.39517038, 0.33346713, 0.80747793, 0.23050029, 0.12329066], [0.93480849, nan, nan, nan, nan], [0.61007844, nan, nan, nan, nan], [0.74423927, nan, nan, nan, nan], [0.43421634, nan, nan, nan, nan], [0.07031545, nan, nan, nan, nan]])
The next step is to choose the locations where the random sources change. If the number of rows is $n$, the number of changes in the first column is $n$, the number in the second column is $n/2$ on average, the number in the third column is $n/4$ on average, etc.
So the total number of changes in the matrix is $2n$ on average; since $n$ of those are in the first column, the other $n$ are in the rest of the matrix.
To place the remaining $n$ changes, we generate random columns from a geometric distribution with $p=0.5$. If we generate a value out of bounds, we set it to 0 (so the first column gets the extras).
p = 0.5
n = nrows
cols = np.random.geometric(p, n)
cols[cols >= ncols] = 0
cols
array([4, 1, 3, 2, 1, 3, 1, 2, 1, 3, 1, 1, 4, 2, 1, 2, 1, 1, 1, 1, 1, 0, 1, 1, 2, 1, 2, 0, 2, 1, 1, 2, 0, 1, 3, 1, 2, 4, 0, 1, 1, 4, 1, 3, 1, 1, 1, 3, 2, 1, 3, 0, 3, 3, 4, 1, 4, 2, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 1, 1, 3, 2, 1, 1, 1, 1, 4, 1, 4, 3, 1, 3, 1, 2, 1, 0, 1, 3, 2, 3, 3, 1, 0, 1, 1, 2, 2, 2, 2])
Within each column, we choose a random row from a uniform distribution. Ideally we would choose without replacement, but it is faster and easier to choose with replacement, and I doubt it matters.
rows = np.random.randint(nrows, size=n)
rows
array([81, 74, 76, 24, 12, 56, 95, 4, 1, 12, 95, 44, 50, 76, 95, 34, 4, 91, 39, 38, 51, 13, 48, 54, 30, 9, 31, 4, 28, 79, 75, 45, 29, 65, 53, 83, 12, 53, 66, 88, 63, 76, 97, 28, 60, 99, 27, 84, 11, 8, 30, 71, 18, 4, 6, 57, 19, 58, 74, 37, 32, 53, 12, 91, 72, 36, 26, 87, 0, 94, 74, 69, 48, 91, 12, 50, 49, 64, 39, 3, 37, 91, 30, 27, 81, 79, 39, 45, 76, 86, 47, 54, 79, 0, 36, 86, 76, 43, 48, 8])
Now we can put random values at rach of the change points.
array[rows, cols] = np.random.random(n)
array[0:6]
array([[0.75218982, 0.33346713, 0.33151149, 0.23050029, 0.12329066], [0.93480849, 0.54326904, nan, nan, nan], [0.61007844, nan, nan, nan, nan], [0.74423927, nan, nan, nan, 0.12887454], [0.32223056, 0.44575854, 0.08476646, 0.88175009, nan], [0.07031545, nan, nan, nan, nan]])
Next we want to do a zero-order hold to fill in the NaNs. NumPy doesn't do that, but Pandas does. So I'll create a DataFrame:
df = pd.DataFrame(array)
df.head()
0 | 1 | 2 | 3 | 4 | |
---|---|---|---|---|---|
0 | 0.752190 | 0.333467 | 0.331511 | 0.23050 | 0.123291 |
1 | 0.934808 | 0.543269 | NaN | NaN | NaN |
2 | 0.610078 | NaN | NaN | NaN | NaN |
3 | 0.744239 | NaN | NaN | NaN | 0.128875 |
4 | 0.322231 | 0.445759 | 0.084766 | 0.88175 | NaN |
And then use fillna
along the columns.
filled = df.fillna(method='ffill', axis=0)
filled.head()
0 | 1 | 2 | 3 | 4 | |
---|---|---|---|---|---|
0 | 0.752190 | 0.333467 | 0.331511 | 0.23050 | 0.123291 |
1 | 0.934808 | 0.543269 | 0.331511 | 0.23050 | 0.123291 |
2 | 0.610078 | 0.543269 | 0.331511 | 0.23050 | 0.123291 |
3 | 0.744239 | 0.543269 | 0.331511 | 0.23050 | 0.128875 |
4 | 0.322231 | 0.445759 | 0.084766 | 0.88175 | 0.128875 |
Finally we add up the rows.
total = filled.sum(axis=1)
total.head()
0 1.770959 1 2.163380 2 1.838650 3 1.978395 4 1.863380 dtype: float64
If we put the results into a Wave, here's what it looks like:
wave = Wave(total.values)
wave.plot()
Here's the whole process in a function:
def voss(nrows, ncols=16):
"""Generates pink noise using the Voss-McCartney algorithm.
nrows: number of values to generate
rcols: number of random sources to add
returns: NumPy array
"""
array = np.empty((nrows, ncols))
array.fill(np.nan)
array[0, :] = np.random.random(ncols)
array[:, 0] = np.random.random(nrows)
# the total number of changes is nrows
n = nrows
cols = np.random.geometric(0.5, n)
cols[cols >= ncols] = 0
rows = np.random.randint(nrows, size=n)
array[rows, cols] = np.random.random(n)
df = pd.DataFrame(array)
df.fillna(method='ffill', axis=0, inplace=True)
total = df.sum(axis=1)
return total.values
To test it I'll generate 11025 values:
ys = voss(11025)
ys
array([ 7.74931667, 8.11992447, 7.81759903, ..., 10.21400096, 9.58848738, 10.08793761])
And make them into a Wave:
wave = Wave(ys)
wave.unbias()
wave.normalize()
Here's what it looks like:
wave.plot()
As expected, it is more random-walk-like than white noise, but more random looking than red noise.
Here's what it sounds like:
wave.make_audio()
And here's the power spectrum:
spectrum = wave.make_spectrum()
spectrum.hs[0] = 0
spectrum.plot_power()
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
The estimated slope is close to -1.
spectrum.estimate_slope().slope
-0.9961871315528988
We can get a better sense of the average power spectrum by generating a longer sample:
seg_length = 64 * 1024
iters = 100
wave = Wave(voss(seg_length * iters))
len(wave)
6553600
And using Barlett's method to compute the average.
spectrum = bartlett_method(wave, seg_length=seg_length, win_flag=False)
spectrum.hs[0] = 0
len(spectrum)
32769
It's pretty close to a straight line, with some curvature at the highest frequencies.
spectrum.plot_power()
decorate(xlabel='Frequency (Hz)',
ylabel='Power',
**loglog)
And the slope is close to -1.
spectrum.estimate_slope().slope
-1.0020572884106542