# Turntable Wow and Flutter "Polar" Plot Script

#### JP

##### Major Contributor
This script creates polar plots from wow and flutter test tracks. I'll provide more background information on it later. The current "released" code is below.

Python:
``````#Version 5C

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.

_FILE = 'wavefile.wav'

info_line = 'PLOT TITLE'

Hz_per_tick = 3
sec_per_rev = 1.8  # 1.8 for 33, 1.33 for 45, .766 for 78.26

stereo_channel = 0    #0 = Left, 1 = Right
filter_freq = 60

#end edit

def instfreq(sig,Fs,filter_freq):
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,filter_freq/(Fs/2), btype='lowpass')
zi = signal.lfilter_zi(b, a) #Initialize the filter to the mean of the leading edge of the data
rawfreq,_ = signal.lfilter(b,a,rawfreq,zi=zi*np.mean(rawfreq[0:2000])) #reduces glitch, first pass

b, a = signal.iirfilter(3,filter_freq/(Fs/2), btype='lowpass')
instfreq = signal.filtfilt(b,a,rawfreq) #back and forth linear phase IIR filter (6 pole)

return (instfreq)

Fs = float(y[0])
if np.size(y[1][0]) == 2:
sig = y[1][:,stereo_channel][0:int(Fs*(sec_per_rev*3))] #Grab 3*sec_per_rev of audio from the specified channel
else:
sig = y[1][0:int(Fs*(sec_per_rev*3))] #mono file

t = np.arange(sec_per_rev,0,-1/Fs)  #Reverse time (theta axis)
theta = t*2*np.pi/sec_per_rev   #Time becomes degrees (1 rev = 2pi radians)
theta = np.roll(theta,int((sec_per_rev*Fs/4)))  #Rotate 90 deg to put 0 on top (sec_per_rev*Fs/4)

freq1 = instfreq(sig,Fs,filter_freq)

freq1 = np.roll(freq1,-int(Fs*.2))#Throw away the first .2sec to guarantee the IIR transient settles
if1 = freq1[0:int(Fs*sec_per_rev)]
if2 = freq1[int(Fs*sec_per_rev):int(2*Fs*sec_per_rev)]

maxf = (max(max(if1),max(if2))+.2) #This shuld be changed to make the maximum frequency an even 10th of a Hz.

r1 = 20.-(maxf-if1)/Hz_per_tick  #20 radial ticks at Hz_per_tick is fixed, adaptive scaling
r2 = 20.-(maxf-if2)/Hz_per_tick  #is an exercise for later

#plt.figure(1)
plt.figure(figsize=(11,11))
ax = plt.subplot(111, projection='polar')
ax.plot(theta,r1)
ax.plot(theta,r2)

dgr = (2*np.pi)/360.

mod_date = datetime.datetime.fromtimestamp(os.path.getmtime(_FILE))
ax.text(226.*dgr, 28.5, 'Mean Rev1 {:4.3f}Hz'.format(np.mean(if1)) + "\n" + \
'Mean Rev2 {:4.3f}Hz'.format(np.mean(if2)) + "\n" + \
_FILE + "\n" + \
mod_date.strftime("%b %d, %Y %H:%M:%S"), fontsize=9)

ax.set_rmax(20)
#Set up the ticks y is radial x is theta, it turns out x and y
#methods still work in polar projection but sometimes do funny things

tick_loc = np.arange(1,21,1)

myticks = []
for x in range(0,20,1):
myticks.append('{:4.2f}Hz'.format(maxf-(19*Hz_per_tick)+x*Hz_per_tick))

ax.set_rgrids(tick_loc, labels = myticks, angle = 90, fontsize = 8)

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(info_line, va='bottom', fontsize=16)

plt.savefig(info_line.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.5)

plt.show()``````

Python installation instructions are in this thread.

Last edited:
OP

#### JP

##### Major Contributor
Below is a "beta" version that includes the ability to select how many revolutions you want to plot, select an offset, and perform a standard x, y plot instead of a polar projection. The impetus for these features was another project I was working on where I borrowed a lot of the functionality of the original script, so I brought features from that project back in to the polar plot script. This is very much a work in progress, and one I don't give much time to.

