python音频处理thinkdsp.zip

from __future__ import print_function, division import array import copy import math import numpy as np import random import scipy import scipy.stats import scipy.fftpack import struct import subprocess import thinkplot import warnings from fractions import gcd from wave import open as open_wave import matplotlib.pyplot as pyplot try: from IPython.display import Audio except: warnings.warn( "Can't import Audio from IPython.display; " "Wave.make_audio() will not work." ) PI2 = math.pi * 2 def random_seed(x): """Initialize the random and np.random generators. x: int seed """ random.seed(x) np.random.seed(x) class UnimplementedMethodException(Exception): """Exception if someone calls a method that should be overridden.""" class WavFileWriter: """Writes wav files.""" def __init__(self, filename="sound.wav", framerate=11025): """Opens the file and sets parameters. filename: string framerate: samples per second """ self.filename = filename self.framerate = framerate self.nchannels = 1 self.sampwidth = 2 self.bits = self.sampwidth * 8 self.bound = 2 ** (self.bits - 1) - 1 self.fmt = "h" self.dtype = np.int16 self.fp = open_wave(self.filename, "w") self.fp.setnchannels(self.nchannels) self.fp.setsampwidth(self.sampwidth) self.fp.setframerate(self.framerate) def write(self, wave): """Writes a wave. wave: Wave """ zs = wave.quantize(self.bound, self.dtype) self.fp.writeframes(zs.tostring()) def close(self, duration=0): """Closes the file. duration: how many seconds of silence to append """ if duration: self.write(rest(duration)) self.fp.close() def read_wave(filename="sound.wav"): """Reads a wave file. filename: string returns: Wave """ fp = open_wave(filename, "r") nchannels = fp.getnchannels() nframes = fp.getnframes() sampwidth = fp.getsampwidth() framerate = fp.getframerate() z_str = fp.readframes(nframes) fp.close() dtype_map = {1: np.int8, 2: np.int16, 3: "special", 4: np.int32} if sampwidth not in dtype_map: raise ValueError("sampwidth %d unknown" % sampwidth) if sampwidth == 3: xs = np.fromstring(z_str, dtype=np.int8).astype(np.int32) ys = (xs[2::3] * 256 + xs[1::3]) * 256 + xs[0::3] else: ys = np.fromstring(z_str, dtype=dtype_map[sampwidth]) # if it's in stereo, just pull out the first channel if nchannels == 2: ys = ys[::2] # ts = np.arange(len(ys)) / framerate wave = Wave(ys, framerate=framerate) wave.normalize() return wave def play_wave(filename="sound.wav", player="aplay"): """Plays a wave file. filename: string player: string name of executable that plays wav files """ cmd = "%s %s" % (player, filename) popen = subprocess.Popen(cmd, shell=True) popen.communicate() def find_index(x, xs): """Find the index corresponding to a given value in an array.""" n = len(xs) start = xs[0] end = xs[-1] i = round((n - 1) * (x - start) / (end - start)) return int(i) class _SpectrumParent: """Contains code common to Spectrum and DCT. """ def __init__(self, hs, fs, framerate, full=False): """Initializes a spectrum. hs: array of amplitudes (real or complex) fs: array of frequencies framerate: frames per second full: boolean to indicate full or real FFT """ self.hs = np.asanyarray(hs) self.fs = np.asanyarray(fs) self.framerate = framerate self.full = full @property def max_freq(self): """Returns the Nyquist frequency for this spectrum.""" return self.framerate / 2 @property def amps(self): """Returns a sequence of amplitudes (read-only property).""" return np.absolute(self.hs) @property def power(self): """Returns a sequence of powers (read-only property).""" return self.amps ** 2 def copy(self): """Makes a copy. Returns: new Spectrum """ return copy.deepcopy(self) def max_diff(self, other): """Computes the maximum absolute difference between spectra. other: Spectrum returns: float """ assert self.framerate == other.framerate assert len(self) == len(other) hs = self.hs - other.hs return np.max(np.abs(hs)) def ratio(self, denom, thresh=1, val=0): """The ratio of two spectrums. denom: Spectrum thresh: values smaller than this are replaced val: with this value returns: new Wave """ ratio_spectrum = self.copy() ratio_spectrum.hs /= denom.hs ratio_spectrum.hs[denom.amps < thresh] = val return ratio_spectrum def invert(self): """Inverts this spectrum/filter. returns: new Wave """ inverse = self.copy() inverse.hs = 1 / inverse.hs return inverse @property def freq_res(self): return self.framerate / 2 / (len(self.fs) - 1) def render_full(self, high=None): """Extracts amps and fs from a full spectrum. high: cutoff frequency returns: fs, amps """ hs = np.fft.fftshift(self.hs) amps = np.abs(hs) fs = np.fft.fftshift(self.fs) i = 0 if high is None else find_index(-high, fs) j = None if high is None else find_index(high, fs) + 1 return fs[i:j], amps[i:j] def plot(self, high=None, **options): """Plots amplitude vs frequency. Note: if this is a full spectrum, it ignores low and high high: frequency to cut off at """ if self.full: fs, amps = self.render_full(high) thinkplot.plot(fs, amps, **options) else: i = None if high is None else find_index(high, self.fs) thinkplot.plot(self.fs[:i], self.amps[:i], **options) def plot_power(self, high=None, **options): """Plots power vs frequency. high: frequency to cut off at """ if self.full: fs, amps = self.render_full(high) thinkplot.plot(fs, amps ** 2, **options) else: i = None if high is None else find_index(high, self.fs) thinkplot.plot(self.fs[:i], self.power[:i], **options) def estimate_slope(self): """Runs linear regression on log power vs log frequency. returns: slope, inter, r2, p, stderr """ x = np.log(self.fs[1:]) y = np.log(self.power[1:]) t = scipy.stats.linregress(x, y) return t def peaks(self): """Finds the highest peaks and their frequencies. returns: sorted list of (amplitude, frequency) pairs """ t = list(zip(self.amps, self.fs)) t.sort(reverse=True) return t class Spectrum(_SpectrumParent): """Represents the spectrum of a signal.""" def __len__(self): """Length of the spectrum.""" return len(self.hs) def __add__(self, other): """Adds two spectrums elementwise. other: Spectrum returns: new Spectrum """ if other == 0: return self.copy() assert all(self.fs == other.fs) hs = self.hs + other.hs return Spectrum(hs, self.fs, self.framerate, self.full) __radd__ = __add__ def __mul__(self, other): """Multiplies two spectrums elementwise.