diff --git a/test_code/aliasing.py b/test_code/aliasing.py new file mode 100644 index 0000000..c51b19a --- /dev/null +++ b/test_code/aliasing.py @@ -0,0 +1,104 @@ +"""This file contains code used in "Think DSP", +by Allen B. Downey, available from greenteapress.com + +Copyright 2013 Allen B. Downey +License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html +""" + +from __future__ import print_function, division + +import thinkdsp +import thinkplot + +FORMATS = ['pdf', 'png'] + +def triangle_example(freq): + """Makes a figure showing a triangle wave. + + freq: frequency in Hz + """ + framerate = 10000 + signal = thinkdsp.TriangleSignal(freq) + + duration = signal.period*3 + segment = signal.make_wave(duration, framerate=framerate) + segment.plot() + thinkplot.save(root='triangle-%d-1' % freq, + xlabel='Time (s)', + axis=[0, duration, -1.05, 1.05]) + + wave = signal.make_wave(duration=0.5, framerate=framerate) + spectrum = wave.make_spectrum() + + thinkplot.preplot(cols=2) + spectrum.plot() + thinkplot.config(xlabel='Frequency (Hz)', + ylabel='Amplitude') + + thinkplot.subplot(2) + spectrum.plot() + thinkplot.config(ylim=[0, 500], + xlabel='Frequency (Hz)') + + thinkplot.save(root='triangle-%d-2' % freq) + + +def square_example(freq): + """Makes a figure showing a square wave. + + freq: frequency in Hz + """ + framerate = 10000 + signal = thinkdsp.SquareSignal(freq) + + duration = signal.period*3 + segment = signal.make_wave(duration, framerate=framerate) + segment.plot() + thinkplot.save(root='square-%d-1' % freq, + xlabel='Time (s)', + axis=[0, duration, -1.05, 1.05]) + + wave = signal.make_wave(duration=0.5, framerate=framerate) + spectrum = wave.make_spectrum() + spectrum.plot() + thinkplot.save(root='square-%d-2' % freq, + xlabel='Frequency (Hz)', + ylabel='Amplitude') + + +def aliasing_example(offset=0.000003): + """Makes a figure showing the effect of aliasing. + """ + framerate = 10000 + + def plot_segment(freq): + signal = thinkdsp.CosSignal(freq) + duration = signal.period*4 + thinkplot.Hlines(0, 0, duration, color='gray') + segment = signal.make_wave(duration, framerate=framerate*10) + segment.plot(linewidth=0.5, color='gray') + segment = signal.make_wave(duration, framerate=framerate) + segment.plot_vlines(label=freq, linewidth=4) + + thinkplot.preplot(rows=2) + plot_segment(4500) + thinkplot.config(axis=[-0.00002, 0.0007, -1.05, 1.05]) + + thinkplot.subplot(2) + plot_segment(5500) + thinkplot.config(axis=[-0.00002, 0.0007, -1.05, 1.05]) + + thinkplot.save(root='aliasing1', + xlabel='Time (s)', + formats=FORMATS) + + +def main(): + triangle_example(freq=200) + # triangle_example(freq=1100) + square_example(freq=100) + # aliasing_example() + + +if __name__ == '__main__': + main() diff --git a/test_code/chirp.py b/test_code/chirp.py new file mode 100644 index 0000000..4939b76 --- /dev/null +++ b/test_code/chirp.py @@ -0,0 +1,224 @@ +"""This file contains code used in "Think DSP", +by Allen B. Downey, available from greenteapress.com + +Copyright 2015 Allen B. Downey +License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html +""" + +from __future__ import print_function, division + +import math +import numpy as np + +import thinkdsp +import thinkplot + +import matplotlib.pyplot as pyplot + +import warnings +warnings.filterwarnings('ignore') + +PI2 = math.pi * 2 + + +def linear_chirp_evaluate(ts, low=440, high=880, amp=1.0): + """Computes the waveform of a linear chirp and prints intermediate values. + + low: starting frequency + high: ending frequency + amp: amplitude + """ + print('ts', ts) + + freqs = np.linspace(low, high, len(ts)-1) + print('freqs', freqs) + + dts = np.diff(ts) + print('dts', dts) + + dphis = np.insert(PI2 * freqs * dts, 0, 0) + print('dphis', dphis) + + phases = np.cumsum(dphis) + print('phases', phases) + + ys = amp * np.cos(phases) + print('ys', ys) + + return ys + + +def discontinuity(num_periods=30, hamming=False): + """Plots the spectrum of a sinusoid with/without windowing. + + num_periods: how many periods to compute + hamming: boolean whether to apply Hamming window + """ + signal = thinkdsp.SinSignal(freq=440) + duration = signal.period * num_periods + wave = signal.make_wave(duration) + + if hamming: + wave.hamming() + + print(len(wave.ys), wave.ys[0], wave.ys[-1]) + spectrum = wave.make_spectrum() + spectrum.plot(high=880) + + +def three_spectrums(): + """Makes a plot showing three spectrums for a sinusoid. + """ + thinkplot.preplot(rows=1, cols=3) + + pyplot.subplots_adjust(wspace=0.3, hspace=0.4, + right=0.95, left=0.1, + top=0.95, bottom=0.1) + + xticks = range(0, 900, 200) + + thinkplot.subplot(1) + thinkplot.config(xticks=xticks) + discontinuity(num_periods=30, hamming=False) + + thinkplot.subplot(2) + thinkplot.config(xticks=xticks, xlabel='Frequency (Hz)') + discontinuity(num_periods=30.25, hamming=False) + + thinkplot.subplot(3) + thinkplot.config(xticks=xticks) + discontinuity(num_periods=30.25, hamming=True) + + # thinkplot.save(root='windowing1') + thinkplot.show() + + +def window_plot(): + """Makes a plot showing a sinusoid, hamming window, and their product. + """ + signal = thinkdsp.SinSignal(freq=440) + duration = signal.period * 10.25 + wave1 = signal.make_wave(duration) + wave2 = signal.make_wave(duration) + + ys = np.hamming(len(wave1.ys)) + ts = wave1.ts + window = thinkdsp.Wave(ys, ts, wave1.framerate) + + wave2.hamming() + + thinkplot.preplot(rows=3, cols=1) + + pyplot.subplots_adjust(wspace=0.3, hspace=0.3, + right=0.95, left=0.1, + top=0.95, bottom=0.05) + + thinkplot.subplot(1) + wave1.plot() + thinkplot.config(axis=[0, duration, -1.07, 1.07]) + + thinkplot.subplot(2) + window.plot() + thinkplot.config(axis=[0, duration, -1.07, 1.07]) + + thinkplot.subplot(3) + wave2.plot() + thinkplot.config(axis=[0, duration, -1.07, 1.07], + xlabel='Time (s)') + + # thinkplot.save(root='windowing2') + thinkplot.show() + + +def chirp_spectrum(): + """Plots the spectrum of a one-second one-octave linear chirp. + """ + signal = thinkdsp.Chirp(start=220, end=440) + wave = signal.make_wave(duration=1) + + thinkplot.preplot(3, cols=3) + duration = 0.01 + wave.segment(0, duration).plot(xfactor=1000) + thinkplot.config(ylim=[-1.05, 1.05]) + + thinkplot.subplot(2) + wave.segment(0.5, duration).plot(xfactor=1000) + thinkplot.config(yticklabels='invisible', + xlabel='Time (ms)') + + thinkplot.subplot(3) + wave.segment(0.9, duration).plot(xfactor=1000) + thinkplot.config(yticklabels='invisible') + + # thinkplot.save('chirp3') + thinkplot.show() + + + spectrum = wave.make_spectrum() + spectrum.plot(high=700) + # thinkplot.save('chirp1', + # xlabel='Frequency (Hz)', + # ylabel='Amplitude') + thinkplot.show() + +def chirp_spectrogram(): + """Makes a spectrogram of a one-second one-octave linear chirp. + """ + signal = thinkdsp.Chirp(start=220, end=440) + wave = signal.make_wave(duration=1, framerate=11025) + spectrogram = wave.make_spectrogram(seg_length=512) + + print('time res', spectrogram.time_res) + print('freq res', spectrogram.freq_res) + print('product', spectrogram.time_res * spectrogram.freq_res) + + spectrogram.plot(high=700) + + # thinkplot.save('chirp2', + # xlabel='Time (s)', + # ylabel='Frequency (Hz)') + thinkplot.show() + +def overlapping_windows(): + """Makes a figure showing overlapping hamming windows. + """ + n = 256 + window = np.hamming(n) + + thinkplot.preplot(num=5) + start = 0 + for i in range(5): + xs = np.arange(start, start+n) + thinkplot.plot(xs, window) + + start += n/2 + + # thinkplot.save(root='windowing3', + # xlabel='Index', + # axis=[0, 800, 0, 1.05]) + thinkplot.show() + +def invert_spectrogram(): + """Tests Spectrogram.make_wave. + """ + signal = thinkdsp.Chirp(start=220, end=440) + wave = signal.make_wave(duration=1, framerate=11025) + spectrogram = wave.make_spectrogram(seg_length=512) + + wave2 = spectrogram.make_wave() + + for i, (y1, y2) in enumerate(zip(wave.ys, wave2.ys)): + if abs(y1 - y2) > 1e-14: + print(i, y1, y2) + + +def main(): + chirp_spectrum() + chirp_spectrogram() + overlapping_windows() + window_plot() + three_spectrums() + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/test_code/noise.py b/test_code/noise.py new file mode 100644 index 0000000..ad78cbc --- /dev/null +++ b/test_code/noise.py @@ -0,0 +1,154 @@ +"""This file contains code used in "Think DSP", +by Allen B. Downey, available from greenteapress.com + +Copyright 2014 Allen B. Downey +License: GNU GPLv3 http://www.gnu.org/licenses/gpl.html +""" + +from __future__ import print_function, division + +import numpy +import thinkstats2 +import thinkdsp +import thinkplot + + +def process_noise(signal, root='white'): + """Plots wave and spectrum for noise signals. + + signal: Signal + root: string used to generate file names + """ + framerate = 11025 + wave = signal.make_wave(duration=0.5, framerate=framerate) + + # 0: waveform + segment = wave.segment(duration=0.1) + segment.plot(linewidth=1, alpha=0.5) + # thinkplot.save(root=root+'noise0', + # xlabel='Time (s)', + # ylim=[-1.05, 1.05]) + thinkplot.show() + + spectrum = wave.make_spectrum() + + # 1: spectrum + spectrum.plot_power(linewidth=1, alpha=0.5) + # thinkplot.save(root=root+'noise1', + # xlabel='Frequency (Hz)', + # ylabel='Power', + # xlim=[0, spectrum.fs[-1]]) + thinkplot.show() + + slope, _, _, _, _ = spectrum.estimate_slope() + print('estimated slope', slope) + + # 2: integrated spectrum + integ = spectrum.make_integrated_spectrum() + integ.plot_power() + # thinkplot.save(root=root+'noise2', + # xlabel='Frequency (Hz)', + # ylabel='Cumulative fraction of total power', + # xlim=[0, framerate/2]) + thinkplot.show() + + # 3: log-log power spectrum + spectrum.hs[0] = 0 + thinkplot.preplot(cols=2) + spectrum.plot_power(linewidth=1, alpha=0.5) + thinkplot.config(xlabel='Frequency (Hz)', + ylabel='Power', + xlim=[0, framerate/2]) + + thinkplot.subplot(2) + spectrum.plot_power(linewidth=1, alpha=0.5) + thinkplot.config(xlabel='Frequency (Hz)', + xscale='log', + yscale='log', + xlim=[0, framerate/2]) + + # thinkplot.save(root=root+'noise3') + thinkplot.show() + +def plot_gaussian_noise(): + """Shows the distribution of the spectrum of Gaussian noise. + """ + thinkdsp.random_seed(18) + signal = thinkdsp.UncorrelatedGaussianNoise() + wave = signal.make_wave(duration=0.5, framerate=11025) + spectrum = wave.make_spectrum() + + thinkplot.preplot(2, cols=2) + thinkstats2.NormalProbabilityPlot(spectrum.real, label='real') + thinkplot.config(xlabel='Normal sample', + ylabel='Amplitude', + ylim=[-250, 250], + loc='lower right') + + thinkplot.subplot(2) + thinkstats2.NormalProbabilityPlot(spectrum.imag, label='imag') + thinkplot.config(xlabel='Normal sample', + ylim=[-250, 250], + loc='lower right') + + # thinkplot.save(root='noise1') + thinkplot.show() + +def plot_pink_noise(): + """Makes a plot showing power spectrums for pink noise. + """ + thinkdsp.random_seed(20) + + duration = 1.0 + framerate = 512 + + def make_spectrum(signal): + wave = signal.make_wave(duration=duration, framerate=framerate) + spectrum = wave.make_spectrum() + spectrum.hs[0] = 0 + return spectrum + + signal = thinkdsp.UncorrelatedUniformNoise() + white = make_spectrum(signal) + + signal = thinkdsp.PinkNoise() + pink = make_spectrum(signal) + + signal = thinkdsp.BrownianNoise() + red = make_spectrum(signal) + + linewidth = 2 + # colorbrewer2.org 4-class sequential OrRd + white.plot_power(label='white', color='#fdcc8a', linewidth=linewidth) + pink.plot_power(label='pink', color='#fc8d59', linewidth=linewidth) + red.plot_power(label='red', color='#d7301f', linewidth=linewidth) + # thinkplot.save(root='noise-triple', + # xlabel='Frequency (Hz)', + # ylabel='Power', + # xscale='log', + # yscale='log', + # xlim=[1, red.fs[-1]]) + thinkplot.show() + +def main(): + thinkdsp.random_seed(17) + plot_pink_noise() + + thinkdsp.random_seed(17) + plot_gaussian_noise() + + thinkdsp.random_seed(20) + signal = thinkdsp.UncorrelatedUniformNoise() + process_noise(signal, root='white') + + thinkdsp.random_seed(20) + signal = thinkdsp.PinkNoise(beta=1.0) + process_noise(signal, root='pink') + + thinkdsp.random_seed(17) + signal = thinkdsp.BrownianNoise() + process_noise(signal, root='red') + + +if __name__ == '__main__': + main() \ No newline at end of file