Python:
``````'''
Hz_per_tick =     How many Hz per radial tick in the polar projection.

sec_per_rev =     How many second per one revolution of the platter:
33PRM = 1.8
45RPM = 1.33
78RPM = .766

revs =            The number of revolutions to plot.

revoffset =       Offset number of revolutions before plotting.

pltpolar =        0 - Plot x,y projection.
1 - Plot polar projection.

pltlegend =       0 - No legend.
1 - Legend with mean frequency per rev.

stereo_channel =  0 - Left Channel.
1 - Right Channel.

filter_freq =     Lowpass filter frequency in Hz.  Default is 60.

**** Version 6 Beta ****

'''

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 adjust input parameters:

_FILE = 'wavefile.wav'

info_line = 'PLOT TITLE'

Hz_per_tick = 3
sec_per_rev = 1.8
revs = 6
revoffset = 0

pltpolar = 1
pltlegend = 0

stereo_channel = 0
filter_freq = 60

#end edit

def instfreq(sig,Fs,filter_freq):
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,filter_freq/(Fs/2), btype='lowpass')
zi = signal.lfilter_zi(b, a) #Initialize the filter to the mean of the leading edge of the data
rawfreq,_ = signal.lfilter(b,a,rawfreq,zi=zi*np.mean(rawfreq[0:2000])) #reduces glitch, first pass

b, a = signal.iirfilter(3,filter_freq/(Fs/2), btype='lowpass')
instfreq = signal.filtfilt(b,a,rawfreq) #back and forth linear phase IIR filter (6 pole)

return (instfreq)

Fs = float(y[0])
if np.size(y[1][0]) == 2:
sig = y[1][:,stereo_channel][0:int(Fs*(sec_per_rev*(.5+revs+revoffset)))] #Grab revs + revoffest + .5 of audio from the specified channel
else:
sig = y[1][0:int(Fs*(sec_per_rev*3))] #mono file

freq1 = instfreq(sig,Fs,filter_freq)

freq1 = np.roll(freq1,-int(Fs*.2))#Throw away the first .2sec to guarantee the IIR transient settles

maxf = (max(freq1[int((Fs*sec_per_rev)*(revoffset)):int((Fs*sec_per_rev)*(revoffset+revs))])+1)

plotdata = {}
for x in range(revs):
plotdata[x] = freq1[int((Fs*sec_per_rev)*(x+revoffset)):int((Fs*sec_per_rev)*(x+1+revoffset))]
if pltpolar == 1:
plotdata[x] = 20.-(maxf-plotdata[x])/Hz_per_tick

if pltpolar == 1:

plt.figure(figsize=(11,11))
ax = plt.subplot(111, projection='polar')

t = np.arange(sec_per_rev,0,-1/Fs)  #Reverse time (theta axis)
theta = t*2*np.pi/sec_per_rev   #Time becomes degrees (1 rev = 2pi radians)
theta = np.roll(theta,int((sec_per_rev*Fs/4)))  #Rotate 90 deg to put 0 on top (sec_per_rev*Fs/4)

for key in plotdata:
ax.plot(theta, plotdata[key], label = 'Rev ' + str(key+1) + ': {:4.3f}Hz'.format(np.mean(plotdata[key])))

dgr = (2*np.pi)/360.

mod_date = datetime.datetime.fromtimestamp(os.path.getmtime(_FILE))
ax.text(226.*dgr, 28.5, _FILE + "\n" + \
mod_date.strftime("%b %d, %Y %H:%M:%S"), fontsize=9)

ax.set_rmax(20)
#Set up the ticks y is radial x is theta, it turns out x and y
#methods still work in polar projection but sometimes do funny things

tick_loc = np.arange(1,21,1)

myticks = []
for x in range(0,20,1):
myticks.append('{:4.2f}Hz'.format(maxf-(19*Hz_per_tick)+x*Hz_per_tick))

ax.set_rgrids(tick_loc, labels = myticks, angle = 90, fontsize = 8)

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}'])

else:

plt.figure(figsize=(20,8))
ax = plt.subplot(111)

t = np.arange(0,sec_per_rev,1/Fs)  #normal time (x axis)

