• 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!

Phono Cartridge Response Measurement Script

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
This thread is intended to consolidate the technical discussion and further development of Scott Wurcer's frequency response measurement script to a central area. I'd like to keep this discussion about the tool rather than specific cartridges.

I'll edit this and the subsequent reserved posts in the coming weeks with relevant info, background, links, etc.

Measurement Library thread is here.
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
The latest code for the script can be found here. I intend to get this on GitHub at some point, but as I'm not a developer I've only read code from repos and have never had to use them for workflow, so I've a learning curve. Patience appreciated.

You will need Python to run this with scipy, matplotlib, numpy, and librosa libraries installed. Note that I work on a Mac so there may be some methods and nomenclatures that are OS dependent, particularly around file handling. In the future I'll test against Windows as well to ensure it is agnostic.

Also note that indents are important with Python, so you want to ensure the formatting doesn't change when you copy/paste.

Python:
'''
PLEASE READ
To use this script you need to edit the user parameters section. 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.

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".

file_0 =        The first (L) or only file to plot

file_1 =        The second (R) file to plot.  If only using one file change to = ''

plotstyle =     1 - traditional
                2 - dual axis (twinx)
                3 - dual plot FR and distortion
                4 - dual plot FR zoom and dual axis (twinx)
                5 - small plot FR only

riaamode =      0 - off
                1 - bass emphasis
                2 - trebel de-emphasis
                3 - both

riaainv =       0 - disable
                1 - inverse RIAA EQ per riaamode setting

str100 =        0 - disable
                1 - enable 6dB/oct correction from 500Hz to 40Hz

xg7001 =        0 - disable
                1 - custom bass filter for Denon XG-7001 sweep (@stereoplay filter)

normalize =     Frequency in Hz to set as 0dB in the plot

file0norm =     0 - normalize both files independently
                1 - normalize both files to file 0 level

onekfstart =    0 - disable
                1 - start plot from 1kHz

endf =          0 - highest frequency to plot
               

'''

swversion = "17.0"



from scipy import signal
from scipy.io.wavfile import read
from pathlib import Path
from matplotlib.legend_handler import HandlerBase
from matplotlib.offsetbox import AnchoredText
from itertools import chain
from datetime import datetime
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
import os
import librosa




#edit user parameters


file_0 = 'P100CMK4_105pF_TRS-1007_L.wav'
file_1 = 'P100CMK4_105pF_TRS-1007_R.wav'

infoline = 'Cart / Load / Record'

equipinfo = 'Arm -> Phonostage -> ADC'

plotstyle = 4
plotdataout = 0
roundlvl = 1

riaamode = 1
riaainv =  0
str100 = 0
xg7001 = 0

normalize = 1000
file0norm = 0
onekfstart = 0
endf = 20000

ovdylim = 0
ovdylimvalue = [-35,5]

topdb = 100
framelength = 1024
hoplength = 256


#end Edit



fileopenidx = 0


def align_yaxis(ax1, ax2):
    y_lims = np.array([ax.get_ylim() for ax in [ax1, ax2]])

    # force 0 to appear on both axes, comment if don't need
    y_lims[:, 0] = y_lims[:, 0].clip(None, 0)
    y_lims[:, 1] = y_lims[:, 1].clip(0, None)

    # normalize both axes
    y_mags = (y_lims[:,1] - y_lims[:,0]).reshape(len(y_lims),1)
    y_lims_normalized = y_lims / y_mags

    # find combined range
    y_new_lims_normalized = np.array([np.min(y_lims_normalized), np.max(y_lims_normalized)])

    # denormalize combined range to get new axes
    new_lim1, new_lim2 = y_new_lims_normalized * y_mags
    return new_lim1, new_lim2


class AnyObjectHandler(HandlerBase):
    def create_artists(self, legend, orig_handle,
                       x0, y0, width, height, fontsize, trans):

        l1 = plt.Line2D([x0,y0+width], [0.7*height,0.7*height],
                           color=orig_handle[0], linestyle=orig_handle[1])

        l2 = plt.Line2D([x0,y0+width], [0.3*height,0.3*height],
                           color=orig_handle[2], linestyle=orig_handle[3])
        return [l1, l2]
 

def ft_window(n):       #Matlab's flat top window
    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



def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx


