from scipy.io import wavfile (rate, tone) = wavfile.read('files/Low_tone.wav') #rate is the sampling frequency samples = len(tone) samples NFFT = 1024 # the length of the windowing segments # Pxx is the segments x freqs array of instantaneous power, freqs is # the frequency vector, bins are the centers of the time bins in which # the power is computed, and im is the matplotlib.image.AxesImage # instance figure(figsize=(10, 8)) ax1 = subplot(211) plot(linspace(0, 1, samples) * samples / rate, tone) xlabel('time (seconds)') subplot(212, sharex=ax1) Pxx, freqs, bins, im = specgram(tone, NFFT=NFFT, Fs=rate, noverlap=900, cmap=cm.gist_heat) xlabel('time (seconds)') 4000/82. print Pxx.shape, freqs.shape energy = sum(Pxx ** 2, axis=0) energy.shape figure(figsize=(10, 4)) plot(bins, energy) xlim(0, 4) xticks(arange(0, 4, 0.25)); grid(True) xlabel("time (s)"); argmax(bins > 0.75) ampl = Pxx[:, 45] plot(freqs, ampl) xlabel("frequency (Hz)") title("Initial frequency content of the low E-string") plot(freqs, log(ampl)) xlabel("frequency (Hz)") title("Initial frequency content of the low E-string") from scipy.signal import hilbert signal = log(ampl) envelope = real(hilbert(signal)) plot(freqs, signal, label='original curve') plot(freqs, envelope, label='hilbert envelope') legend() plot(freqs, log(ampl)) xlabel("frequency (Hz)") title("Initial frequency content of the low E-string, zoomed on lower frequencies") xlim(0, 600) df = freqs[1] - freqs[0] df from scipy.signal import cspline1d, cspline1d_eval harmonics_freqs = arange(82.5, 4000, 82.5) cspline = cspline1d(ampl) harmonics_ampl = cspline1d_eval(cspline, harmonics_freqs, dx=df,x0=freqs[0]) plot(harmonics_freqs, log(harmonics_ampl), 'o') plot(freqs, log(ampl)) def gen_harmonic_sound(t, harmonics_freqs, harmonics_ampl): #normalizing step harmonics_ampl = harmonics_ampl / max(harmonics_ampl) signal = zeros(t.shape) for i in range(len(harmonics_freqs)): signal += harmonics_ampl[i] * sin(2 * pi * harmonics_freqs[i] * t) return signal t = arange(0, 5, 1/8000.) signal = gen_harmonic_sound(t, harmonics_freqs, harmonics_ampl) plot(t, signal) xlim(0, 0.5) import StringIO import base64 import struct from IPython.core.display import HTML def wavPlayer(data, rate): """ will display html 5 player for compatible browser The browser need to know how to play wav through html5. there is no autoplay to prevent file playing when the browser opens Adapted from SciPy.io. """ buffer = StringIO.StringIO() buffer.write(b'RIFF') buffer.write(b'\x00\x00\x00\x00') buffer.write(b'WAVE') buffer.write(b'fmt ') if data.ndim == 1: noc = 1 else: noc = data.shape[1] bits = data.dtype.itemsize * 8 sbytes = rate*(bits // 8)*noc ba = noc * (bits // 8) buffer.write(struct.pack('' or (data.dtype.byteorder == '=' and sys.byteorder == 'big'): data = data.byteswap() buffer.write(data.tostring()) # return buffer.getvalue() # Determine file size and place it in correct # position at start of the file. size = buffer.tell() buffer.seek(4) buffer.write(struct.pack(' Simple Test

