CoCalc -- systems.py

Path: examples/think-dsp/code / systems.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 thinkdsp
11
import thinkplot
12

13
import numpy as np
14

15

16
def plot_filter():
17
    """Plots the filter that corresponds to a 2-element moving average.
18
    """
19
    impulse = np.zeros(8)
20
    impulse[0] = 1
21
    wave = thinkdsp.Wave(impulse, framerate=8)
22
    
23
    impulse_spectrum = wave.make_spectrum(full=True)
24
    window_array = np.array([0.5, 0.5, 0, 0, 0, 0, 0, 0,])
25
    window = thinkdsp.Wave(window_array, framerate=8)
26
    filtr = window.make_spectrum(full=True)
27

28
    filtr.plot()
29
    thinkplot.config(xlabel='Frequency', ylabel='Amplitude')
30
    thinkplot.save('systems1')
31

32

33
def read_response():
34
    """Reads the impulse response file and removes the initial silence.
35
    """
36
    response = thinkdsp.read_wave('180960__kleeb__gunshot.wav')
37
    start = 0.26
38
    response = response.segment(start=start)
39
    response.shift(-start)
40
    response.normalize()
41
    return response
42

43

44
def plot_response(response):
45
    """Plots an input wave and the corresponding output.
46
    """
47
    thinkplot.preplot(cols=2)
48
    response.plot()
49
    thinkplot.config(xlabel='Time (s)',
50
                     xlim=[0.26, response.end],
51
                     ylim=[-1.05, 1.05])
52

53
    thinkplot.subplot(2)
54
    transfer = response.make_spectrum()
55
    transfer.plot()
56
    thinkplot.config(xlabel='Frequency (Hz)',
57
                     xlim=[0, 22500],
58
                     ylabel='Amplitude')
59
    thinkplot.save(root='systems6')
60

61
    violin = thinkdsp.read_wave('92002__jcveliz__violin-origional.wav')
62

63
    start = 0.11
64
    violin = violin.segment(start=start)
65
    violin.shift(-start)
66

67
    violin.truncate(len(response))
68
    violin.normalize()
69
    spectrum = violin.make_spectrum()
70

71
    output = (spectrum * transfer).make_wave()
72
    output.normalize()
73

74
    thinkplot.preplot(rows=2)
75
    violin.plot(label='input')
76
    thinkplot.config(xlim=[0, violin.end],
77
                     ylim=[-1.05, 1.05])
78

79
    thinkplot.subplot(2)
80
    output.plot(label='output')
81
    thinkplot.config(xlabel='Time (s)',
82
                     xlim=[0, violin.end],
83
                     ylim=[-1.05, 1.05])
84

85
    thinkplot.save(root='systems7')
86

87

88
def shifted_scaled(wave, shift, factor):
89
    res = wave.copy()
90
    res.shift(shift)
91
    res.scale(factor)
92
    return res
93

94

95
def plot_convolution(response):
96
    """Plots the impulse response and a shifted, scaled copy.
97
    """
98
    shift = 1
99
    factor = 0.5
100
    
101
    gun2 = response + shifted_scaled(response, shift, factor)
102
    gun2.plot()
103
    thinkplot.config(xlabel='Time (s)',
104
                     ylim=[-1.05, 1.05],
105
                     legend=False)
106
    thinkplot.save(root='systems8')
107

108

109
def plot_sawtooth(response):
110
    signal = thinkdsp.SawtoothSignal(freq=441)
111
    wave = signal.make_wave(duration=0.1, framerate=response.framerate)
112

113
    total = 0
114
    for t, y in zip(wave.ts, wave.ys):
115
        total += shifted_scaled(response, t, y)
116

117
    total.normalize()
118

119
    high = 5000
120
    wave.make_spectrum().plot(high=high, color='0.7', label='original')
121
    segment = total.segment(duration=0.2)
122
    segment.make_spectrum().plot(high=high, label='convolved')
123
    thinkplot.config(xlabel='Frequency (Hz)', ylabel='Amplitude')
124
    thinkplot.save(root='systems9')
125

126

127
def main():
128
    plot_filter()
129

130
    response = read_response()
131
    plot_response(response)
132
    plot_convolution(response)
133

134

135
if __name__ == '__main__':
136
    main()
137

138