def createplotdata(insig, Fs):
    fout = []
    aout = []
    foutx = []
    aoutx = []
    fout2 = []
    aout2 = []
    fout3 = []
    aout3 = []
    global norm


    def interpolate(f, a, minf, maxf, fstep):
        f_out = []
        a_out = []
        amp = 0
        count = 0
        for x in range(minf,(maxf)+1,fstep):
            for y in range(0,len(f)):
                if f[y] == x:
                    amp = amp + a[y]
                    count = count + 1
            if count != 0:
                f_out.append(x)
                a_out.append(20*np.log10(amp/count))
            amp = 0
            count = 0
        return f_out, a_out

 
    def rfft(insig, Fs, minf, maxf, fstep):
        freq = []
        amp = []
        freqx = []
        ampx = []
        freq2h = []
        amp2h = []
        freq3h = []
        amp3h = []

        F = int(Fs/fstep)
        win = ft_window(F)

        if chinfile == 1:
            for x in range(0,len(insig)-F,F):
                y = abs(np.fft.rfft(insig[x:x+F]*win))
                f = np.argmax(y) #use largest bin
                if f >=minf/fstep and f <=maxf/fstep:
                    freq.append(f*fstep)
                    amp.append(y[f])
                if 2*f<F/2-2 and f > minf/fstep and f < maxf/fstep:
                    f2 = np.argmax(y[(2*f)-2:(2*f)+2])
                    freq2h.append(f*fstep)
                    amp2h.append(y[2*f-2+f2])
                if 3*f<F/2-2 and f > minf/fstep and f < maxf/fstep:
                    f3 = np.argmax(y[(3*f)-2:(3*f)+2])
                    freq3h.append(f*fstep)
                    amp3h.append(y[3*f-2+f3])


        else:
            for x in range(0,len(insig[0])-F,F):
                y0 = abs(np.fft.rfft(insig[0,x:x+F]*win))
                y1 = abs(np.fft.rfft(insig[1,x:x+F]*win))
                f0 = np.argmax(y0) #use largest bin
                f1 = np.argmax(y1) #use largest bin
                if f0 >=minf/fstep and f0 <=maxf/fstep:
                    freq.append(f0*fstep)
                    freqx.append(f1*fstep)
                    amp.append(y0[f0])
                    ampx.append(y1[f1])
                if 2*f0<F/2-2 and f0 > minf/fstep and f0 < maxf/fstep:
                    f2 = np.argmax(y0[(2*f0)-2:(2*f0)+2])
                    freq2h.append(f0*fstep)
                    amp2h.append(y0[2*f0-2+f2])
                if 3*f0<F/2-2 and f0 > minf/fstep and f0 < maxf/fstep:
                    f3 = np.argmax(y0[(3*f0)-2:(3*f0)+2])
                    freq3h.append(f0*fstep)
                    amp3h.append(y0[3*f0-2+f3])

        return freq, amp, freqx, ampx, freq2h, amp2h, freq3h, amp3h

 

    def normstr100(f, a):
        fmin = 40
        fmax = 500
        slope = -6.02
        for x in range(find_nearest(f, fmin), (find_nearest(f, fmax))):
            a[x] = a[x] + 20*np.log10(1*((f[x])/fmax)**((slope/20)/np.log10(2)))
        return a


    def chunk(insig, Fs, fmin, fmax, step, offset):
        f, a, fx, ax, f2, a2, f3, a3 = rfft(insig, Fs, fmin, fmax, step)
        f, a = interpolate(f, a, fmin, fmax, step)
        fx, ax = interpolate(fx, ax, fmin, fmax, step)
        f2, a2 = interpolate(f2, a2, fmin, fmax, step)
        f3, a3 = interpolate(f3, a3, fmin, fmax, step)
        a = [x - offset for x in a]
        ax = [x - offset for x in ax]
        a2 = [x - offset for x in a2]
        a3 = [x - offset for x in a3]

        return f, a, fx, ax, f2, a2, f3, a3
 

    def concat(f, a, fx, ax, f2, a2, f3, a3, fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3):
        fout = fout + f
        aout = aout + a
        foutx = foutx + fx
        aoutx = aoutx + ax
        fout2 = fout2 + f2
        aout2 = aout2 + a2
        fout3 = fout3 + f3
        aout3 = aout3 + a3

        return fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3
 
    if onekfstart == 0:


        f, a, fx, ax, f2, a2, f3, a3 = chunk(insig, Fs, 20, 45, 5, 26.03)
        fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3 = concat(f, a, fx, ax, f2, a2, f3, a3, fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3)

        f, a, fx, ax, f2, a2, f3, a3 = chunk(insig, Fs, 50, 90, 10, 19.995)
        fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3 = concat(f, a, fx, ax, f2, a2, f3, a3, fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3)

        f, a, fx, ax, f2, a2, f3, a3 = chunk(insig, Fs, 100, 980, 20, 13.99)
        fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3 = concat(f, a, fx, ax, f2, a2, f3, a3, fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3)

    f, a, fx, ax, f2, a2, f3, a3 = chunk(insig, Fs, 1000, endf, 100, 0)
    fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3 = concat(f, a, fx, ax, f2, a2, f3, a3, fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3)

 
    if str100 == 1:
        aout = normstr100(fout, aout)
        aout2 = normstr100(fout2, aout2)
        aout3 = normstr100(fout3, aout3)
        if chinfile == 2:
            aoutx = normstr100(foutx, aoutx)


    if file0norm == 1 and fileopenidx == 1:
        i = find_nearest(fout, normalize)
        norm = aout[i]
    elif file0norm == 0:
        i = find_nearest(fout, normalize)
        norm = aout[i]


    aout = aout-norm #amplitude is in dB so normalize by subtraction at [i]
    aoutx = aoutx-norm
    aout2 = aout2-norm
    aout3 = aout3-norm
 
    sos = signal.iirfilter(3,.5, btype='lowpass', output='sos') #filter some noise
    aout = signal.sosfiltfilt(sos,aout)
    aout2 = signal.sosfiltfilt(sos,aout2)
    aout3 = signal.sosfiltfilt(sos,aout3)

    if chinfile == 2 and len(aoutx) >1:
        aoutx = signal.sosfiltfilt(sos,aoutx)

    return fout, aout, foutx, aoutx, fout2, aout2, fout3, aout3

 

def ordersignal(sig, Fs):
    F = int(Fs/100)
    win = ft_window(F)

    if chinfile == 1:
        y = abs(np.fft.rfft(sig[0:F]*win))
        minf = np.argmax(y)
        y = abs(np.fft.rfft(sig[len(sig)-F:len(sig)]*win))
        maxf = np.argmax(y)
    else:
        y = abs(np.fft.rfft(sig[0,0:F]*win))
        minf = np.argmax(y)
        y = abs(np.fft.rfft(sig[0][len(sig[0])-F:len(sig[0])]*win))
        maxf = np.argmax(y)

    if maxf < minf:
        maxf,minf = minf,maxf
        sig = np.flipud(sig)


 
    return sig, minf, maxf



def riaaiir(sig, Fs, mode, inv):
    if Fs == 96000:
        at = [1, -0.66168391, -0.18158841]
        bt = [0.1254979638905360, 0.0458786797031512, 0.0018820452752401]
        ars = [1, -0.60450091, -0.39094593]
        brs = [0.90861261463964900, -0.52293147388301200, -0.34491369168550900]
    if inv == 1:
        at,bt = bt,at
        ars,brs = brs,ars
    if mode == 1:
        sig = signal.lfilter(brs,ars,sig)
    if mode == 2:
        sig = signal.lfilter(bt,at,sig)
    if mode == 3:
        sig = signal.lfilter(bt,at,sig)
        sig = signal.lfilter(brs,ars,sig)
    return sig


def normxg7001(sig, Fs):
    if Fs == 96000:
        b = [1.0080900, -0.9917285, 0]
        a = [1, -0.9998364, 0]
        sig = signal.lfilter(b,a,sig)
    return sig


def openaudio(_FILE):
    global chinfile
    global fileopenidx
    chinfile = 1

    srinfile = librosa.get_samplerate(_FILE)
 
    audio, Fs = librosa.load(_FILE, sr=None, mono=False)


    if len(audio.shape) == 2:
        chinfile = 2
        filelength = audio.shape[1] / Fs
    else:
        filelength = audio.shape[0] / Fs

    print('Input File:   ' + str(_FILE))
    print('Sample Rate:  ' + str("{:,}".format(srinfile) + 'Hz'))

    if Fs <96000:
        print('              Resampling to 96,000Hz')
        audio = librosa.resample(audio, orig_sr=Fs, target_sr=96000)
        Fs = 96000
 
    print('Channels:     ' + str(chinfile))
    print(f"Length:       {filelength}s")


    if riaamode != 0:
        audio = riaaiir(audio, Fs, riaamode, riaainv)

    if xg7001 == 1:
        audio = normxg7001(audio, Fs)


    audio, index = librosa.effects.trim(audio, top_db=topdb, frame_length=framelength, hop_length=hoplength)
 
    print(f"In/Out (s):   {index / Fs}")


    audio, minf, maxf = ordersignal(audio, Fs)

    print('Min Freq:     ' + str("{:,}".format(minf * 100) + 'Hz'))
    print('Max Freq:     ' + str("{:,}".format(maxf * 100) + 'Hz\n'))
 
    fileopenidx +=1
 
    return audio, Fs, minf, maxf