""".format(base64=base64.encodestring(val)) display(HTML(src)) wavPlayer(2**13 * signal.astype(np.int16), 8000) figure(figsize=(10, 4)) plot(bins, energy) xlim(0, 4) xticks(arange(0, 4, 0.25)); grid(True) xlabel("time (s)"); from scipy import optimize selection = bins > 0.75 fit_time = bins[selection] fit_points = energy[selection] fitfunc = lambda p, x: p[0]*exp(-p[1] * x) # Target function errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target function p0 = [1e11, -1.] # Initial guess for the parameters p1, success = optimize.leastsq(errfunc, p0[:], args=(fit_time, fit_points)) time = linspace(fit_time[0], fit_time[-1], num=50) plot(fit_time, fit_points, label='data') # Plot of the data plot(time, fitfunc(p1, time), label = 'fit') # and the fit legend() p1 t = arange(0, 5, 1/8000.) signal = gen_harmonic_sound(t, harmonics_freqs, harmonics_ampl) signal = signal * fitfunc(p1, t) signal = signal / norm(signal, ord=inf) * 32000 signal = signal.astype(np.int16) wavPlayer(signal, 8000) plot(t, signal) NFFT = 1024 # the length of the windowing segments # Pxx is the segments x freqs array of instantaneous power, freqs is # the frequency vector, bins are the centers of the time bins in which # the power is computed, and im is the matplotlib.image.AxesImage # instance figure(figsize=(10, 8)) ax1 = subplot(211) plot(linspace(0, 1, samples) * samples / rate, tone) xlabel('time (seconds)') subplot(212, sharex=ax1) Pxx, freqs, bins, im = specgram(tone, NFFT=NFFT, Fs=rate, noverlap=900, cmap=cm.gist_heat) xlabel('time (seconds)') figure(figsize=(10, 4)) plot(freqs, ampl) xlabel("frequency (Hz)") title("Initial frequency content of the low E-string") xlim(0, 1000) xticks(arange(0, 1000, 50)); plot(arange(82, 1000, 82), 100000 * ones(arange(82, 1000, 82).shape), 'o') argmax(freqs >= 330) argmax(freqs >= 2* 82.5) plot(Pxx[22, :], label='2nd harmonic') plot(Pxx[43, :], label='4th harmonic') legend() def plot_decay(harmonic): ind = argmax(freqs >= harmonic* 82.5) decay = Pxx[ind, :] decay = decay / max(decay) label = "harmonic %s" % harmonic plot(bins, decay, label=label) plot_decay(2) plot_decay(4) plot_decay(5) plot_decay(7) plot_decay(10) legend() xlabel("time (s)") xlim(0, 8) p1 p2 = array([1., 1.]) normed_fitfunc = lambda p, x: fitfunc(p, x) / p[0] figure(figsize=(16, 4)) fit_values = normed_fitfunc(p1, bins) plot(bins+0.75, fit_values, label='exp decay t=2.17 s') plot(bins+0.75, normed_fitfunc(p2, bins), label='exp decay t=1 s') plot_decay(2) plot_decay(4) plot_decay(5) plot_decay(7) plot_decay(10) plot_decay(20) legend() xlabel("time (s)") xlim(0, 8) def gen_harmonic_sound_simple_decay(t, harmonics_freqs, harmonics_ampl): #normalizing step harmonics_ampl = harmonics_ampl / max(harmonics_ampl) signal = zeros(t.shape) for i in range(len(harmonics_freqs)): if i < 5: signal += harmonics_ampl[i] * sin(2 * pi * harmonics_freqs[i] * t) * exp(-t / 2.17) else: signal += harmonics_ampl[i] * sin(2 * pi * harmonics_freqs[i] * t) * exp(-t / 1.0) return signal t = arange(0, 5, 1/8000.) signal = gen_harmonic_sound_simple_decay(t, harmonics_freqs, harmonics_ampl) signal = signal * fitfunc(p1, t) signal = signal / norm(signal, ord=inf) * 32000 signal = signal.astype(np.int16) wavPlayer(signal, 8000) figure(figsize=(16, 4)) for harmonic in [1, 2, 4, 5, 10]: plot_decay(harmonic) decay_time = 2.17 / 5. * harmonic label = "exp decay tau=%.2f" % (decay_time) plot(bins+0.75, normed_fitfunc([1.0, decay_time], bins), label=label) legend() print decay_time label = "tau %.2f" % (decay_time) def gen_harmonic_sound_linear_decay(t, harmonics_freqs, harmonics_ampl): #normalizing step harmonics_ampl = harmonics_ampl / max(harmonics_ampl) signal = zeros(t.shape) for i in range(len(harmonics_freqs)): decay_time = 2.17 / 4.5 * i signal += harmonics_ampl[i] * sin(2 * pi * harmonics_freqs[i] * t) * exp(-t / decay_time) return signal t = arange(0, 5, 1/8000.) signal = gen_harmonic_sound_linear_decay(t, harmonics_freqs, harmonics_ampl) plot(t, signal) signal = (15000*signal).astype(np.int16) wavPlayer(signal, 8000)