• WANTED: Happy members who like to discuss audio and other topics related to our interest. Desire to learn and share knowledge of science required. There are daily reviews of audio hardware and expert members to help answer your questions. Click here to have your audio equipment measured for free!

Turntable Wow and Flutter "Polar" Plot Script

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
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)


y = read(_FILE)
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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
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)


y = read(_FILE)
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
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
Reserved
 
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
Reserved
 
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
Reserved
 

ban25

Addicted to Fun and Learning
Joined
Nov 5, 2022
Messages
693
Likes
672
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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
There was no audio data for that revolution - file isn’t long enough.
 

Balle Clorin

Major Contributor
Joined
Dec 26, 2017
Messages
1,237
Likes
1,137
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.


1700375410057.png



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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
It should detect mono - that part hasn’t changed for years.
 

morillon

Major Contributor
Joined
Apr 19, 2022
Messages
1,277
Likes
251
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
? :eek:
have you identified the problem on the record or turntable side?
 
Last edited:

Balle Clorin

Major Contributor
Joined
Dec 26, 2017
Messages
1,237
Likes
1,137
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
1700409308984.png


1700409350119.png

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
Joined
Jan 17, 2024
Messages
39
Likes
14
@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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
That's just a warning - benign. If you're using the non-beta code it'll only plot two revolutions.
 

mfalcon

Member
Joined
Jan 17, 2024
Messages
39
Likes
14
@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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
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
Joined
May 31, 2023
Messages
28
Likes
0
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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,250
Likes
2,433
Location
Brookfield, CT
The nominal f being bang-on is rather trusting as well.
 

Balle Clorin

Major Contributor
Joined
Dec 26, 2017
Messages
1,237
Likes
1,137
The nominal frequency is off too sometimes, or most of the time ..
 
Top Bottom