## MAT 201A Winter 2016
## HW 4
## Mark Hirsch
## my solution to the homework draws from notes taken during the help session provided by Owen for this assignment.
%pylab inline
from __future__ import print_function
from __future__ import division
import wave
from IPython.display import Audio
from scipy.io import wavfile
import glob
import IPython.display
Populating the interactive namespace from numpy and matplotlib
##load first sound file
rate, x= wavfile.read('glockenspiel.wav')
##get size of file and setup bins to include Nyquist freq
N = x.size
numBins = int(N/2+1)
x = x/2.0**15
##plot and play file
plot(x);
title('glockenspiel')
Audio('glockenspiel.wav')
##do the real DFT, get the magnitude spectrum & frequency (and plot them)
x = fft.rfft(x)
magSpec = abs(x)
freqs = rate * array([i / float(N) for i in range(numBins)])
plot(freqs,magSpec);
xlabel('freq')
ylabel('magSpec')
<matplotlib.text.Text at 0x10ea44810>
##establish a cut off point/freq (to highlight area of 95% of energy)
energy = sum(magSpec)
cutOff = where(cumsum(magSpec) >= 0.95 * energy)[0][0]
cutOffFreq = freqs[cutOff]
##show the region in which 95% of the energy is
plot(freqs,magSpec);
xlabel('freq')
ylabel('mag')
##nifty formatting trick Owen showed us
axvline(cutOffFreq, color='g',linewidth=4.0,alpha=0.25,label='{0:.2f} Hz'.format(cutOffFreq))
legend();
##limit the plot to this 95% energy area
plot(freqs,magSpec)
xlabel('freq')
ylabel('mag')
xlim(0,cutOffFreq);
##get the "most prominent peak" and show with green line
peak = argmax(magSpec)
assert(peak <= cutOff)
peak = freqs[peak]
plot(freqs,magSpec)
xlim(0,cutOffFreq)
xlabel('freq')
ylabel('mag')
axvline(peak,color='g',linewidth=4.0,alpha=0.25,label='{0:.2f} Hz'.format(cutOffFreq))
legend();
##as per recommendation by Owen, make a function of this process and run
##it for other two audio files
def getBigPeak(soundFile):
##load first sound file
rate, x= wavfile.read(soundFile)
##get size of file and setup bins to include Nyquist freq
N = x.size
numBins = int(N/2+1)
x = x/2.0**15
##plot and play file
plot(x);
title(soundFile)
Audio(soundFile)
##do the real DFT, get the magnitude spectrum & frequency (and plot them)
x = fft.rfft(x)
magSpec = abs(x)
freqs = rate * array([i / float(N) for i in range(numBins)])
plot(freqs,magSpec);
xlabel('freq')
ylabel('magSpec')
##establish a cut off point/freq (to highlight area of 95% of energy)
energy = sum(magSpec)
cutOff = where(cumsum(magSpec) >= 0.95 * energy)[0][0]
cutOffFreq = freqs[cutOff]
##show the region in which 95% of the energy is
plot(freqs,magSpec);
xlabel('freq')
ylabel('mag')
##nifty formatting trick Owen showed us
axvline(cutOffFreq, color='g',linewidth=4.0,alpha=0.25,label='{0:.2f} Hz'.format(cutOffFreq))
legend();
##limit the plot to this 95% energy area
plot(freqs,magSpec)
xlabel('freq')
ylabel('mag')
xlim(0,cutOffFreq);
##get the "most prominent peak" and show with green line
peak = argmax(magSpec)
assert(peak <= cutOff)
peak = freqs[peak]
plot(freqs,magSpec)
xlim(0,cutOffFreq)
xlabel('freq')
ylabel('mag')
axvline(peak,color='g',linewidth=4.0,alpha=0.25,label='{0:.2f} Hz'.format(cutOffFreq))
legend();
getBigPeak('piano.wav')
getBigPeak('tom.wav')