# -*- coding: utf-8 -*-
"""
Created on Fri Dec 23 13:36:45 2022
@author: xxxx
"""
from scipy import signal
from scipy.io.wavfile import read
import matplotlib.pyplot as plt
import numpy as np
#edit here to add a HOME directory
#HOME = "/home/david/" #Linux example
HOME = 'C:/Users/xxxx/Documents/data/Hifi/Michell/RPM WF/wf011.wav' #Windows example
#edit here for 45 rpm samples
Period = 1.8 # 33rpm
#Period = 4/3 # 45 rpm
#end edit
def instfreq(sig,Fs):
z = signal.hilbert(sig)
rawfreq = Fs/(2*np.pi)*np.diff(np.unwrap(np.angle(z)))
rawfreq = np.append(rawfreq,rawfreq[len(rawfreq)-1]) #np.diff drops one end point
b, a = signal.iirfilter(1,100./(Fs/2), btype='lowpass')
a = np.convolve(a,a) #make an ordinary 100Hz 4 pole filter.
a = np.convolve(a,a)
b = np.convolve(b,b)
b = np.convolve(b,b)
instfreq = signal.lfilter(b,a,rawfreq)
return (instfreq)
_FILE = input('wf011.wav') #file needs to be >= 4 sec.
y = read(HOME + _FILE)
Fs = float(y[0])
if np.size(y[1][0]) == 2:
sig = y[1][:,0][0:int(Fs*4)] #Grab 4 sec of a stereo file from the first (left?) channel
else:
sig = y[1][0:int(Fs*4)] #mono file
t = np.arange(Period,0,-1/Fs) #Reverse time (theta axis)
theta = t*2*np.pi/Period #Time becomes degrees (Periodsec = 2pi radians)
theta = np.roll(theta,int(Fs*.45)) #Rotate 90 deg to put 0 on top (Period*Fs/4)
freq1 = instfreq(sig,Fs)
glitch = 2000 #remove leading glitch from output (probably should depend on Fs)
if1 = freq1[glitch:glitch+int(Fs*Period)]
if2 = freq1[glitch+int(Fs*Period):glitch+int(2*Fs*Period)]
maxf = int(max(max(if1),max(if2))+.5)
r1 = 20.-(maxf-if1)/3. #20 radial ticks at 3Hz is fixed, adaptive scaling
r2 = 20.-(maxf-if2)/3. #is an exercise for later
ax = plt.subplot(111, projection='polar')
ax.plot(theta,r1)
ax.plot(theta,r2)
ax.set_rmax(20)
#Set up the ticks y is radial x is theta, it turns out x and y
#methods work in polar projection but sometimes do funny things
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(8)
plt.yticks(np.arange(1,21,1))
ax.set_rlabel_position(90)
myticks = []
for x in range(0,20,1):
myticks.append(str(maxf-57+x*3))
ax.set_yticklabels(myticks)
ax.set_xticklabels(['90'+u'\N{DEGREE SIGN}','45'+u'\N{DEGREE SIGN}','0'+u'\N{DEGREE SIGN}',\
'315'+u'\N{DEGREE SIGN}','270'+u'\N{DEGREE SIGN}','225'+u'\N{DEGREE SIGN}',\
'180'+u'\N{DEGREE SIGN}','135'+u'\N{DEGREE SIGN}'])
ax.grid(True)
ax.set_title("Instantaneous Frequency "+ _FILE, va='bottom')
plt.show()