%pylab inline
#rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
import wave
import glob
from scipy.io import wavfile
from scipy.signal import freqs
Paths=glob.glob('C:/Users/Meera/AppData/Local/Temp/media/*.wav')
print(Paths) #Reading paths of audio in this order - Goleckenspiel, Piano, Tom
['C:/Users/Meera/AppData/Local/Temp/media\\glockenspiel.wav', 'C:/Users/Meera/AppData/Local/Temp/media\\piano.wav', 'C:/Users/Meera/AppData/Local/Temp/media\\tom.wav']
#Function to calculate the period(in terms of lag) between the peaks in the autocorrelated signal
def returnperiod(sr,samples):
peak=0
peakarray=[]
for i in range(1,len(samples)): #checking for peak starting with lag=1 (0 lag is the same signal)
if (c1[i-1]<c1[i] and c1[i]>c1[i+1]) and (peak<2) :
print("Lag index of peak ", peak +1, "is", i)
peakarray.append(i)
plt.axvline(i,color='r',label='{:d}'.format(i))
legend()
peak=peak+1
if peak==2: #we only need 2 peaks
break
return peakarray
def findfreq (sr,samples): #Defining a function to calculate the frequency
N=len(samples)
num=(N/2)+1
X=fft.rfft(samples,len(samples)) #Computing one real FFT with frame size equal to the number of samples for each file.
Magspec=abs(X)/(N/2)
nyquist = sr/2.0
freq = linspace(0, nyquist, num, endpoint=True) #Converting to frequency scale
Totalenergy=sum(Magspec)
lags, c, lines, line = acorr(samples, maxlags=None);
Cutoffbin=where(cumsum(Magspec)>=0.95*Totalenergy)[0][0] # creating a cut off to restrict signal to 95% of total energy
cutofffreq=Cutoffbin*sr/len(samples)
subplot(121)
plot(freq, Magspec)
xlabel('frequency')
ylabel('magnitude')
plt.axvline(cutofffreq,color='r',label='{:2f}Hz'.format(cutofffreq)) #Red line indicated the cut off frequency
legend()
subplot(122)
plot(freq,20*log10(Magspec/Totalenergy)) #Converting to db scale
xlim(0,cutofffreq)
xlabel('frequency')
ylabel('magnitude-db')
gcf().set_figwidth(18)
peakbin=argmax(Magspec) #Finding maximum argument
peakfreq=peakbin*sr/len(samples) #Converting to frequency in Hz
print('The peak frequency is {:f}Hz'.format(peakfreq))
sr, samples = wavfile.read(Paths[0])
samples=samples[:].astype(double)
print(Paths[0])
lags, c, lines, line = acorr(samples, maxlags=None,usevlines=True); #Performing autocorrelation of the signal
grid()
gcf().set_figwidth(18)
c1=c[samples.shape[0]:] #We do not need the negative lags, hence only considering the posive lags
subplot(121)
plot(c1)
xlabel('lag index')
ylabel('magnitude')
subplot(122)
plot(c1[0:200]) #zooming into the autocorrelated signal
xlabel('lag index')
ylabel('magnitude')
gcf().set_figwidth(18)
peakarray=returnperiod(sr,samples) #to find the period, we subtract the index of first peak from second peak
period=peakarray[1]-peakarray[0]
fundfreq=sr/period #To find the fundamental frequency-> fs/period
print("The fundamental frequency is",fundfreq, "Hz")
C:/Users/Meera/AppData/Local/Temp/media\glockenspiel.wav Lag index of peak 1 is 32 Lag index of peak 2 is 65 The fundamental frequency is 1336.36363636 Hz
findfreq(sr,samples)
The peak frequency is 1323.342835Hz
sr, samples = wavfile.read(Paths[2])
samples=samples[:].astype(double)
print(Paths[2])
lags, c, lines, line = acorr(samples, maxlags=None,usevlines=True);
grid()
gcf().set_figwidth(18)
c1=c[samples.shape[0]:] #We do not need the negative lags, hence only considering the posive lags
subplot(121)
plot(c1)
xlabel('lag index')
ylabel('magnitude')
subplot(122)
plot(c1[0:1000]) #zooming into the autocorrelated signal
xlabel('lag index')
ylabel('magnitude')
gcf().set_figwidth(18)
peakarray=returnperiod(sr,samples)
period=peakarray[1]-peakarray[0] #to find the period, we subtract the index of first peak from second
fundfreq=sr/period #To find the fundamental frequency-> fs/period
print("The fundamental frequency is",fundfreq , "Hz")
C:/Users/Meera/AppData/Local/Temp/media\tom.wav Lag index of peak 1 is 467 Lag index of peak 2 is 940 The fundamental frequency is 93.2346723044 Hz
findfreq(sr,samples)
The peak frequency is 92.565431Hz
from scipy.signal import lfilter
sr, samples = wavfile.read(Paths[1])
samples=samples[:].astype(double)
print(Paths[1])
lags, c, lines, line = acorr(samples, maxlags=None,usevlines=True);grid()
gcf().set_figwidth(18)
c1=c[samples.shape[0]:]
subplot(131)
plot(c1)
xlabel('lag index')
ylabel('magnitude')
subplot(132)
plot(c1[0:1200])
xlabel('lag index')
ylabel('magnitude')
from scipy import signal
b, a = signal.butter(1, 0.004, 'low', analog=False) #Using a butterworth low pass filter, which passes frequencies upto around 88Hz
c1 = filtfilt(b, a, c1)
subplot(133)
plot(c1[0:1200])
title('after low pass filtering the autocorrelated signal')
gcf().set_figwidth(18)
peakarray=returnperiod(sr,samples) #After filtering, it is easier to calculate the fundamental frequency (Fig 3)
period=peakarray[1]-peakarray[0]
fundfreq=sr/period
print("The fundamental frequency is",fundfreq)
C:/Users/Meera/AppData/Local/Temp/media\piano.wav Lag index of peak 1 is 566 Lag index of peak 2 is 1133 The fundamental frequency is 77.7777777778
findfreq(sr,samples)
The peak frequency is 78.057861Hz