Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 7307
1
"""This file contains code used in "Think DSP",
2
by Allen B. Downey, available from greenteapress.com
3
4
Copyright 2013 Allen B. Downey
5
License: MIT License (https://opensource.org/licenses/MIT)
6
"""
7
8
import copy
9
import math
10
11
import numpy as np
12
import random
13
import scipy
14
import scipy.stats
15
import scipy.fftpack
16
import subprocess
17
import warnings
18
19
from wave import open as open_wave
20
from scipy.io import wavfile
21
22
import matplotlib.pyplot as plt
23
24
try:
25
from IPython.display import Audio
26
except:
27
warnings.warn(
28
"Can't import Audio from IPython.display; " "Wave.make_audio() will not work."
29
)
30
31
PI2 = math.pi * 2
32
33
34
def random_seed(x):
35
"""Initialize the random and np.random generators.
36
37
x: int seed
38
"""
39
random.seed(x)
40
np.random.seed(x)
41
42
43
class UnimplementedMethodException(Exception):
44
"""Exception if someone calls a method that should be overridden."""
45
46
47
class WavFileWriter:
48
"""Writes wav files."""
49
50
def __init__(self, filename="sound.wav", framerate=11025):
51
"""Opens the file and sets parameters.
52
53
filename: string
54
framerate: samples per second
55
"""
56
self.filename = filename
57
self.framerate = framerate
58
self.nchannels = 1
59
self.sampwidth = 2
60
self.bits = self.sampwidth * 8
61
self.bound = 2 ** (self.bits - 1) - 1
62
63
self.fmt = "h"
64
self.dtype = np.int16
65
66
self.fp = open_wave(self.filename, "w")
67
self.fp.setnchannels(self.nchannels)
68
self.fp.setsampwidth(self.sampwidth)
69
self.fp.setframerate(self.framerate)
70
71
def write(self, wave):
72
"""Writes a wave.
73
74
wave: Wave
75
"""
76
zs = wave.quantize(self.bound, self.dtype)
77
self.fp.writeframes(zs.tostring())
78
79
def close(self, duration=0):
80
"""Closes the file.
81
82
duration: how many seconds of silence to append
83
"""
84
if duration:
85
self.write(rest(duration))
86
87
self.fp.close()
88
89
90
def read_wave(filename="sound.wav"):
91
"""Reads a wave file.
92
93
filename: string
94
95
returns: Wave
96
"""
97
fp = open_wave(filename, "r")
98
99
nchannels = fp.getnchannels()
100
nframes = fp.getnframes()
101
sampwidth = fp.getsampwidth()
102
framerate = fp.getframerate()
103
104
z_str = fp.readframes(nframes)
105
106
fp.close()
107
108
dtype_map = {1: np.int8, 2: np.int16, 3: "special", 4: np.int32}
109
if sampwidth not in dtype_map:
110
raise ValueError("sampwidth %d unknown" % sampwidth)
111
112
if sampwidth == 3:
113
xs = np.fromstring(z_str, dtype=np.int8).astype(np.int32)
114
ys = (xs[2::3] * 256 + xs[1::3]) * 256 + xs[0::3]
115
else:
116
ys = np.fromstring(z_str, dtype=dtype_map[sampwidth])
117
118
# if it's in stereo, just pull out the first channel
119
if nchannels == 2:
120
ys = ys[::2]
121
122
# ts = np.arange(len(ys)) / framerate
123
wave = Wave(ys, framerate=framerate)
124
wave.normalize()
125
return wave
126
127
128
def read_wave_with_scipy(filename):
129
"""Reads a wave file.
130
131
filename: string
132
133
returns: Wave
134
"""
135
# TODO: Check back later and see if this works on 24-bit data,
136
# and runs without throwing warnings.
137
framerate, ys = wavfile.read(filename)
138
139
# if it's in stereo, just pull out the first channel
140
if ys.ndim == 2:
141
ys = ys[:, 0]
142
143
# ts = np.arange(len(ys)) / framerate
144
wave = Wave(ys, framerate=framerate)
145
wave.normalize()
146
return wave
147
148
149
def play_wave(filename="sound.wav", player="aplay"):
150
"""Plays a wave file.
151
152
filename: string
153
player: string name of executable that plays wav files
154
"""
155
cmd = "%s %s" % (player, filename)
156
popen = subprocess.Popen(cmd, shell=True)
157
popen.communicate()
158
159
160
def find_index(x, xs):
161
"""Find the index corresponding to a given value in an array."""
162
n = len(xs)
163
start = xs[0]
164
end = xs[-1]
165
i = round((n - 1) * (x - start) / (end - start))
166
return int(i)
167
168
169
class _SpectrumParent:
170
"""Contains code common to Spectrum and DCT.
171
"""
172
173
def __init__(self, hs, fs, framerate, full=False):
174
"""Initializes a spectrum.
175
176
hs: array of amplitudes (real or complex)
177
fs: array of frequencies
178
framerate: frames per second
179
full: boolean to indicate full or real FFT
180
"""
181
self.hs = np.asanyarray(hs)
182
self.fs = np.asanyarray(fs)
183
self.framerate = framerate
184
self.full = full
185
186
@property
187
def max_freq(self):
188
"""Returns the Nyquist frequency for this spectrum."""
189
return self.framerate / 2
190
191
@property
192
def amps(self):
193
"""Returns a sequence of amplitudes (read-only property)."""
194
return np.absolute(self.hs)
195
196
@property
197
def power(self):
198
"""Returns a sequence of powers (read-only property)."""
199
return self.amps ** 2
200
201
def copy(self):
202
"""Makes a copy.
203
204
Returns: new Spectrum
205
"""
206
return copy.deepcopy(self)
207
208
def max_diff(self, other):
209
"""Computes the maximum absolute difference between spectra.
210
211
other: Spectrum
212
213
returns: float
214
"""
215
assert self.framerate == other.framerate
216
assert len(self) == len(other)
217
218
hs = self.hs - other.hs
219
return np.max(np.abs(hs))
220
221
def ratio(self, denom, thresh=1, val=0):
222
"""The ratio of two spectrums.
223
224
denom: Spectrum
225
thresh: values smaller than this are replaced
226
val: with this value
227
228
returns: new Wave
229
"""
230
ratio_spectrum = self.copy()
231
ratio_spectrum.hs /= denom.hs
232
ratio_spectrum.hs[denom.amps < thresh] = val
233
return ratio_spectrum
234
235
def invert(self):
236
"""Inverts this spectrum/filter.
237
238
returns: new Wave
239
"""
240
inverse = self.copy()
241
inverse.hs = 1 / inverse.hs
242
return inverse
243
244
@property
245
def freq_res(self):
246
return self.framerate / 2 / (len(self.fs) - 1)
247
248
def render_full(self, high=None):
249
"""Extracts amps and fs from a full spectrum.
250
251
high: cutoff frequency
252
253
returns: fs, amps
254
"""
255
hs = np.fft.fftshift(self.hs)
256
amps = np.abs(hs)
257
fs = np.fft.fftshift(self.fs)
258
i = 0 if high is None else find_index(-high, fs)
259
j = None if high is None else find_index(high, fs) + 1
260
return fs[i:j], amps[i:j]
261
262
def plot(self, high=None, **options):
263
"""Plots amplitude vs frequency.
264
265
Note: if this is a full spectrum, it ignores low and high
266
267
high: frequency to cut off at
268
"""
269
if self.full:
270
fs, amps = self.render_full(high)
271
plt.plot(fs, amps, **options)
272
else:
273
i = None if high is None else find_index(high, self.fs)
274
plt.plot(self.fs[:i], self.amps[:i], **options)
275
276
def plot_power(self, high=None, **options):
277
"""Plots power vs frequency.
278
279
high: frequency to cut off at
280
"""
281
if self.full:
282
fs, amps = self.render_full(high)
283
plt.plot(fs, amps ** 2, **options)
284
else:
285
i = None if high is None else find_index(high, self.fs)
286
plt.plot(self.fs[:i], self.power[:i], **options)
287
288
def estimate_slope(self):
289
"""Runs linear regression on log power vs log frequency.
290
291
returns: slope, inter, r2, p, stderr
292
"""
293
x = np.log(self.fs[1:])
294
y = np.log(self.power[1:])
295
t = scipy.stats.linregress(x, y)
296
return t
297
298
def peaks(self):
299
"""Finds the highest peaks and their frequencies.
300
301
returns: sorted list of (amplitude, frequency) pairs
302
"""
303
t = list(zip(self.amps, self.fs))
304
t.sort(reverse=True)
305
return t
306
307
308
class Spectrum(_SpectrumParent):
309
"""Represents the spectrum of a signal."""
310
311
def __len__(self):
312
"""Length of the spectrum."""
313
return len(self.hs)
314
315
def __add__(self, other):
316
"""Adds two spectrums elementwise.
317
318
other: Spectrum
319
320
returns: new Spectrum
321
"""
322
if other == 0:
323
return self.copy()
324
325
assert all(self.fs == other.fs)
326
hs = self.hs + other.hs
327
return Spectrum(hs, self.fs, self.framerate, self.full)
328
329
__radd__ = __add__
330
331
def __mul__(self, other):
332
"""Multiplies two spectrums elementwise.
333
334
other: Spectrum
335
336
returns: new Spectrum
337
"""
338
assert all(self.fs == other.fs)
339
hs = self.hs * other.hs
340
return Spectrum(hs, self.fs, self.framerate, self.full)
341
342
def convolve(self, other):
343
"""Convolves two Spectrums.
344
345
other: Spectrum
346
347
returns: Spectrum
348
"""
349
assert all(self.fs == other.fs)
350
if self.full:
351
hs1 = np.fft.fftshift(self.hs)
352
hs2 = np.fft.fftshift(other.hs)
353
hs = np.convolve(hs1, hs2, mode="same")
354
hs = np.fft.ifftshift(hs)
355
else:
356
# not sure this branch would mean very much
357
hs = np.convolve(self.hs, other.hs, mode="same")
358
359
return Spectrum(hs, self.fs, self.framerate, self.full)
360
361
@property
362
def real(self):
363
"""Returns the real part of the hs (read-only property)."""
364
return np.real(self.hs)
365
366
@property
367
def imag(self):
368
"""Returns the imaginary part of the hs (read-only property)."""
369
return np.imag(self.hs)
370
371
@property
372
def angles(self):
373
"""Returns a sequence of angles (read-only property)."""
374
return np.angle(self.hs)
375
376
def scale(self, factor):
377
"""Multiplies all elements by the given factor.
378
379
factor: what to multiply the magnitude by (could be complex)
380
"""
381
self.hs *= factor
382
383
def low_pass(self, cutoff, factor=0):
384
"""Attenuate frequencies above the cutoff.
385
386
cutoff: frequency in Hz
387
factor: what to multiply the magnitude by
388
"""
389
self.hs[abs(self.fs) > cutoff] *= factor
390
391
def high_pass(self, cutoff, factor=0):
392
"""Attenuate frequencies below the cutoff.
393
394
cutoff: frequency in Hz
395
factor: what to multiply the magnitude by
396
"""
397
self.hs[abs(self.fs) < cutoff] *= factor
398
399
def band_stop(self, low_cutoff, high_cutoff, factor=0):
400
"""Attenuate frequencies between the cutoffs.
401
402
low_cutoff: frequency in Hz
403
high_cutoff: frequency in Hz
404
factor: what to multiply the magnitude by
405
"""
406
# TODO: test this function
407
fs = abs(self.fs)
408
indices = (low_cutoff < fs) & (fs < high_cutoff)
409
self.hs[indices] *= factor
410
411
def pink_filter(self, beta=1):
412
"""Apply a filter that would make white noise pink.
413
414
beta: exponent of the pink noise
415
"""
416
denom = self.fs ** (beta / 2.0)
417
denom[0] = 1
418
self.hs /= denom
419
420
def differentiate(self):
421
"""Apply the differentiation filter.
422
423
returns: new Spectrum
424
"""
425
new = self.copy()
426
new.hs *= PI2 * 1j * new.fs
427
return new
428
429
def integrate(self):
430
"""Apply the integration filter.
431
432
returns: new Spectrum
433
"""
434
new = self.copy()
435
zero = (new.fs == 0)
436
new.hs[~zero] /= PI2 * 1j * new.fs[~zero]
437
new.hs[zero] = np.inf
438
return new
439
440
def make_integrated_spectrum(self):
441
"""Makes an integrated spectrum.
442
"""
443
cs = np.cumsum(self.power)
444
cs /= cs[-1]
445
return IntegratedSpectrum(cs, self.fs)
446
447
def make_wave(self):
448
"""Transforms to the time domain.
449
450
returns: Wave
451
"""
452
if self.full:
453
ys = np.fft.ifft(self.hs)
454
else:
455
ys = np.fft.irfft(self.hs)
456
457
# NOTE: whatever the start time was, we lose it when
458
# we transform back; we could fix that by saving start
459
# time in the Spectrum
460
# ts = self.start + np.arange(len(ys)) / self.framerate
461
return Wave(ys, framerate=self.framerate)
462
463
464
class IntegratedSpectrum:
465
"""Represents the integral of a spectrum."""
466
467
def __init__(self, cs, fs):
468
"""Initializes an integrated spectrum:
469
470
cs: sequence of cumulative amplitudes
471
fs: sequence of frequencies
472
"""
473
self.cs = np.asanyarray(cs)
474
self.fs = np.asanyarray(fs)
475
476
def plot_power(self, low=0, high=None, expo=False, **options):
477
"""Plots the integrated spectrum.
478
479
low: int index to start at
480
high: int index to end at
481
"""
482
cs = self.cs[low:high]
483
fs = self.fs[low:high]
484
485
if expo:
486
cs = np.exp(cs)
487
488
plt.plot(fs, cs, **options)
489
490
def estimate_slope(self, low=1, high=-12000):
491
"""Runs linear regression on log cumulative power vs log frequency.
492
493
returns: slope, inter, r2, p, stderr
494
"""
495
# print self.fs[low:high]
496
# print self.cs[low:high]
497
x = np.log(self.fs[low:high])
498
y = np.log(self.cs[low:high])
499
t = scipy.stats.linregress(x, y)
500
return t
501
502
503
class Dct(_SpectrumParent):
504
"""Represents the spectrum of a signal using discrete cosine transform."""
505
506
@property
507
def amps(self):
508
"""Returns a sequence of amplitudes (read-only property).
509
510
Note: for DCTs, amps are positive or negative real.
511
"""
512
return self.hs
513
514
def __add__(self, other):
515
"""Adds two DCTs elementwise.
516
517
other: DCT
518
519
returns: new DCT
520
"""
521
if other == 0:
522
return self
523
524
assert self.framerate == other.framerate
525
hs = self.hs + other.hs
526
return Dct(hs, self.fs, self.framerate)
527
528
__radd__ = __add__
529
530
def make_wave(self):
531
"""Transforms to the time domain.
532
533
returns: Wave
534
"""
535
N = len(self.hs)
536
ys = scipy.fftpack.idct(self.hs, type=2) / 2 / N
537
# NOTE: whatever the start time was, we lose it when
538
# we transform back
539
# ts = self.start + np.arange(len(ys)) / self.framerate
540
return Wave(ys, framerate=self.framerate)
541
542
543
class Spectrogram:
544
"""Represents the spectrum of a signal."""
545
546
def __init__(self, spec_map, seg_length):
547
"""Initialize the spectrogram.
548
549
spec_map: map from float time to Spectrum
550
seg_length: number of samples in each segment
551
"""
552
self.spec_map = spec_map
553
self.seg_length = seg_length
554
555
def any_spectrum(self):
556
"""Returns an arbitrary spectrum from the spectrogram."""
557
index = next(iter(self.spec_map))
558
return self.spec_map[index]
559
560
@property
561
def time_res(self):
562
"""Time resolution in seconds."""
563
spectrum = self.any_spectrum()
564
return float(self.seg_length) / spectrum.framerate
565
566
@property
567
def freq_res(self):
568
"""Frequency resolution in Hz."""
569
return self.any_spectrum().freq_res
570
571
def times(self):
572
"""Sorted sequence of times.
573
574
returns: sequence of float times in seconds
575
"""
576
ts = sorted(iter(self.spec_map))
577
return ts
578
579
def frequencies(self):
580
"""Sequence of frequencies.
581
582
returns: sequence of float freqencies in Hz.
583
"""
584
fs = self.any_spectrum().fs
585
return fs
586
587
def plot(self, high=None, **options):
588
"""Make a pseudocolor plot.
589
590
high: highest frequency component to plot
591
"""
592
fs = self.frequencies()
593
i = None if high is None else find_index(high, fs)
594
fs = fs[:i]
595
ts = self.times()
596
597
# make the array
598
size = len(fs), len(ts)
599
array = np.zeros(size, dtype=np.float)
600
601
# copy amplitude from each spectrum into a column of the array
602
for j, t in enumerate(ts):
603
spectrum = self.spec_map[t]
604
array[:, j] = spectrum.amps[:i]
605
606
underride(options, cmap='inferno_r', shading='auto')
607
plt.pcolormesh(ts, fs, array, **options)
608
609
def get_data(self, high=None, **options):
610
"""Returns spectogram as 2D numpy array
611
612
high: highest frequency component to return
613
"""
614
fs = self.frequencies()
615
i = None if high is None else find_index(high, fs)
616
fs = fs[:i]
617
ts = self.times()
618
619
# make the array
620
size = len(fs), len(ts)
621
array = np.zeros(size, dtype=np.float)
622
623
# copy amplitude from each spectrum into a column of the array
624
for j, t in enumerate(ts):
625
spectrum = self.spec_map[t]
626
array[:, j] = spectrum.amps[:i]
627
628
return array
629
630
def make_wave(self):
631
"""Inverts the spectrogram and returns a Wave.
632
633
returns: Wave
634
"""
635
res = []
636
for t, spectrum in sorted(self.spec_map.items()):
637
wave = spectrum.make_wave()
638
n = len(wave)
639
640
window = 1 / np.hamming(n)
641
wave.window(window)
642
643
i = wave.find_index(t)
644
start = i - n // 2
645
end = start + n
646
res.append((start, end, wave))
647
648
starts, ends, waves = zip(*res)
649
low = min(starts)
650
high = max(ends)
651
652
ys = np.zeros(high - low, np.float)
653
for start, end, wave in res:
654
ys[start:end] = wave.ys
655
656
# ts = np.arange(len(ys)) / self.framerate
657
return Wave(ys, framerate=wave.framerate)
658
659
660
class Wave:
661
"""Represents a discrete-time waveform.
662
663
"""
664
665
def __init__(self, ys, ts=None, framerate=None):
666
"""Initializes the wave.
667
668
ys: wave array
669
ts: array of times
670
framerate: samples per second
671
"""
672
self.ys = np.asanyarray(ys)
673
self.framerate = framerate if framerate is not None else 11025
674
675
if ts is None:
676
self.ts = np.arange(len(ys)) / self.framerate
677
else:
678
self.ts = np.asanyarray(ts)
679
680
def copy(self):
681
"""Makes a copy.
682
683
Returns: new Wave
684
"""
685
return copy.deepcopy(self)
686
687
def __len__(self):
688
return len(self.ys)
689
690
@property
691
def start(self):
692
return self.ts[0]
693
694
@property
695
def end(self):
696
return self.ts[-1]
697
698
@property
699
def duration(self):
700
"""Duration (property).
701
702
returns: float duration in seconds
703
"""
704
return len(self.ys) / self.framerate
705
706
def __add__(self, other):
707
"""Adds two waves elementwise.
708
709
other: Wave
710
711
returns: new Wave
712
"""
713
if other == 0:
714
return self
715
716
assert self.framerate == other.framerate
717
718
# make an array of times that covers both waves
719
start = min(self.start, other.start)
720
end = max(self.end, other.end)
721
n = int(round((end - start) * self.framerate)) + 1
722
ys = np.zeros(n)
723
ts = start + np.arange(n) / self.framerate
724
725
def add_ys(wave):
726
i = find_index(wave.start, ts)
727
728
# make sure the arrays line up reasonably well
729
diff = ts[i] - wave.start
730
dt = 1 / wave.framerate
731
if (diff / dt) > 0.1:
732
warnings.warn(
733
"Can't add these waveforms; their " "time arrays don't line up."
734
)
735
736
j = i + len(wave)
737
ys[i:j] += wave.ys
738
739
add_ys(self)
740
add_ys(other)
741
742
return Wave(ys, ts, self.framerate)
743
744
__radd__ = __add__
745
746
def __or__(self, other):
747
"""Concatenates two waves.
748
749
other: Wave
750
751
returns: new Wave
752
"""
753
if self.framerate != other.framerate:
754
raise ValueError("Wave.__or__: framerates do not agree")
755
756
ys = np.concatenate((self.ys, other.ys))
757
# ts = np.arange(len(ys)) / self.framerate
758
return Wave(ys, framerate=self.framerate)
759
760
def __mul__(self, other):
761
"""Multiplies two waves elementwise.
762
763
Note: this operation ignores the timestamps; the result
764
has the timestamps of self.
765
766
other: Wave
767
768
returns: new Wave
769
"""
770
# the spectrums have to have the same framerate and duration
771
assert self.framerate == other.framerate
772
assert len(self) == len(other)
773
774
ys = self.ys * other.ys
775
return Wave(ys, self.ts, self.framerate)
776
777
def max_diff(self, other):
778
"""Computes the maximum absolute difference between waves.
779
780
other: Wave
781
782
returns: float
783
"""
784
assert self.framerate == other.framerate
785
assert len(self) == len(other)
786
787
ys = self.ys - other.ys
788
return np.max(np.abs(ys))
789
790
def convolve(self, other):
791
"""Convolves two waves.
792
793
Note: this operation ignores the timestamps; the result
794
has the timestamps of self.
795
796
other: Wave or NumPy array
797
798
returns: Wave
799
"""
800
if isinstance(other, Wave):
801
assert self.framerate == other.framerate
802
window = other.ys
803
else:
804
window = other
805
806
ys = np.convolve(self.ys, window, mode="full")
807
# ts = np.arange(len(ys)) / self.framerate
808
return Wave(ys, framerate=self.framerate)
809
810
def diff(self):
811
"""Computes the difference between successive elements.
812
813
returns: new Wave
814
"""
815
ys = np.diff(self.ys)
816
ts = self.ts[1:].copy()
817
return Wave(ys, ts, self.framerate)
818
819
def cumsum(self):
820
"""Computes the cumulative sum of the elements.
821
822
returns: new Wave
823
"""
824
ys = np.cumsum(self.ys)
825
ts = self.ts.copy()
826
return Wave(ys, ts, self.framerate)
827
828
def quantize(self, bound, dtype):
829
"""Maps the waveform to quanta.
830
831
bound: maximum amplitude
832
dtype: numpy data type or string
833
834
returns: quantized signal
835
"""
836
return quantize(self.ys, bound, dtype)
837
838
def apodize(self, denom=20, duration=0.1):
839
"""Tapers the amplitude at the beginning and end of the signal.
840
841
Tapers either the given duration of time or the given
842
fraction of the total duration, whichever is less.
843
844
denom: float fraction of the segment to taper
845
duration: float duration of the taper in seconds
846
"""
847
self.ys = apodize(self.ys, self.framerate, denom, duration)
848
849
def hamming(self):
850
"""Apply a Hamming window to the wave.
851
"""
852
self.ys *= np.hamming(len(self.ys))
853
854
def window(self, window):
855
"""Apply a window to the wave.
856
857
window: sequence of multipliers, same length as self.ys
858
"""
859
self.ys *= window
860
861
def scale(self, factor):
862
"""Multplies the wave by a factor.
863
864
factor: scale factor
865
"""
866
self.ys *= factor
867
868
def shift(self, shift):
869
"""Shifts the wave left or right in time.
870
871
shift: float time shift
872
"""
873
# TODO: track down other uses of this function and check them
874
self.ts += shift
875
876
def roll(self, roll):
877
"""Rolls this wave by the given number of locations.
878
"""
879
self.ys = np.roll(self.ys, roll)
880
881
def truncate(self, n):
882
"""Trims this wave to the given length.
883
884
n: integer index
885
"""
886
self.ys = truncate(self.ys, n)
887
self.ts = truncate(self.ts, n)
888
889
def zero_pad(self, n):
890
"""Trims this wave to the given length.
891
892
n: integer index
893
"""
894
self.ys = zero_pad(self.ys, n)
895
self.ts = self.start + np.arange(n) / self.framerate
896
897
def normalize(self, amp=1.0):
898
"""Normalizes the signal to the given amplitude.
899
900
amp: float amplitude
901
"""
902
self.ys = normalize(self.ys, amp=amp)
903
904
def unbias(self):
905
"""Unbiases the signal.
906
"""
907
self.ys = unbias(self.ys)
908
909
def find_index(self, t):
910
"""Find the index corresponding to a given time."""
911
n = len(self)
912
start = self.start
913
end = self.end
914
i = round((n - 1) * (t - start) / (end - start))
915
return int(i)
916
917
def segment(self, start=None, duration=None):
918
"""Extracts a segment.
919
920
start: float start time in seconds
921
duration: float duration in seconds
922
923
returns: Wave
924
"""
925
if start is None:
926
start = self.ts[0]
927
i = 0
928
else:
929
i = self.find_index(start)
930
931
j = None if duration is None else self.find_index(start + duration)
932
return self.slice(i, j)
933
934
def slice(self, i, j):
935
"""Makes a slice from a Wave.
936
937
i: first slice index
938
j: second slice index
939
"""
940
ys = self.ys[i:j].copy()
941
ts = self.ts[i:j].copy()
942
return Wave(ys, ts, self.framerate)
943
944
def make_spectrum(self, full=False):
945
"""Computes the spectrum using FFT.
946
947
full: boolean, whethere to compute a full FFT
948
(as opposed to a real FFT)
949
950
returns: Spectrum
951
"""
952
n = len(self.ys)
953
d = 1 / self.framerate
954
955
if full:
956
hs = np.fft.fft(self.ys)
957
fs = np.fft.fftfreq(n, d)
958
else:
959
hs = np.fft.rfft(self.ys)
960
fs = np.fft.rfftfreq(n, d)
961
962
return Spectrum(hs, fs, self.framerate, full)
963
964
def make_dct(self):
965
"""Computes the DCT of this wave.
966
"""
967
N = len(self.ys)
968
hs = scipy.fftpack.dct(self.ys, type=2)
969
fs = (0.5 + np.arange(N)) / 2
970
return Dct(hs, fs, self.framerate)
971
972
def make_spectrogram(self, seg_length, win_flag=True):
973
"""Computes the spectrogram of the wave.
974
975
seg_length: number of samples in each segment
976
win_flag: boolean, whether to apply hamming window to each segment
977
978
returns: Spectrogram
979
"""
980
if win_flag:
981
window = np.hamming(seg_length)
982
i, j = 0, seg_length
983
step = int(seg_length // 2)
984
985
# map from time to Spectrum
986
spec_map = {}
987
988
while j < len(self.ys):
989
segment = self.slice(i, j)
990
if win_flag:
991
segment.window(window)
992
993
# the nominal time for this segment is the midpoint
994
t = (segment.start + segment.end) / 2
995
spec_map[t] = segment.make_spectrum()
996
997
i += step
998
j += step
999
1000
return Spectrogram(spec_map, seg_length)
1001
1002
def get_xfactor(self, options):
1003
try:
1004
xfactor = options["xfactor"]
1005
options.pop("xfactor")
1006
except KeyError:
1007
xfactor = 1
1008
return xfactor
1009
1010
def plot(self, **options):
1011
"""Plots the wave.
1012
1013
If the ys are complex, plots the real part.
1014
1015
"""
1016
xfactor = self.get_xfactor(options)
1017
plt.plot(self.ts * xfactor, np.real(self.ys), **options)
1018
1019
def plot_vlines(self, **options):
1020
"""Plots the wave with vertical lines for samples.
1021
1022
"""
1023
xfactor = self.get_xfactor(options)
1024
plt.vlines(self.ts * xfactor, 0, self.ys, **options)
1025
1026
def corr(self, other):
1027
"""Correlation coefficient two waves.
1028
1029
other: Wave
1030
1031
returns: float coefficient of correlation
1032
"""
1033
corr = np.corrcoef(self.ys, other.ys)[0, 1]
1034
return corr
1035
1036
def cov_mat(self, other):
1037
"""Covariance matrix of two waves.
1038
1039
other: Wave
1040
1041
returns: 2x2 covariance matrix
1042
"""
1043
return np.cov(self.ys, other.ys)
1044
1045
def cov(self, other):
1046
"""Covariance of two unbiased waves.
1047
1048
other: Wave
1049
1050
returns: float
1051
"""
1052
total = sum(self.ys * other.ys) / len(self.ys)
1053
return total
1054
1055
def cos_cov(self, k):
1056
"""Covariance with a cosine signal.
1057
1058
freq: freq of the cosine signal in Hz
1059
1060
returns: float covariance
1061
"""
1062
n = len(self.ys)
1063
factor = math.pi * k / n
1064
ys = [math.cos(factor * (i + 0.5)) for i in range(n)]
1065
total = 2 * sum(self.ys * ys)
1066
return total
1067
1068
def cos_transform(self):
1069
"""Discrete cosine transform.
1070
1071
returns: list of frequency, cov pairs
1072
"""
1073
n = len(self.ys)
1074
res = []
1075
for k in range(n):
1076
cov = self.cos_cov(k)
1077
res.append((k, cov))
1078
1079
return res
1080
1081
def write(self, filename="sound.wav"):
1082
"""Write a wave file.
1083
1084
filename: string
1085
"""
1086
print("Writing", filename)
1087
wfile = WavFileWriter(filename, self.framerate)
1088
wfile.write(self)
1089
wfile.close()
1090
1091
def play(self, filename="sound.wav"):
1092
"""Plays a wave file.
1093
1094
filename: string
1095
"""
1096
self.write(filename)
1097
play_wave(filename)
1098
1099
def make_audio(self):
1100
"""Makes an IPython Audio object.
1101
"""
1102
audio = Audio(data=self.ys.real, rate=self.framerate)
1103
return audio
1104
1105
1106
def unbias(ys):
1107
"""Shifts a wave array so it has mean 0.
1108
1109
ys: wave array
1110
1111
returns: wave array
1112
"""
1113
return ys - ys.mean()
1114
1115
1116
def normalize(ys, amp=1.0):
1117
"""Normalizes a wave array so the maximum amplitude is +amp or -amp.
1118
1119
ys: wave array
1120
amp: max amplitude (pos or neg) in result
1121
1122
returns: wave array
1123
"""
1124
high, low = abs(max(ys)), abs(min(ys))
1125
return amp * ys / max(high, low)
1126
1127
1128
def shift_right(ys, shift):
1129
"""Shifts a wave array to the right and zero pads.
1130
1131
ys: wave array
1132
shift: integer shift
1133
1134
returns: wave array
1135
"""
1136
res = np.zeros(len(ys) + shift)
1137
res[shift:] = ys
1138
return res
1139
1140
1141
def shift_left(ys, shift):
1142
"""Shifts a wave array to the left.
1143
1144
ys: wave array
1145
shift: integer shift
1146
1147
returns: wave array
1148
"""
1149
return ys[shift:]
1150
1151
1152
def truncate(ys, n):
1153
"""Trims a wave array to the given length.
1154
1155
ys: wave array
1156
n: integer length
1157
1158
returns: wave array
1159
"""
1160
return ys[:n]
1161
1162
1163
def quantize(ys, bound, dtype):
1164
"""Maps the waveform to quanta.
1165
1166
ys: wave array
1167
bound: maximum amplitude
1168
dtype: numpy data type of the result
1169
1170
returns: quantized signal
1171
"""
1172
if max(ys) > 1 or min(ys) < -1:
1173
warnings.warn("Warning: normalizing before quantizing.")
1174
ys = normalize(ys)
1175
1176
zs = (ys * bound).astype(dtype)
1177
return zs
1178
1179
1180
def apodize(ys, framerate, denom=20, duration=0.1):
1181
"""Tapers the amplitude at the beginning and end of the signal.
1182
1183
Tapers either the given duration of time or the given
1184
fraction of the total duration, whichever is less.
1185
1186
ys: wave array
1187
framerate: int frames per second
1188
denom: float fraction of the segment to taper
1189
duration: float duration of the taper in seconds
1190
1191
returns: wave array
1192
"""
1193
# a fixed fraction of the segment
1194
n = len(ys)
1195
k1 = n // denom
1196
1197
# a fixed duration of time
1198
k2 = int(duration * framerate)
1199
1200
k = min(k1, k2)
1201
1202
w1 = np.linspace(0, 1, k)
1203
w2 = np.ones(n - 2 * k)
1204
w3 = np.linspace(1, 0, k)
1205
1206
window = np.concatenate((w1, w2, w3))
1207
return ys * window
1208
1209
1210
class Signal:
1211
"""Represents a time-varying signal."""
1212
1213
def __add__(self, other):
1214
"""Adds two signals.
1215
1216
other: Signal
1217
1218
returns: Signal
1219
"""
1220
if other == 0:
1221
return self
1222
return SumSignal(self, other)
1223
1224
__radd__ = __add__
1225
1226
@property
1227
def period(self):
1228
"""Period of the signal in seconds (property).
1229
1230
Since this is used primarily for purposes of plotting,
1231
the default behavior is to return a value, 0.1 seconds,
1232
that is reasonable for many signals.
1233
1234
returns: float seconds
1235
"""
1236
return 0.1
1237
1238
def plot(self, framerate=11025):
1239
"""Plots the signal.
1240
1241
The default behavior is to plot three periods.
1242
1243
framerate: samples per second
1244
"""
1245
duration = self.period * 3
1246
wave = self.make_wave(duration, start=0, framerate=framerate)
1247
wave.plot()
1248
1249
def make_wave(self, duration=1, start=0, framerate=11025):
1250
"""Makes a Wave object.
1251
1252
duration: float seconds
1253
start: float seconds
1254
framerate: int frames per second
1255
1256
returns: Wave
1257
"""
1258
n = round(duration * framerate)
1259
ts = start + np.arange(n) / framerate
1260
ys = self.evaluate(ts)
1261
return Wave(ys, ts, framerate=framerate)
1262
1263
1264
def infer_framerate(ts):
1265
"""Given ts, find the framerate.
1266
1267
Assumes that the ts are equally spaced.
1268
1269
ts: sequence of times in seconds
1270
1271
returns: frames per second
1272
"""
1273
# TODO: confirm that this is never used and remove it
1274
dt = ts[1] - ts[0]
1275
framerate = 1.0 / dt
1276
return framerate
1277
1278
1279
class SumSignal(Signal):
1280
"""Represents the sum of signals."""
1281
1282
def __init__(self, *args):
1283
"""Initializes the sum.
1284
1285
args: tuple of signals
1286
"""
1287
self.signals = args
1288
1289
@property
1290
def period(self):
1291
"""Period of the signal in seconds.
1292
1293
Note: this is not correct; it's mostly a placekeeper.
1294
1295
But it is correct for a harmonic sequence where all
1296
component frequencies are multiples of the fundamental.
1297
1298
returns: float seconds
1299
"""
1300
return max(sig.period for sig in self.signals)
1301
1302
def evaluate(self, ts):
1303
"""Evaluates the signal at the given times.
1304
1305
ts: float array of times
1306
1307
returns: float wave array
1308
"""
1309
ts = np.asarray(ts)
1310
return sum(sig.evaluate(ts) for sig in self.signals)
1311
1312
1313
class Sinusoid(Signal):
1314
"""Represents a sinusoidal signal."""
1315
1316
def __init__(self, freq=440, amp=1.0, offset=0, func=np.sin):
1317
"""Initializes a sinusoidal signal.
1318
1319
freq: float frequency in Hz
1320
amp: float amplitude, 1.0 is nominal max
1321
offset: float phase offset in radians
1322
func: function that maps phase to amplitude
1323
"""
1324
self.freq = freq
1325
self.amp = amp
1326
self.offset = offset
1327
self.func = func
1328
1329
@property
1330
def period(self):
1331
"""Period of the signal in seconds.
1332
1333
returns: float seconds
1334
"""
1335
return 1.0 / self.freq
1336
1337
def evaluate(self, ts):
1338
"""Evaluates the signal at the given times.
1339
1340
ts: float array of times
1341
1342
returns: float wave array
1343
"""
1344
ts = np.asarray(ts)
1345
phases = PI2 * self.freq * ts + self.offset
1346
ys = self.amp * self.func(phases)
1347
return ys
1348
1349
1350
def CosSignal(freq=440, amp=1.0, offset=0):
1351
"""Makes a cosine Sinusoid.
1352
1353
freq: float frequency in Hz
1354
amp: float amplitude, 1.0 is nominal max
1355
offset: float phase offset in radians
1356
1357
returns: Sinusoid object
1358
"""
1359
return Sinusoid(freq, amp, offset, func=np.cos)
1360
1361
1362
def SinSignal(freq=440, amp=1.0, offset=0):
1363
"""Makes a sine Sinusoid.
1364
1365
freq: float frequency in Hz
1366
amp: float amplitude, 1.0 is nominal max
1367
offset: float phase offset in radians
1368
1369
returns: Sinusoid object
1370
"""
1371
return Sinusoid(freq, amp, offset, func=np.sin)
1372
1373
1374
def Sinc(freq=440, amp=1.0, offset=0):
1375
"""Makes a Sinc function.
1376
1377
freq: float frequency in Hz
1378
amp: float amplitude, 1.0 is nominal max
1379
offset: float phase offset in radians
1380
1381
returns: Sinusoid object
1382
"""
1383
return Sinusoid(freq, amp, offset, func=np.sinc)
1384
1385
1386
class ComplexSinusoid(Sinusoid):
1387
"""Represents a complex exponential signal."""
1388
1389
def evaluate(self, ts):
1390
"""Evaluates the signal at the given times.
1391
1392
ts: float array of times
1393
1394
returns: float wave array
1395
"""
1396
ts = np.asarray(ts)
1397
phases = PI2 * self.freq * ts + self.offset
1398
ys = self.amp * np.exp(1j * phases)
1399
return ys
1400
1401
1402
class SquareSignal(Sinusoid):
1403
"""Represents a square signal."""
1404
1405
def evaluate(self, ts):
1406
"""Evaluates the signal at the given times.
1407
1408
ts: float array of times
1409
1410
returns: float wave array
1411
"""
1412
ts = np.asarray(ts)
1413
cycles = self.freq * ts + self.offset / PI2
1414
frac, _ = np.modf(cycles)
1415
ys = self.amp * np.sign(unbias(frac))
1416
return ys
1417
1418
1419
class SawtoothSignal(Sinusoid):
1420
"""Represents a sawtooth signal."""
1421
1422
def evaluate(self, ts):
1423
"""Evaluates the signal at the given times.
1424
1425
ts: float array of times
1426
1427
returns: float wave array
1428
"""
1429
ts = np.asarray(ts)
1430
cycles = self.freq * ts + self.offset / PI2
1431
frac, _ = np.modf(cycles)
1432
ys = normalize(unbias(frac), self.amp)
1433
return ys
1434
1435
1436
class ParabolicSignal(Sinusoid):
1437
"""Represents a parabolic signal."""
1438
1439
def evaluate(self, ts):
1440
"""Evaluates the signal at the given times.
1441
1442
ts: float array of times
1443
1444
returns: float wave array
1445
"""
1446
ts = np.asarray(ts)
1447
cycles = self.freq * ts + self.offset / PI2
1448
frac, _ = np.modf(cycles)
1449
ys = (frac - 0.5) ** 2
1450
ys = normalize(unbias(ys), self.amp)
1451
return ys
1452
1453
1454
class CubicSignal(ParabolicSignal):
1455
"""Represents a cubic signal."""
1456
1457
def evaluate(self, ts):
1458
"""Evaluates the signal at the given times.
1459
1460
ts: float array of times
1461
1462
returns: float wave array
1463
"""
1464
ys = ParabolicSignal.evaluate(self, ts)
1465
ys = np.cumsum(ys)
1466
ys = normalize(unbias(ys), self.amp)
1467
return ys
1468
1469
1470
class GlottalSignal(Sinusoid):
1471
"""Represents a periodic signal that resembles a glottal signal."""
1472
1473
def evaluate(self, ts):
1474
"""Evaluates the signal at the given times.
1475
1476
ts: float array of times
1477
1478
returns: float wave array
1479
"""
1480
ts = np.asarray(ts)
1481
cycles = self.freq * ts + self.offset / PI2
1482
frac, _ = np.modf(cycles)
1483
ys = frac ** 2 * (1 - frac)
1484
ys = normalize(unbias(ys), self.amp)
1485
return ys
1486
1487
1488
class TriangleSignal(Sinusoid):
1489
"""Represents a triangle signal."""
1490
1491
def evaluate(self, ts):
1492
"""Evaluates the signal at the given times.
1493
1494
ts: float array of times
1495
1496
returns: float wave array
1497
"""
1498
ts = np.asarray(ts)
1499
cycles = self.freq * ts + self.offset / PI2
1500
frac, _ = np.modf(cycles)
1501
ys = np.abs(frac - 0.5)
1502
ys = normalize(unbias(ys), self.amp)
1503
return ys
1504
1505
from scipy.integrate import cumtrapz
1506
1507
class Chirp(Signal):
1508
"""Represents a signal with variable frequency."""
1509
1510
def __init__(self, start=440, end=880, amp=1.0):
1511
"""Initializes a linear chirp.
1512
1513
start: float frequency in Hz
1514
end: float frequency in Hz
1515
amp: float amplitude, 1.0 is nominal max
1516
"""
1517
self.start = start
1518
self.end = end
1519
self.amp = amp
1520
1521
@property
1522
def period(self):
1523
"""Period of the signal in seconds.
1524
1525
returns: float seconds
1526
"""
1527
return ValueError("Non-periodic signal.")
1528
1529
def evaluate(self, ts):
1530
"""Evaluates the signal at the given times.
1531
1532
Note: This version is a little more complicated than the one
1533
in the book because it handles the case where the ts are
1534
not equally spaced.
1535
1536
ts: float array of times
1537
1538
returns: float wave array
1539
"""
1540
def interpolate(ts, f0, f1):
1541
t0, t1 = ts[0], ts[-1]
1542
return f0 + (f1 - f0) * (ts - t0) / (t1 - t0)
1543
1544
# compute the frequencies
1545
ts = np.asarray(ts)
1546
freqs = interpolate(ts, self.start, self.end)
1547
1548
# compute the time intervals
1549
dts = np.diff(ts, append=ts[-1])
1550
1551
# compute the changes in phase
1552
dphis = PI2 * freqs * dts
1553
dphis = np.roll(dphis, 1)
1554
1555
# compute phase
1556
phases = np.cumsum(dphis)
1557
1558
# compute the amplitudes
1559
ys = self.amp * np.cos(phases)
1560
return ys
1561
1562
1563
class ExpoChirp(Chirp):
1564
"""Represents a signal with varying frequency."""
1565
1566
def evaluate(self, ts):
1567
"""Evaluates the signal at the given times.
1568
1569
ts: float array of times
1570
1571
returns: float wave array
1572
"""
1573
f0, f1 = np.log10(self.start), np.log10(self.end)
1574
freqs = np.logspace(f0, f1, len(ts))
1575
dts = np.diff(ts, prepend=0)
1576
dphis = PI2 * freqs * dts
1577
phases = np.cumsum(dphis)
1578
ys = self.amp * np.cos(phases)
1579
return ys
1580
1581
1582
class SilentSignal(Signal):
1583
"""Represents silence."""
1584
1585
def evaluate(self, ts):
1586
"""Evaluates the signal at the given times.
1587
1588
ts: float array of times
1589
1590
returns: float wave array
1591
"""
1592
return np.zeros(len(ts))
1593
1594
1595
class Impulses(Signal):
1596
"""Represents silence."""
1597
1598
def __init__(self, locations, amps=1):
1599
self.locations = np.asanyarray(locations)
1600
self.amps = amps
1601
1602
def evaluate(self, ts):
1603
"""Evaluates the signal at the given times.
1604
1605
ts: float array of times
1606
1607
returns: float wave array
1608
"""
1609
ys = np.zeros(len(ts))
1610
indices = np.searchsorted(ts, self.locations)
1611
ys[indices] = self.amps
1612
return ys
1613
1614
1615
class Noise(Signal):
1616
"""Represents a noise signal (abstract parent class)."""
1617
1618
def __init__(self, amp=1.0):
1619
"""Initializes a white noise signal.
1620
1621
amp: float amplitude, 1.0 is nominal max
1622
"""
1623
self.amp = amp
1624
1625
@property
1626
def period(self):
1627
"""Period of the signal in seconds.
1628
1629
returns: float seconds
1630
"""
1631
return ValueError("Non-periodic signal.")
1632
1633
1634
class UncorrelatedUniformNoise(Noise):
1635
"""Represents uncorrelated uniform noise."""
1636
1637
def evaluate(self, ts):
1638
"""Evaluates the signal at the given times.
1639
1640
ts: float array of times
1641
1642
returns: float wave array
1643
"""
1644
ys = np.random.uniform(-self.amp, self.amp, len(ts))
1645
return ys
1646
1647
1648
class UncorrelatedGaussianNoise(Noise):
1649
"""Represents uncorrelated gaussian noise."""
1650
1651
def evaluate(self, ts):
1652
"""Evaluates the signal at the given times.
1653
1654
ts: float array of times
1655
1656
returns: float wave array
1657
"""
1658
ys = np.random.normal(0, self.amp, len(ts))
1659
return ys
1660
1661
1662
class BrownianNoise(Noise):
1663
"""Represents Brownian noise, aka red noise."""
1664
1665
def evaluate(self, ts):
1666
"""Evaluates the signal at the given times.
1667
1668
Computes Brownian noise by taking the cumulative sum of
1669
a uniform random series.
1670
1671
ts: float array of times
1672
1673
returns: float wave array
1674
"""
1675
dys = np.random.uniform(-1, 1, len(ts))
1676
# ys = scipy.integrate.cumtrapz(dys, ts)
1677
ys = np.cumsum(dys)
1678
ys = normalize(unbias(ys), self.amp)
1679
return ys
1680
1681
1682
class PinkNoise(Noise):
1683
"""Represents Brownian noise, aka red noise."""
1684
1685
def __init__(self, amp=1.0, beta=1.0):
1686
"""Initializes a pink noise signal.
1687
1688
amp: float amplitude, 1.0 is nominal max
1689
"""
1690
self.amp = amp
1691
self.beta = beta
1692
1693
def make_wave(self, duration=1, start=0, framerate=11025):
1694
"""Makes a Wave object.
1695
1696
duration: float seconds
1697
start: float seconds
1698
framerate: int frames per second
1699
1700
returns: Wave
1701
"""
1702
signal = UncorrelatedUniformNoise()
1703
wave = signal.make_wave(duration, start, framerate)
1704
spectrum = wave.make_spectrum()
1705
1706
spectrum.pink_filter(beta=self.beta)
1707
1708
wave2 = spectrum.make_wave()
1709
wave2.unbias()
1710
wave2.normalize(self.amp)
1711
return wave2
1712
1713
1714
def rest(duration):
1715
"""Makes a rest of the given duration.
1716
1717
duration: float seconds
1718
1719
returns: Wave
1720
"""
1721
signal = SilentSignal()
1722
wave = signal.make_wave(duration)
1723
return wave
1724
1725
1726
def make_note(midi_num, duration, sig_cons=CosSignal, framerate=11025):
1727
"""Make a MIDI note with the given duration.
1728
1729
midi_num: int MIDI note number
1730
duration: float seconds
1731
sig_cons: Signal constructor function
1732
framerate: int frames per second
1733
1734
returns: Wave
1735
"""
1736
freq = midi_to_freq(midi_num)
1737
signal = sig_cons(freq)
1738
wave = signal.make_wave(duration, framerate=framerate)
1739
wave.apodize()
1740
return wave
1741
1742
1743
def make_chord(midi_nums, duration, sig_cons=CosSignal, framerate=11025):
1744
"""Make a chord with the given duration.
1745
1746
midi_nums: sequence of int MIDI note numbers
1747
duration: float seconds
1748
sig_cons: Signal constructor function
1749
framerate: int frames per second
1750
1751
returns: Wave
1752
"""
1753
freqs = [midi_to_freq(num) for num in midi_nums]
1754
signal = sum(sig_cons(freq) for freq in freqs)
1755
wave = signal.make_wave(duration, framerate=framerate)
1756
wave.apodize()
1757
return wave
1758
1759
1760
def midi_to_freq(midi_num):
1761
"""Converts MIDI note number to frequency.
1762
1763
midi_num: int MIDI note number
1764
1765
returns: float frequency in Hz
1766
"""
1767
x = (midi_num - 69) / 12.0
1768
freq = 440.0 * 2 ** x
1769
return freq
1770
1771
1772
def sin_wave(freq, duration=1, offset=0):
1773
"""Makes a sine wave with the given parameters.
1774
1775
freq: float cycles per second
1776
duration: float seconds
1777
offset: float radians
1778
1779
returns: Wave
1780
"""
1781
signal = SinSignal(freq, offset=offset)
1782
wave = signal.make_wave(duration)
1783
return wave
1784
1785
1786
def cos_wave(freq, duration=1, offset=0):
1787
"""Makes a cosine wave with the given parameters.
1788
1789
freq: float cycles per second
1790
duration: float seconds
1791
offset: float radians
1792
1793
returns: Wave
1794
"""
1795
signal = CosSignal(freq, offset=offset)
1796
wave = signal.make_wave(duration)
1797
return wave
1798
1799
1800
def mag(a):
1801
"""Computes the magnitude of a numpy array.
1802
1803
a: numpy array
1804
1805
returns: float
1806
"""
1807
return np.sqrt(np.dot(a, a))
1808
1809
1810
def zero_pad(array, n):
1811
"""Extends an array with zeros.
1812
1813
array: numpy array
1814
n: length of result
1815
1816
returns: new NumPy array
1817
"""
1818
res = np.zeros(n)
1819
res[: len(array)] = array
1820
return res
1821
1822
1823
def decorate(**options):
1824
"""Decorate the current axes.
1825
1826
Call decorate with keyword arguments like
1827
1828
decorate(title='Title',
1829
xlabel='x',
1830
ylabel='y')
1831
1832
The keyword arguments can be any of the axis properties
1833
1834
https://matplotlib.org/api/axes_api.html
1835
1836
In addition, you can use `legend=False` to suppress the legend.
1837
1838
And you can use `loc` to indicate the location of the legend
1839
(the default value is 'best')
1840
"""
1841
loc = options.pop("loc", "best")
1842
if options.pop("legend", True):
1843
legend(loc=loc)
1844
1845
plt.gca().set(**options)
1846
plt.tight_layout()
1847
1848
1849
def legend(**options):
1850
"""Draws a legend only if there is at least one labeled item.
1851
1852
options are passed to plt.legend()
1853
https://matplotlib.org/api/_as_gen/matplotlib.plt.legend.html
1854
1855
"""
1856
underride(options, loc="best", frameon=False)
1857
1858
ax = plt.gca()
1859
handles, labels = ax.get_legend_handles_labels()
1860
if handles:
1861
ax.legend(handles, labels, **options)
1862
1863
1864
def remove_from_legend(bad_labels):
1865
"""Removes some labels from the legend.
1866
1867
bad_labels: sequence of strings
1868
"""
1869
ax = plt.gca()
1870
handles, labels = ax.get_legend_handles_labels()
1871
handle_list, label_list = [], []
1872
for handle, label in zip(handles, labels):
1873
if label not in bad_labels:
1874
handle_list.append(handle)
1875
label_list.append(label)
1876
ax.legend(handle_list, label_list)
1877
1878
1879
def underride(d, **options):
1880
"""Add key-value pairs to d only if key is not in d.
1881
1882
If d is None, create a new dictionary.
1883
1884
d: dictionary
1885
options: keyword args to add to d
1886
"""
1887
if d is None:
1888
d = {}
1889
1890
for key, val in options.items():
1891
d.setdefault(key, val)
1892
1893
return d
1894
1895
1896
1897
1898
def main():
1899
1900
cos_basis = cos_wave(440)
1901
sin_basis = sin_wave(440)
1902
1903
wave = cos_wave(440, offset=math.pi / 2)
1904
cos_cov = cos_basis.cov(wave)
1905
sin_cov = sin_basis.cov(wave)
1906
print(cos_cov, sin_cov, mag((cos_cov, sin_cov)))
1907
return
1908
1909
wfile = WavFileWriter()
1910
for sig_cons in [
1911
SinSignal,
1912
TriangleSignal,
1913
SawtoothSignal,
1914
GlottalSignal,
1915
ParabolicSignal,
1916
SquareSignal,
1917
]:
1918
print(sig_cons)
1919
sig = sig_cons(440)
1920
wave = sig.make_wave(1)
1921
wave.apodize()
1922
wfile.write(wave)
1923
wfile.close()
1924
return
1925
1926
signal = GlottalSignal(440)
1927
signal.plot()
1928
plt.show()
1929
return
1930
1931
wfile = WavFileWriter()
1932
for m in range(60, 0, -1):
1933
wfile.write(make_note(m, 0.25))
1934
wfile.close()
1935
return
1936
1937
wave1 = make_note(69, 1)
1938
wave2 = make_chord([69, 72, 76], 1)
1939
wave = wave1 | wave2
1940
1941
wfile = WavFileWriter()
1942
wfile.write(wave)
1943
wfile.close()
1944
return
1945
1946
sig1 = CosSignal(freq=440)
1947
sig2 = CosSignal(freq=523.25)
1948
sig3 = CosSignal(freq=660)
1949
sig4 = CosSignal(freq=880)
1950
sig5 = CosSignal(freq=987)
1951
sig = sig1 + sig2 + sig3 + sig4
1952
1953
# wave = Wave(sig, duration=0.02)
1954
# wave.plot()
1955
1956
wave = sig.make_wave(duration=1)
1957
# wave.normalize()
1958
1959
wfile = WavFileWriter(wave)
1960
wfile.write()
1961
wfile.close()
1962
1963
1964
if __name__ == "__main__":
1965
main()
1966
1967