CoCalc -- noise.py

Path: examples/think-dsp/code / noise.py
Views: ⁷³³⁶
1
"""This file contains code used in "Think DSP",
2
by Allen B. Downey, available from greenteapress.com
3

4
Copyright 2014 Allen B. Downey
5
License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html
6
"""
7

8
from __future__ import print_function, division
9

10
import numpy
11
import thinkstats2
12
import thinkdsp
13
import thinkplot
14

15

16
def process_noise(signal, root='white'):
17
    """Plots wave and spectrum for noise signals.
18

19
    signal: Signal
20
    root: string used to generate file names
21
    """
22
    framerate = 11025
23
    wave = signal.make_wave(duration=0.5, framerate=framerate)
24

25
    # 0: waveform
26
    segment = wave.segment(duration=0.1)
27
    segment.plot(linewidth=1, alpha=0.5)
28
    thinkplot.save(root=root+'noise0',
29
                   xlabel='Time (s)',
30
                   ylim=[-1.05, 1.05])
31

32
    spectrum = wave.make_spectrum()
33

34
    # 1: spectrum
35
    spectrum.plot_power(linewidth=1, alpha=0.5)
36
    thinkplot.save(root=root+'noise1',
37
                   xlabel='Frequency (Hz)',
38
                   ylabel='Power',
39
                   xlim=[0, spectrum.fs[-1]])
40

41
    slope, _, _, _, _ = spectrum.estimate_slope()
42
    print('estimated slope', slope)
43

44
    # 2: integrated spectrum
45
    integ = spectrum.make_integrated_spectrum()
46
    integ.plot_power()
47
    thinkplot.save(root=root+'noise2',
48
                   xlabel='Frequency (Hz)',
49
                   ylabel='Cumulative fraction of total power',
50
                   xlim=[0, framerate/2])
51

52
    # 3: log-log power spectrum
53
    spectrum.hs[0] = 0
54
    thinkplot.preplot(cols=2)
55
    spectrum.plot_power(linewidth=1, alpha=0.5)
56
    thinkplot.config(xlabel='Frequency (Hz)',
57
                     ylabel='Power',
58
                     xlim=[0, framerate/2])
59

60
    thinkplot.subplot(2)
61
    spectrum.plot_power(linewidth=1, alpha=0.5)
62
    thinkplot.config(xlabel='Frequency (Hz)',                  
63
                   xscale='log',
64
                   yscale='log',
65
                   xlim=[0, framerate/2])
66

67
    thinkplot.save(root=root+'noise3')
68

69

70
def plot_gaussian_noise():
71
    """Shows the distribution of the spectrum of Gaussian noise.
72
    """
73
    thinkdsp.random_seed(18)
74
    signal = thinkdsp.UncorrelatedGaussianNoise()
75
    wave = signal.make_wave(duration=0.5, framerate=11025)
76
    spectrum = wave.make_spectrum()
77

78
    thinkplot.preplot(2, cols=2)
79
    thinkstats2.NormalProbabilityPlot(spectrum.real, label='real')
80
    thinkplot.config(xlabel='Normal sample',
81
                     ylabel='Amplitude',
82
                     ylim=[-250, 250],
83
                     loc='lower right')
84

85
    thinkplot.subplot(2)
86
    thinkstats2.NormalProbabilityPlot(spectrum.imag, label='imag')
87
    thinkplot.config(xlabel='Normal sample',
88
                     ylim=[-250, 250],
89
                     loc='lower right')
90

91
    thinkplot.save(root='noise1')
92

93

94
def plot_pink_noise():
95
    """Makes a plot showing power spectrums for pink noise.
96
    """
97
    thinkdsp.random_seed(20)
98

99
    duration = 1.0
100
    framerate = 512
101

102
    def make_spectrum(signal):
103
        wave = signal.make_wave(duration=duration, framerate=framerate)
104
        spectrum = wave.make_spectrum()
105
        spectrum.hs[0] = 0
106
        return spectrum
107

108
    signal = thinkdsp.UncorrelatedUniformNoise()
109
    white = make_spectrum(signal)
110

111
    signal = thinkdsp.PinkNoise()
112
    pink = make_spectrum(signal)
113

114
    signal = thinkdsp.BrownianNoise()
115
    red = make_spectrum(signal)
116

117
    linewidth = 2
118
    # colorbrewer2.org 4-class sequential OrRd
119
    white.plot_power(label='white', color='#fdcc8a', linewidth=linewidth)
120
    pink.plot_power(label='pink', color='#fc8d59', linewidth=linewidth)
121
    red.plot_power(label='red', color='#d7301f', linewidth=linewidth)
122
    thinkplot.save(root='noise-triple',
123
                   xlabel='Frequency (Hz)',
124
                   ylabel='Power',
125
                   xscale='log',
126
                   yscale='log',
127
                   xlim=[1, red.fs[-1]])
128

129

130
def main():
131
    thinkdsp.random_seed(17)
132
    plot_pink_noise()
133

134
    thinkdsp.random_seed(17)
135
    plot_gaussian_noise()
136

137
    thinkdsp.random_seed(20)
138
    signal = thinkdsp.UncorrelatedUniformNoise()
139
    process_noise(signal, root='white')
140

141
    thinkdsp.random_seed(20)
142
    signal = thinkdsp.PinkNoise(beta=1.0)
143
    process_noise(signal, root='pink')
144

145
    thinkdsp.random_seed(17)
146
    signal = thinkdsp.BrownianNoise()
147
    process_noise(signal, root='red')
148

149

150
if __name__ == '__main__':
151
    main()
152

153