for key in plotdata:
ax.plot(t, plotdata[key], label = 'Rev ' + str(key+1) + ': {:4.3f}Hz'.format(np.mean(plotdata[key])))

if pltlegend == 1:
ax.legend(loc='lower center', ncol=10)

ax.set_ylabel("Frequency (Hz)")
ax.set_xlabel("Time (s)")

if pltlegend == 1:
ax.legend(loc='lower center', ncol=10)

ax.grid(True)

ax.set_title(info_line, va='bottom', fontsize=16)

plt.savefig(info_line.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.5)

plt.show()``````

Last edited:
OP

Reserved

OP

Reserved

OP

Reserved

Forum Donor
Uh, why?

#### thegeton

##### Addicted to Fun and Learning
Forum Donor

Why did OP reserve 3 sequential post entries, or why do a TT W&F polar plot?

#### ban25

##### Addicted to Fun and Learning
Interesting. I have a belt-drive table that has pretty severe W+F and another direct-drive table that is flawless. Will have to take a measurement this weekend to compare.

OP

#### JP

##### Major Contributor
There was no audio data for that revolution - file isn’t long enough.

#### Balle Clorin

##### Major Contributor
Clearaudio has a test record with a 15 minute 3150Hz track.. here is 100 revolutions on that on a belt drive, yes it had drift that day... DC motor no tachometer control.

I found my error I have to have L&R file, single track does not work

OP

#### JP

##### Major Contributor
It should detect mono - that part hasn’t changed for years.

#### morillon

##### Major Contributor
Clearaudio has a test record with a 15 minute 3150Hz track.. here is 100 revolutions on that on a belt drive, yes it had drift that day... DC motor no tachometer control.

View attachment 327712

I found my error I have to have L&R file, single track does not work
?
have you identified the problem on the record or turntable side?

Last edited:

#### Balle Clorin

##### Major Contributor
Wow and Flutter. Cannot remember , old file. The WF track is not very good, + I may have adjusted the speed too.

Here is a new recording today 30 revolutions B&K QR2010 on Michell Gyro SE, belt DC motor

The best I have seen so far,
it even beats the w&f in this video for Technics top TT with the DS Audio ES-001 Eccentricity device…

I am not saying my TT is better than the Technics, but the test record limits the accuracy. Test records can only be used very carefully or else the results are misleading

Last edited:

#### mfalcon

##### Member
@JP I'm getting the below error
/Volumes/NAS/plots/wow and flutter/polarplot.py:96: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.

ax.set_xticklabels(['90'+u'\N{DEGREE SIGN}','45'+u'\N{DEGREE SIGN}','0'+u'\N{DEGREE SIGN}',\
I'm using Python 3.11.7.1. The OS is MacOS 10.13. This is my regular Python installation.

It's weird, I have like a minute of audio but it plots two revolutions, and locks up after throwing this warning. Then I have to control-c to get out of the terminal. Have you seen anything like this? Thanks in advance.

OP

#### JP

##### Major Contributor
That's just a warning - benign. If you're using the non-beta code it'll only plot two revolutions.

#### mfalcon

##### Member
@JP I started using the beta. Still get the warning but also good plots. Now if I can find a good test record for W&F.

thanks for all your work

OP

#### JP

##### Major Contributor
There's no really good test records, and different pressings means there aren't even known good matrix numbers necessarily. One method you can use to separate what is the record from the 'table is to zoom in to something in the waveform that can "index" the start of the track, and two a second plot with the record rotated 180deg on the platter.

#### n3mmr

##### Member
It seems the script assumes one complete turn takes precisely (60/33⅓) = 1.8 seconds.
Isn't that a bit trusting?
If you get enough data, and you know the nominal frequency you should be able to calculate a more exact real time for one rev, or?
Or did I miss something in the script?

OP

#### JP

##### Major Contributor
The nominal f being bang-on is rather trusting as well.

#### Balle Clorin

##### Major Contributor
The nominal frequency is off too sometimes, or most of the time ..

Replies
11
Views
1K
Replies
1K
Views
134K
Replies
7
Views
1K
Replies
20
Views
2K
Replies
13
Views
2K