if __name__ == "__main__":


    input_sig, Fs, minf, maxf = openaudio(file_0)
    fo0, ao0, fox0, aox0, fo2h0, ao2h0, fo3h0, ao3h0 = createplotdata(input_sig, Fs)

    deltaadj = ao0[find_nearest(fo0, normalize)]
    deltah0 = round((max(ao0 - deltaadj)), roundlvl)
    deltal0 = abs(round((min(ao0 - deltaadj)), roundlvl))

    if aox0.size > 0:
        print('X-talk @1kHz: ' + (str(round(aox0[find_nearest(fox0, 1000)], 2))) + 'dB\n\n')
   

    if file_1:
        input_sig, Fs, minf, maxf = openaudio(file_1)
        fo1, ao1, fox1, aox1, fo2h1, ao2h1, fo3h1, ao3h1 = createplotdata(input_sig, Fs)
 
        deltaadj = ao1[find_nearest(fo1, normalize)]
        deltah1 = round((max(ao1 - deltaadj)), roundlvl)
        deltal1 = abs(round((min(ao1 - deltaadj)), roundlvl))

        if aox1.size > 0:
            print('X-talk @1kHz: ' + (str(round(aox1[find_nearest(fox1, 1000)], 2))) + 'dB\n\n')



    if plotdataout == 1:

        dao0 = [*ao0, *[''] * (len(fo0) - len(ao0))]
        daox0 = [*aox0, *[''] * (len(fo0) - len(aox0))]
        dao2h0 = [*ao2h0, *[''] * (len(fo0) - len(ao2h0))]
        dao3h0 = [*ao3h0, *[''] * (len(fo0) - len(ao3h0))]

        print('\n\nFile 0 Plot Data: (freq, ampl, x-talk, 2h, 3h)\n\n')

        dataout = list(zip(fo0, dao0, daox0, dao2h0, dao3h0))
        for fo, ao, aox, ao2, ao3 in dataout:
            print(fo, ao, aox, ao2, ao3, sep=', ')

        if file_1:
            dao1 = [*ao1, *[''] * (len(fo1) - len(ao1))]
            daox1 = [*aox1, *[''] * (len(fo1) - len(aox1))]
            dao2h1 = [*ao2h1, *[''] * (len(fo1) - len(ao2h1))]
            dao3h1 = [*ao3h1, *[''] * (len(fo1) - len(ao3h1))]

            print('\n\nFile 1 Plot Data: (freq, ampl, x-talk, 2h, 3h)\n\n')

            dataout = list(zip(fo1, dao1, daox1, dao2h1, dao3h1))
            for fo, ao, aox, ao2, ao3 in dataout:
                print(fo, ao, aox, ao2, ao3, sep=', ')



    plt.rcParams["xtick.minor.visible"] =  True
    plt.rcParams["ytick.minor.visible"] =  True

    if plotstyle == 1:
        fig, axs = plt.subplots(1, 1, figsize=(14,6))
        axs = np.ravel([axs])


        axs[0].semilogx(fo0,ao0, color = '#0000ff', label = 'Freq Response')

        axs[0].semilogx(fo2h0,ao2h0,color = '#0080ff', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
        axs[0].semilogx(fo3h0,ao3h0,color = '#00dfff', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

        axs[0].semilogx(fox0,aox0,color = '#0000ff', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')
 


        if file_1:
            axs[0].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')

            axs[0].semilogx(fo2h1,ao2h1,color = '#ff8000', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
            axs[0].semilogx(fo3h1,ao3h1,color = '#ffdf00', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

            axs[0].semilogx(fox1,aox1,color = '#ff0000', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')

            plt.legend([("#0000ff", "-", "#ff0000", "-"), ("#0000ff", (0, (3, 1, 1, 1)), "#ff0000", (0, (3, 1, 1, 1))),
                        ("#0080ff", "-", "#ff8000", "-"), ("#00dfff", "-", "#ffdf00", "-")],
                       ['Freq Response', 'Crosstalk', '2ⁿᵈ Harmonic', '3ʳᵈ Harmonic'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)

            axs[0].set_ylim((min(chain(aox0, aox1)) -2), (max(chain(ao0, ao1)) +2))

        else:
     
            plt.legend(loc=4)

        axs[0].set_ylabel("Amplitude (dB)")
        axs[0].set_xlabel("Frequency (Hz)")

        plt.autoscale(enable=True, axis='y')

        if ovdylim == 1:
            axs[0].set_ylim(*ovdylimvalue)
 


    if plotstyle == 2:
        fig, axs = plt.subplots(1, sharex=True, figsize=(14,6))
        axs = np.ravel([axs])
        axtwin = axs[0].twinx()


        axs[0].set_ylim(-5,5)



        if max(ao0) <7:
            axs[0].set_ylim(-25, 7)

        if max(ao0) < 4:
            axs[0].set_ylim(-25,5)
     
        if max(ao0) < 2:
            axs[0].set_ylim(-29,3)

        if max(ao0) < 0.5:
            axs[0].set_ylim(-30,2)


        if aox0.size > 0:
            if file_1:
                axs[0].set_ylim((min(chain(aox0, aox1)) -2), (max(chain(ao0, ao1)) +2))
            else:
                axs[0].set_ylim((min(aox0) -2), (max(ao0) +2))
 
 
        if ovdylim == 1:
            axs[0].set_ylim(*ovdylimvalue)

 

        axs[0].semilogx(fo0,ao0, color = '#0000ff', label = 'Freq Response')

        axtwin.semilogx(fo2h0,ao2h0,color = '#0080ff', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
        axtwin.semilogx(fo3h0,ao3h0,color = '#00dfff', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

        axs[0].semilogx(fox0,aox0,color = '#0000ff', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')
 


        if file_1:
            axs[0].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')

            axtwin.semilogx(fo2h1,ao2h1,color = '#ff8000', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
            axtwin.semilogx(fo3h1,ao3h1,color = '#ffdf00', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

            axs[0].semilogx(fox1,aox1,color = '#ff0000', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')

            plt.legend([("#0000ff", "-", "#ff0000", "-"), ("#0000ff", (0, (3, 1, 1, 1)), "#ff0000", (0, (3, 1, 1, 1))),
                        ("#0080ff", "-", "#ff8000", "-"), ("#00dfff", "-", "#ffdf00", "-")],
                       ['Freq Response', 'Crosstalk', '2ⁿᵈ Harmonic', '3ʳᵈ Harmonic'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)


            if aox0.size > 0  and aox1.size > 0:
                axs[0].set_ylim((min(chain(aox0, aox1)) -2), (max(chain(ao0, ao1)) +2))

        else:
     
            lines1, labels1 = axs[0].get_legend_handles_labels()
            lines2, labels2 = axtwin.get_legend_handles_labels()
            plt.legend(lines1 + lines2, labels1 + labels2, loc=4)
     

        new_lim1, new_lim2 = align_yaxis(axs[0], axtwin)
        axs[0].set_ylim(new_lim1)
        axtwin.set_ylim(new_lim2)


        axs[0].set_ylabel("Amplitude (dB)")
        axtwin.set_ylabel("Distortion (dB)")
        axs[0].set_xlabel("Frequency (Hz)")



    if plotstyle == 3:
        fig, axs = plt.subplots(2, 1, sharex=True, figsize=(14,6))


        axs[0].set_ylim(-5,5)


        if file_1:
            if (min(chain(ao0, ao1)) <-5) or (max(chain(ao0, ao1)) >5):
                axs[0].autoscale(enable=True, axis='y')
        elif (min(ao0) <-5) or (max(ao0) >5):
                axs[0].autoscale(enable=True, axis='y')

        if ovdylim == 1:
            axs[0].set_ylim(*ovdylimvalue)


        axs[0].semilogx(fo0,ao0,color = '#0000ff', label = 'Freq Response')
        axs[1].semilogx(fo2h0,ao2h0,color = '#0080ff', label = '2nd Harmonic')
        axs[1].semilogx(fo3h0,ao3h0,color = '#00dfff', label = '3rd Harmonic')


        if file_1:
            axs[0].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')

            axs[1].semilogx(fo2h1,ao2h1,color = '#ff8000', label = '2ⁿᵈ Harmonic')
            axs[1].semilogx(fo3h1,ao3h1,color = '#ffdf00', label = '3ʳᵈ Harmonic')


            axs[0].legend([("#0000ff", "-", "#ff0000", "-"),],
                       ['Freq Response'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)
     
            axs[1].legend([("#0080ff", "-", "#ff8000", "-"), ("#00dfff", "-", "#ffdf00", "-")],
                       ['2ⁿᵈ Harmonic', '3ʳᵈ Harmonic'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)


        else:
            axs[0].legend(loc=4)
            axs[1].legend(loc=4)


        axs[0].set_ylabel("Amplitude (dB)")
        axs[1].set_ylabel("Distortion (dB)")
        axs[1].set_xlabel("Frequency (Hz)")

        

    if plotstyle == 4:
        fig, axs = plt.subplots(2, 1, sharex=True, figsize=(14,10))
        axtwin = axs[1].twinx()


        axs[0].set_ylim(-5,5)


        if max(ao0) <7:
            axs[1].set_ylim(-25, 7)

        if max(ao0) < 4:
            axs[1].set_ylim(-25,5)
     
        if max(ao0) < 2:
            axs[1].set_ylim(-29,3)

        if max(ao0) < 0.5:
            axs[1].set_ylim(-30,2)


        if aox0.size > 0:
            if file_1:
                axs[1].set_ylim((min(chain(aox0, aox1)) -2), (max(chain(ao0, ao1)) +2))
            else:
                axs[1].set_ylim((min(aox0) -2), (max(ao0) +2))
 
 
        if ovdylim == 1:
            axs[1].set_ylim(*ovdylimvalue)

 

        axs[1].semilogx(fo0,ao0, color = '#0000ff', label = 'Freq Response')

        axtwin.semilogx(fo2h0,ao2h0,color = '#0080ff', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
        axtwin.semilogx(fo3h0,ao3h0,color = '#00dfff', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

        axs[1].semilogx(fox0,aox0,color = '#0000ff', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')
 


        if file_1:
            axs[1].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')

            axtwin.semilogx(fo2h1,ao2h1,color = '#ff8000', label = '2ⁿᵈ Harmonic', alpha = 1, linewidth = 0.75)
            axtwin.semilogx(fo3h1,ao3h1,color = '#ffdf00', label = '3ʳᵈ Harmonic', alpha = 1, linewidth = 0.75)

            axs[1].semilogx(fox1,aox1,color = '#ff0000', linestyle = (0, (3, 1, 1, 1)), label = 'Crosstalk')

            plt.legend([("#0000ff", "-", "#ff0000", "-"), ("#0000ff", (0, (3, 1, 1, 1)), "#ff0000", (0, (3, 1, 1, 1))),
                        ("#0080ff", "-", "#ff8000", "-"), ("#00dfff", "-", "#ffdf00", "-")],
                       ['Freq Response', 'Crosstalk', '2ⁿᵈ Harmonic', '3ʳᵈ Harmonic'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)


            if aox0.size > 0  and aox1.size > 0:
                axs[1].set_ylim((min(chain(aox0, aox1)) -2), (max(chain(ao0, ao1)) +2))

        else:
     
            lines1, labels1 = axs[1].get_legend_handles_labels()
            lines2, labels2 = axtwin.get_legend_handles_labels()
            plt.legend(lines1 + lines2, labels1 + labels2, loc=4)
     

        new_lim1, new_lim2 = align_yaxis(axs[1], axtwin)
        axs[1].set_ylim(new_lim1)
        axtwin.set_ylim(new_lim2)



        if file_1:
            if (min(chain(ao0, ao1)) <-5) or (max(chain(ao0, ao1)) >5):
                axs[0].autoscale(enable=True, axis='y')
        elif (min(ao0) <-5) or (max(ao0) >5):
                axs[0].autoscale(enable=True, axis='y')

        if ovdylim == 1:
            axs[0].set_ylim(*ovdylimvalue)


        axs[0].semilogx(fo0,ao0,color = '#0000ff', label = 'Freq Response')
     


        if file_1:
            axs[0].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')

        


            axs[0].legend([("#0000ff", "-", "#ff0000", "-"),],
                       ['Freq Response'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)
     
           

        else:
            axs[0].legend(loc=4)     


        axs[0].set_ylabel("Amplitude (dB)")
        axs[1].set_ylabel("Amplitude (dB)")
        axtwin.set_ylabel("Distortion (dB)")
        axs[1].set_xlabel("Frequency (Hz)")

        gs = GridSpec(2, 1, height_ratios=[1, 2])
        axs[0].set_position(gs[0].get_position(fig))
        axs[1].set_position(gs[1].get_position(fig))




    if plotstyle == 5:
        fig, axs = plt.subplots(1, 1, figsize=(14,3))
        axs = np.ravel([axs])


        axs[0].set_ylim(-5,5)


        if file_1:
            if (min(chain(ao0, ao1)) <-5) or (max(chain(ao0, ao1)) >5):
                axs[0].autoscale(enable=True, axis='y')
        elif (min(ao0) <-5) or (max(ao0) >5):
                axs[0].autoscale(enable=True, axis='y')

        if ovdylim == 1:
            axs[0].set_ylim(*ovdylimvalue)


        axs[0].semilogx(fo0,ao0,color = '#0000ff', label = 'Freq Response')


        if file_1:
            axs[0].semilogx(fo1,ao1, color = '#ff0000', label = 'Freq Response')



            axs[0].legend([("#0000ff", "-", "#ff0000", "-"),],
                       ['Freq Response'],
                       handler_map={tuple: AnyObjectHandler()},loc=4)


        else:
            axs[0].legend(loc=4)
    

        axs[0].set_ylabel("Amplitude (dB)")
        axs[0].set_xlabel("Frequency (Hz)")




    for i, ax in enumerate(axs.flat):
 
        if file0norm == 1:
            if  not (plotstyle == 3 and i != 0):
                ax.axline((normalize, 0), (normalize, 1), color = 'm', lw = 1)
        elif file0norm == 0:
            if not (plotstyle == 3 and i != 0):
                ax.plot(normalize, 0, marker = 'x', color = 'm')
        
        anchored_text = AnchoredText('SJ', 
                            frameon=False, borderpad=0, pad=0.03, 
                            loc=1, bbox_transform=plt.gca().transAxes,
                            prop={'color':'m','fontsize':25,'alpha':.4,
                            'style':'oblique'})
        ax.add_artist(anchored_text)
        
        ax.grid(True, which="major", axis="both", ls="-", color="black")
        ax.grid(True, which="minor", axis="both", ls="-", color="gainsboro")

        ax.set_xticks([0,20,50,100,500,1000,5000,10000,20000,50000,100000])
        ax.set_xticklabels(['0','20','50','100','500','1k','5k','10k','20k','50k','100k'])



    bbox_args = dict(boxstyle="round", color='b', fc='w', ec='b', alpha=1, pad=.15)
    axs[0].annotate('+' + str(deltah0) + ', ' + u"\u2212" + str(deltal0) + ' dB',color = 'b',\
             xy=(fo0[0],(ao0[0]-1)), xycoords='data', \
             xytext=(-10, -20), textcoords='offset points', \
             ha="left", va="center", bbox=bbox_args)

    if file_1:
        bbox_args = dict(boxstyle="round", color='b', fc='w', ec='r', alpha=1, pad=.15)
        axs[0].annotate('+' + str(deltah1) + ', ' + u"\u2212" + str(deltal1) + ' dB',color = 'r',\
                 xy=(fo0[0],(ao0[0]-1)), xycoords='data', \
                 xytext=(-10, -34.5), textcoords='offset points', \
                 ha="left", va="center", bbox=bbox_args)

     


    plt.autoscale(enable=True, axis='x')

    axs[0].set_title(infoline + "\n", fontsize=16)


    now = datetime.now()

    if file_1:
        plt.figtext(.17, .118, "SJPlot v" + swversion + "\n" + file_0 + "\n" + file_1 + "\n" + \
            now.strftime("%b %d, %Y %H:%M"), fontsize=6)
    else:
        plt.figtext(.17, .118, "SJPlot v" + swversion + "\n" + file_0 + "\n" + \
            now.strftime("%b %d, %Y %H:%M"), fontsize=6)



    plt.figtext(.125, 0, equipinfo, alpha=.75, fontsize=8)

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

    plt.show()

    print('\nDone!')
-Version 17.0, Dec 17, 2023-


This is an example plot:

Version 16.3 Example.png


These instructions are for the official Python distribution, which is what I use. You can use whatever distribution you like (Brew, Anaconda, etc.) however I won’t be able to help you with any issues. On the other hand, I’m not a heavy Python user so my help would likely be of limited value anyway.

As of this writing I recommend Python 3.10.9 as 3.11.1 is not supported by Numba, which is a dependency for librosa. You should do your own research on which version you want to run as I won’t necessarily always keep this information current. The good news is you can simply remove and reinstall a different version if you run into package compatibility issues.

To get 3.10.9, go to https://python.org/downloads. Scroll to the section titled “Looking for a specific release?” and click the download link for Python 3.10.9.
Mac-1.png


On the next page scroll down to “Files” and choose the macOS 64-bit universal2 installer.
Mac-2.png


Double-click the installer package from your Downloads folder and follow the defaults to complete installation.
Mac-3.png
It’s a good idea to keep all the files in the same directory as the script so that you don’t have to worry about specifying paths. For our purposes here we’ll name this folder Plots and place it in the Documents folder. You can name it and put it wherever suits you.
Mac-4.png



Next we need to install a few libraries. To do this, open Terminal and run the following three commands:


Code:
pip3 install --upgrade pip

pip3 install librosa

pip3 install matplotlib


You’ll want to ensure these processes complete successfully. Installing librosa will install two other required libraries, NumPy and SciPy, so there’s no need to do these separately.

Next launch IDLE, and select File, New File.

Mac-5.png




Copy and paste the script code from this post in to the “untitled” window, and select File, Save As.

Mac-6.png




Name the script something that makes sense to you, such as FR14.py – FR for Frequency Response, and 14 as example of the version. It doesn’t really matter, but make the extension .py, and be sure to save it in the folder you created above, which is ~Documents/Plots in our example.

Mac-7.png




At this point everything should be set for your first plot. Details of the file requirements for what you want to measure will be at the top of the script file, as will be an explanation of all the settings.

Copy your wav file in to the Plots folder and specify its name in the _FILE variable of the script. Don’t forget to update the Infoline to describe your plot. Then simply select Run, Run Module. You’ll be prompted to save the modified script file, after which the script will run.

MAc-8.png



The next time you run IDLE, simply select File, Open, and select the saved script file.

Windows Installation Instructions.
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
  • Fix axis label for amplitude difference. For channels that aren't normalized to 0, the signs may be incorrect.
  • Add (I)RIAA filters for 192kHz.
  • Add filter for Denon XG-7001 (in 17 Beta).
  • Auto-detect channel order.
  • Do not show "crosstalk" in the legend for mono files.
  • Consider showing phase correlation between the fundamental and crosstalk.
  • Add the ability to process two files for plots sets of both left and right channels on the same figure.
  • Improve y-axis scaling logic.

17.0
  • Added two new plot styles: (4) a dual plot with FR zoom and dual-axis, and (5) a small FR-only plot
  • Added EQ filter for Denon XG-7001 (filter courtesy of @stereoplay)
  • Corrected +/- annotation to be relative to the normalization frequency for the respective file rather than file_1 being relative to file_0. IOW, this no longer factors channel balance in to file_1 calcs.
  • If file0norm = 0, there will be a magenta "X" at 0dB on the plot at the normalization frequency. This denotes that each file was normalized individually, and what the normalization frequency was.
  • If file0norm = 1, there will be a vertical magenta line at the normalization frequency to denote that file_1 was normalized to file_0 level at the normalization frequency, and what the normalization frequency was.
  • (Internal) refactored the plot axes to make them iterable.

16.5
  • Added version to plot image.

16.4
  • Fixed plot scaling error with mono files.

16.3 Beta
  • Removed pre-filter.

16.2 Beta
  • Fixed copy/paste error with plot style 1

16.1 Beta
  • Added equipment info line to describe the capture signal chain.
  • Crosstalk @1kHz outputs on the console line.
  • Adjusted font sizes.

16 Beta
  • Added the ability to process two files, each file representing a channel. File_0 is left, and if specified file_1 is right.
  • Added the ability to normalize each file independently, or to normalize both to file_0.
  • Fixed cosmetic and layout issues.
  • Plotdataout will output data for both files.
  • Date/Time on plot is now the date/time the plot was generated instead of the input file date/time.

15.3 Alpha
  • Added the ability for an Butterworth HP pre-filter. Variable prefilterfreq is the filter frequency and prefilterorder is the order of the filter. If prefilterfreq = 0 no filter will be applied. Be careful to ensure the filter isn't attenuating frequencies of interest.

15.2 Alpha
  • Let matplotlib decide the location for the legend box in plot style 2 instead of hard-coding.

15.1 Alpha
  • Fixed axis colors for plot style 1.

15 Alpha
  • Detects mono or stereo file. If stereo the script will assume the right channel is crosstalk and will plot it as such. If you want to plot R<L, you'll need to swap channels prior to plotting.
  • Script will SRC files <96kHz to 96kHz, but will not SRC files above 96kHz. The only (I)RIAA filters in the script presently are for 96kHz.
  • Variable onekfstart will star the plot at 1kHz if set to 1. Variable endf is the highest frequency to plot. The basement for this setting is around 3000.
  • Plot style 3 will not plot crosstalk.
14 Alpha - Internal only.

13 Beta
  • Full range plots.
  • Compensation for STR-100.
  • RIAA and inverse RIAA implemented separately for bass and treble.
  • SRC all files to 96kHz (temporary to make filtering easier).

Code:
 
Last edited:
OP
JP

JP

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

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
APPENDIX

Comparison plots of all STR-100 versions.

B&K QR2009 and QR2010 comparison plots.

Code:
[a0, a1, a2]
[b0, b1, b2]

Sample Rate = 96k

     10kHz Attenuation (-13.73dB)

          [1., -0.66168391, -0.18158841]
          [1., 0.3655731, 0.01499662]

          ^ Gain at 1kHz: 18.0272664dB

          Unity gain at 1kHz:

               Standard

                    [1., -0.66168391, -0.18158841]
                    [0.125497963890536, 0.0458786797031512, 0.0018820452752401]

                    (biquad-m *track* 1.2549796389053600E-01 4.5878679703151200E-02 1.8820452752400900E-03 1.0000000000000000E+00 -6.6168391000000000E-01 -1.8158841000000000E-01)

               Inverse

                    [1., 0.3655731, 0.01499662]
                    [7.96825676687663, 5.27246729339089, 1.44694307676887]
 
                    (biquad-m *track* 7.9682567668766300E+00 -5.2724672933908900E+00 -1.4469430767688700E+00 1.0000000000000000E+00 3.6557310000000000E-01 1.4996620000000000E-02)

 
      Roll-off and Shelf (50.05Hz and 500.5Hz)

           [1., -0.60450091, -0.39094593]
           [1., -0.57552742, -0.37960478]

           ^ Gain at 1kHz: 0.83242476dB

           Unity gain at 1kHz:

                Standard

                     [1., -0.60450091, -0.39094593]
                     [0.908612614639649, -0.522931473883012, -0.344913691685509]

                     (biquad-m *track* 9.08612614639649000E-01 -5.22931473883012000E-01 -3.44913691685509000E-01 1.00000000000000000E+00 -6.04500910000000000E-01 -3.90945930000000000E-01)

                Inverse

                     [1., -0.57552742, -0.37960478]
                     [1.1005790409333, -0.665301031771105, -0.430266896696176]

                     (biquad-m *track* 1.10057904093330000E+00 -6.65301031771105000E-01 -4.30266896696176000E-01 1.00000000000000000E+00 -5.75527420000000000E-01 -3.79604780000000000E-01)


To chain biquads in Audacity Nyquist prompt:

(setf *track* (biquad-m *track* b0 b1 b2 a0 a1 a2)) (biquad-m *track* b0 b1 b2 a0 a1 a2)

Chain in the order above - 10kHz attenuation, and then roll-off and shelf.
 
Last edited:

USER

Addicted to Fun and Learning
Forum Donor
Joined
Mar 30, 2019
Messages
865
Likes
1,357
How to run the script in Windows

3/26/23 NOTE: There may be an issue with the latest version of Anaconda released a few days ago. Avoid the 2023.03 release of Anaconda for now and use an earlier, proven version (the previous release was 2022.10). There is likely a compatibility issue involving the updated version of the soundfile library that comes with it.

Use the 2022.10 release located here: https://repo.anaconda.com/archive/


1. Download and install Anaconda from here: https://www.anaconda.com/products/distribution
  • Install with default, recommended settings
  • Latest versions of Anaconda come with most of the required libraries (e.g. mathplotlib, scipy, numpy)
    • However, you will need to manually install librosa
Untitled.jpg


2. Open Anaconda Navigator and install librosa library
  • Anaconda Navigator is in your apps list on the Start menu if you are using a recent version of Windows
    • Usually in a folder named "Anaconda3 (64-bit)"
  • Open CMD.exe Prompt from within Anaconda Navigator (Home) as seen in the image above and run the following command:
conda install -c conda-forge librosa

Untitled3.jpg

  • Remember, especially if you are only using python to run this script: if it ain't broke, don't fix it
    • Don't upgrade components for the sake of upgrading
      • You are just fine if you are using, say, Python 10.9 and Spyder 5.2.2
    • Updates can lead to library incompatibility issues (e.g. missing .dll error) and you may have to re-install everything
3. Run sub-program Spyder
  • You can open Spyder from within Anaconda Navigator (Home) but once all set up you can run Spyder directly as a standalone Windows app
Untitleds.jpg


4. Open the script in Spyder and customize with minor edits
  • If you have the script and sweep recording files in different folders, find and adjust the following variables used for left and right channel sweeps, respectively:
file_0 = 'wavefile1.wav'
file_1 = 'wavefile2.wav'
  • Use full path names
file_0 = ''/Users/USER/Desktop/Folder/wavefile1.wav'
file_1 = ''/Users/USER/Desktop/Folder/wavefile2.wav'
  • Provide detailed information about your set-up with the following variables:
infoline = 'Cart / Load / Record'
equipinfo = 'Arm -> Phonostage -> ADC'
  • infoline: to make the measurements useful to others please provide as much information as possible
    • Loading information (e.g. total capacitance and resistance) is particularly important
      • Accounting for cable capacitance is strongly encouraged as it affects the frequency response of MM cartridges
    • The DER EE DE-5000 LCR meter is highly recommended for measuring your cartridge channels as well as cable capacitance
      • It is an extremely useful tool with which to check if your cartridge meets spec or is functioning properly
      • Regular multimeters exaggerate cable capacitance
  • The following is an example of the information I would provide
infoline = 'Shure V15 V-MR | L & R | 1.25g (No Brush) | 47kΩ | 210pF | CBS STR-100'
equipinfo = 'Denon DP-35F - Wayne Kirkwood Flat MM - E1DA Cosmos ADC'
  • Adjust additional variables, including the following
riaamode
riaainv
str100
  • Your test record and recording characteristics (e.g. RIAA equalization) determine your settings
    • Please see the script itself for information regarding these settings and options
  • If you want output larger images, increase the dpi variable in the following line toward the end of the script (e.g. 96 to 300):
plt.savefig(infoline.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.75, dpi=96)
5. Run the script when ready
  • Press the green play button on top (F5)

-----------------------------------------

So you got cute and tried updating or changing something you shouldn't have, and uninstalling and reinstalling doesn't get you back.
The most sure-fire thing to do is a perform a full uninstall of Anaconda. Here's how to do it:
  • Open Anaconda Navigator and open CMD.exe Prompt
    • Then type the following:
conda install anaconda-clean
  • When that is finished running, type the following:
anaconda-clean --yes
  • This sets-up a thorough Anaconda uninstall
  • After this, uninstall Anaconda as you normally would
  • Restart your computer
  • Great. Now get rid of as many of the leftover files that are causing problems as possible
    • Go to your equivalent of these two folders:
C:\Users\NAME
C:\Users\NAME\AppData
  • In "C:\Users\NAME" delete anything associated with "anaconda" "conda" "spider" etc.
Untitled75.jpg

  • In this example, I deleted:
.conda
.ipython
.matplotlib
.spyder-py3
anaconda3
and the file .condarc
  • In "C:\Users\NAME\AppData" go to both "local" and "roaming"
    • You may see folders like the following:
Local: conda, spyder, pip
Roaming: jupyter
  • Delete those
  • Restart your computer
  • Install Anaconda and librosa EXACTLY as the instructions say
    • Hopefully you are up and running again
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
Presently working on implemeting:
  • Full-range plots
  • SRC to 96k to simplify the math and logic. Later I may add SRC to 192k for files >96k as we're plotting 2H and 3H
  • The ability to apply RIAA, inverse RIAA, and for both the ability only apply bass turnover and shelf
  • General restructuring in to functions for hygiene and to extend utility, e.g. plotting multiple axis

These next three plots are:
  1. Full-range plot of the (simulated) recording characteristic of JVC TRS-1007, which is RIAA bass turnover and shelf without treble emphasis. This applies to many test records. FILE
  2. Cartridge FR plot of a flat layback from a TRS-1007 sweep track. FILE
  3. The same layback as #2 but with RIAA bass turnover and shelf applied via the script.

TRS-1007 Recording Characteristic.png

20Hz-20kHz Test Flat.png

20Hz-20kHz Test 318 3180.png



*** The IIR filters I'm using were done by Scott. Background and method are published here: https://linearaudio.net/record-replay-riaa-correction-digital-domain.
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
There may be a better solution but it was the first thing that worked for me when I was trying to run it. Inelegant, sure, but at least it doesn't get in the way of the rest of the code. Anyways, it's our problem, not yours, lol.

Well, I would like it to be as easy as possible for folks so that more people use it; all in the pursuit of standardized measurements. I can make it detect the OS, etc. but when I get time to work on it, I want to work on the fun stuff.
 

al2002

Active Member
Forum Donor
Joined
May 18, 2016
Messages
263
Likes
200
This thread is intended to consolidate the technical discussion and further development of Scott Wurcer's frequency response measurement script to a central area. I'd like to keep this discussion about the tool rather than specific cartridges.

I'll edit this and the subsequent reserved posts in the coming weeks with relevant info, background, links, etc.
Where is the original discussion of this topic? Can you provide a link? Is crosstalk measured too?

Inter-channel crosstalk data is essential in evaluating phono cartridges. Minor resonances clearly show up in crosstalk.
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
Where is the original discussion of this topic? Ca nyou provide a link? Is crosstalk measured too?

Inter-channel crosstalk data is essential in evaluating phono cartridges. Minor resonances clearly show up in crosstalk.

I'll add the background info and links to the original discussion on diyAudio when I can. WYSIWYG.
 

Thomas_A

Major Contributor
Forum Donor
Joined
Jun 20, 2019
Messages
3,358
Likes
2,367
Location
Sweden
Great thread. Who can aid the step by step procedure for Mac OS? I have tried to learn more of Python and Terminal stuff, but I realise I need a good book for this.
 
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
Great thread. Who can aid the step by step procedure for Mac OS? I have tried to learn more of Python and Terminal stuff, but I realise I need a good book for this.

I can after I finish up this version. Maybe someone else can jump in interim.
 

USER

Addicted to Fun and Learning
Forum Donor
Joined
Mar 30, 2019
Messages
865
Likes
1,357
This space is reserved for a visually-aided overview of how to use the CBS STR-100 test record with Audacity in order to create a file ready for the script.

Most of the test records being used on this site (the "eBay" ones) are issue 3 (new), which is the only version that is recommended. The record has a highly visible sweep oscillator switching at around 5kHz. In addition, it is not a good record to use for crosstalk measurements. Distortion measurements are limited compared to the very best test records, but they go down to about -40dB for the 2nd harmonic, which is more than good enough to tell you about the general quality of the cartridge. Some batches of the CBS STR-100 test records have pressing issues that affect the results above 10kHz. This is compounded by the fact that a lot of copies are slightly warped as they have been stored in shrink wrap for decades. Record flatteners can certainly help and are recommended, but expect issues to remain. Despite all of this, the test record can provide you with excellent and useful measurements and is one of the very best still available today.

Worst case scenario:

Here is a comparison of the first batch I bought from eBay using a Shure V15 V-MR. #1 has the nastiest issue. However, it is still very much usable because the issue is obvious. Nonetheless, even if you have an off copy, the differences are generally within 0.5dB. Set-up alone can have a similar effect.
CBS STR-100 COMPARISONgif.gif

This all being said, this new batch look good, so YMMV. (These are not to be compared to those above.)
CBS STR-100 GOOD COPIES.gif

The CBS STR-100 issue 3 (new) (the "eBay" version) test record compares favorably with the most accurate--and adjusted--test records, with the biggest difference being a 0.5 dB bell dip from around 5-10kHz. (And the addition of whatever quirk your record may have above 10kHz, of course.) Expect a corrective EQ curve in the future.

You will need to record the first track of the CBS STR-100, which is comprised of a left channel sweep and a right channel sweep. The important thing to do is to remove the initial calibration tone before each respective sweep and cut at the end of the sweep. Cut right where the sweep begins and ends, not necessarily where the calibration tone ends. There should be no silence, only the sweep. Normalizing the file can help you find the spots. Then separate and save the left and right channel sweeps as separate .wav files to load onto the script, making sure you swap the channels for the right channel sweep so that it is the main track on top.

It is strongly advised that you keep a copy of the full, stereo, two-channel recordings as a backup.

There is a setting in the script (riaainv) that enables the processing your file with an inverse-RIAA filter in the likely case that your phono preamp uses RIAA. However, in case they are useful to you, RIAA and Inverse-RIAA nyquist filters for 96k recordings on Audacity are attached here. They are based on the biquads posted above. These filters are to replace ALL previously distributed ones, which have minor errors at best and are way off at worst. These are the most precise filters available on the web.
 

Attachments

  • RIAA & IRIAA Filters.zip
    590 bytes · Views: 109
Last edited:

USER

Addicted to Fun and Learning
Forum Donor
Joined
Mar 30, 2019
Messages
865
Likes
1,357
What is the hardware setup for this? Any old turntable/alignment of the cartridge will do???
Yes, any turntable and cartridge will do. Detailed information including loading settings and cable capacitance is necessary to ensure useful results. Loads and loads of them in this thread, starting around here:


Note that these will need to be re-adjusted with the changes JP is making in order to optimize the results, including the new IRIAA filter. There will be a new thread dedicated to cartridge measurements in the near future. It's going to be a terrific collection. Let's hope ASR can become a good, hospitable place for turntable and cartridge discussion.

The test record makes or breaks this as it usually does. JP uses a nearly impossible to find JVC TRS-1007. We plebs use the CBS STR-100 found on eBay for $17. (If this thread takes off expect prices to sky rocket once they sell out). A project is slowly underway to compare the JVC and CBS records properly, but they are both pretty good even with their quirks.

Ideally you want to record at 96k to get the full 2nd harmonic. 192k or so for 3rd. You may have to run a IRIAA filter depending on the record and/or phono preamp.
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
Is it this thread?


No - the development of it started around here: https://www.diyaudio.com/community/thre ... st-5749320
 
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
In regard to the CBS STR 100 test record, it is indeed what most people use as they are still easy to get, particularly from the eBay seller 92hdcma. This record worked reasonably well when we were plotting from 1kHz, but full-range plots will present some interesting hurdles:

CBS STR 100 Recording Characteristic.png


"The resultant curve for a typical high-fidelity pickup is shown in Fig. 1. Note that in the STR 100 record the 500-Hz break is sharp; the modulation below 500-Hz is precisely at constant amplitude; above 500-Hz at constant velocity. Thus, ideally response can be verified if the trace follows two straight lines - one rising at 6dB per octave up to 500Hz and the other horizontal above 500Hz."

Doc - Jan 18 2023 - 8-44 PM[11].jpeg


STR 100 Test Flat.png

^ FILE

There's a lot of LF modulation in this record, such that if I try to correct it the LF "grunge" is the largest bin in each time slice which renders the method useless. Rather than pre-filtering and figuring out a transfer function (which I've never done before), it may be easier, at least for now, to just apply an amplitude correction to each bin between 40Hz and 500Hz. That math is easy.

Below illustrates the problem. Top two are flat, and the bottom two I applied RIAA base turnover and shelf.

Screen Shot 2023-01-19 at 12.31.16 PM.png
 
Last edited:
OP
JP

JP

Major Contributor
Joined
Jul 4, 2018
Messages
2,246
Likes
2,421
Location
Brookfield, CT
What is the hardware setup for this? Any old turntable/alignment of the cartridge will do???

As @USER said the hardware really isn't critical as long as it's not broken. For full range measurements folks will want to pay attention to arm/cartridge matching as that can definitely affect performance <30Hz or so. I should work on adding an LF FFT inset so people can see where that resonance is and eyeball the Q. This doesn't matter with the currently released version as it plots from 1kHz. The most difficult part for people is being able to determine what the actual load their setup presents to an MM and also MI types, as that can affect HF performance to a large degree.

I'll write more later on how it works if I can rope that @SIY guy in to being my technical editor so I don't wind up saying anything really stupid. In a nutshell it's time slicing the signal and returning the largest bin for a given slice, which are then interpolated/averaged. Slices are Fs/100 in the current version but small slices don't work for LF and large slices don't work for HF. This next version uses ranges of slices that gives us 5Hz bins (smallest that'll work) for 0-50Hz, 10Hz for 50-100, 20Hz for 100-1kHz, and 100Hz for 1kHz+.

I've even used it as a bit of a verification tool by playing spot frequencies and doing a controlled ramp down of the platter speed. "Controlled" meaning I needed the end speed to be a certain percentage of normal speed. This is a 20kHz spot ramped down 50%, so the "sweep" is 20kHz to 10kHz. A perfectly flat cartridge would adhere to a 6dB/oct line.

150MLX_150pF 47k_STR-170 Fake Sweep.png
 
Top Bottom