"""
PLEASE READ
To use this script you need to edit HOME to the directory where the .wav file is.
The .wav file can be any monotonic frequency sweep either up or down in frequency but
it must be trimmed at both ends to remove any leading silence. Frequencies below 1kHz are
ignored since virually all cartridges are well behaved below 1kHz.
Python does not read 24bit packed .wav's 16 bit is more than enough here.
The info_line should be alpha-numeric with entries separated by " / " only. The script
will save a .png file that is named from the info line, replacing " / " with "_". As
example "this / is / a / test" will create a file named "this_is_a_test.png"
"""
from scipy import signal
from scipy.io.wavfile import read
import matplotlib.pyplot as plt
import numpy as np
import datetime
import os
#edit here to add a HOME directory, etc.
HOME = 'Users\16073\Documents\Audacity'
_FILE = 'oc9xml.wav'
info_line = 'oc9xml'
#end edit
try:
input = raw_input
except NameError:
pass
#_FILE = input('Enter a .wav file (no quotes): ') #file needs to be >= sec_per_rev*2
#Needed to steal Matlab's flat top window
def ft_window(n):
w = []
a0 = 0.21557895
a1 = 0.41663158
a2 = 0.277263158
a3 = 0.083578947
a4 = 0.006947368
pi = np.pi
for x in range(0,n):
w.append(a0 - a1*np.cos(2*pi*x/(n-1)) + a2*np.cos(4*pi*x/(n-1)) - a3*np.cos(6*pi*x/(n-1)) + a4*np.cos(8*pi*x/(n-1)))
return w
try:
input = raw_input
except NameError:
pass
#_FILE = input('Enter a .wav file to process ')
#File needs trimming with no leading or trailing silence (undefined f)
y = read(HOME + _FILE)
Fs = float(y[0])
input_sig = y[1]
frequency = []
amplitude = []
frequency2h = []
amplitude2h = []
frequency3h = []
amplitude3h = []
freqout = []
ampout = []
freqout2h = []
ampout2h = []
freqout3h = []
ampout3h = []
F = int(Fs/100)
win = ft_window(F)
#For Fs sampling use small Fs/100 FFT's so frequency is f/100
#All common sampling rates are divisible by 100
#First try, no overlap and flat top window
#Find the minimum and maximum frequencies and flip the data if the sweep
#is going down in frequency.
y = abs(np.fft.rfft(input_sig[0:F]*win))
minf = np.argmax(y)
y = abs(np.fft.rfft(input_sig[len(input_sig)-F:len(input_sig)]*win))
maxf = np.argmax(y)
if maxf < minf:
maxf,minf = minf,maxf
input_sig = np.flipud(input_sig)
for x in range(0,len(input_sig)-F,F):
y = abs(np.fft.rfft(input_sig[x:x+F]*win))
f = np.argmax(y) #use largest bin
if f >=10:
frequency.append(f*100)
amplitude.append(y[f])
if 2*f<F/2-2 and f >= 10:
f2 = np.argmax(y[(2*f)-2:(2*f)+2])
frequency2h.append(f*100)
amplitude2h.append(y[2*f-2+f2])
if 3*f<F/2-2 and f >= 10:
f3 = np.argmax(y[(3*f)-2:(3*f)+2])
frequency3h.append(f*100)
amplitude3h.append(y[3*f-2+f3])
amp = 0
count = 0
x = minf*100
for x in range(minf*100,(maxf+1)*100,100):
for y in range(0,len(frequency)):
if frequency[y] == x:
amp = amp + amplitude[y]
count = count + 1
if count != 0:
freqout.append(x)
ampout.append(20*np.log10(amp/count))
amp = 0
count = 0
amp = 0
count = 0
x = minf*100
for x in range(minf*100,(maxf+1)*100,100):
for y in range(0,len(frequency2h)):
if frequency2h[y] == x:
amp = amp + amplitude2h[y]
count = count + 1
if count != 0:
freqout2h.append(x)
ampout2h.append(20*np.log10(amp/count))
amp = 0
count = 0
amp = 0
count = 0
x = minf*100
for x in range(minf*100,(maxf+1)*100,100):
for y in range(0,len(frequency3h)):
if frequency3h[y] == x:
amp = amp + amplitude3h[y]
count = count + 1
if count != 0:
freqout3h.append(x)
ampout3h.append(20*np.log10(amp/count))
amp = 0
count = 0
b, a = signal.iirfilter(3,.5, btype='lowpass') #filter some noise
ampout = signal.filtfilt(b,a,ampout)
ampout2h = signal.filtfilt(b,a,ampout2h)
ampout3h = signal.filtfilt(b,a,ampout3h)
norm = ampout[0]
ampout = ampout-norm #amplitude is in dB so normalize by subtraction at [0]
ampout2h = ampout2h-norm
ampout3h = ampout3h-norm
f, plt.figure(figsize=(14,6))
plt.semilogx(freqout,ampout,color = 'b')
plt.semilogx(freqout2h,ampout2h,color = 'g')
plt.semilogx(freqout3h,ampout3h,color = 'r')
plt.grid(True, which="both", axis="x", ls="-", color="black")
plt.grid(True, which="major", axis="both", ls="-", color="black")
plt.grid(True, which="minor", axis="both", ls="-", color="gainsboro")
#plt.grid(True, which="both", axis="x", ls="-", color="black")
#plt.xaxis.grid(True, which="minor", ls="-", color="gainsboro")
plt.minorticks_on()
plt.title("Instantaneous Frequency vs. Amplitude and Distortion: "+ _FILE + "\n", fontsize=16)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude (dB)")
plt.figtext(.17,.7,'- 2nd harmonic -',color = 'g')
plt.figtext(.17,.67,'- 3rd harmonic -',color = 'r')
plt.tick_params(which='both',
top='off',
left='off',
right='off',
bottom='off')
mod_date = datetime.datetime.fromtimestamp(os.path.getmtime(_FILE))
plt.figtext(.17, .13, info_line + "\n" + \
mod_date.strftime("%b %d, %Y %H:%M:%S"), fontsize=10)
plt.savefig(info_line.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.75, dpi=96)
#plt.show()