Balle Clorin
Major Contributor
- Joined
- Dec 26, 2017
- Messages
- 2,129
- Likes
- 1,955
/Users/ersol/Documents/data/HiFi/Michell/RPM WFWFtestwav.wav
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: #!!!!!!!!!!!!!!!!!!!! 10 = 1kHz. Change to 0 to start in the basement.
frequency.append(f*100)
amplitude.append(y[f])
"""
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.
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"
plotstyle =
1 - traditional
2 - dual axis (twinx)
3 - dual plot
VERSION 8C
"""
from scipy import signal
from scipy.io.wavfile import read
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import datetime
import os
import librosa
#edit here to add a HOME directory, etc.
#HOME = '/Users/user/Documents/' #THIS IS NOT NEEDED IF RAN FROM A FILE
_FILE = 'wavefile.wav'
infoline = 'Cart / Load / Record'
roundlvl = 2
plotstyle = 2
plotdataout = 0
ovdylim = 0
ovdylimvalue = [-25,5]
topdb = 20
framelength = 512
hoplength = 128
#end edit
print('Input File: ' + str(_FILE))
print('Info Line: ' + str(infoline))
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
ax1.set_ylim(new_lim1)
ax2.set_ylim(new_lim2)
#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
audio, Fs = librosa.load(_FILE, sr=None)
filelength = audio.shape[0] / Fs
print('Sample Rate: ' + str("{:,}".format(Fs) + 'Hz'))
print(f"Length: {filelength}s")
trimmedaudio, index = librosa.effects.trim(audio, top_db=topdb, frame_length=framelength, hop_length=hoplength)
input_sig = (trimmedaudio* 32767).astype(int) #scale samples and convert to int
print(f"In/Out (s): {index / Fs}")
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
deltah = round((max(ampout)), roundlvl)
deltal = abs(round((min(ampout)), roundlvl))
print('Min Freq: ' + str("{:,}".format(minf * 100) + 'Hz'))
print('Max Freq: ' + str("{:,}".format(maxf * 100) + 'Hz'))
print('Min Ampl: ' + str(round((min(ampout)), 4)) + 'dB')
print('Max Ampl: ' + str(round((max(ampout)), 4)) + 'dB')
print('Delta: ' + str(round((max(ampout) - min(ampout)), 4)) + 'dB')
print('Frequency: ' + str(len(freqout)) + ' list elements')
print('Amplitude: ' + str(len(ampout)) + ' list elements')
print('Amplitude 2h: ' + str(len(ampout2h)) + ' list elements')
print('Amplitude 3h: ' + str(len(ampout3h)) + ' list elements')
if plotdataout == 1:
dampout = [*ampout, *[''] * (len(freqout) - len(ampout))]
dampout2h = [*ampout2h, *[''] * (len(freqout) - len(ampout2h))]
dampout3h = [*ampout3h, *[''] * (len(freqout) - len(ampout3h))]
print('\nPlot Data: (freq, ampl, 2h, 3h)\n')
dataout = list(zip(freqout, dampout, dampout2h, dampout3h))
for fo, ao, ao2, ao3 in dataout:
print(fo, ao, ao2, ao3, sep=', ')
plt.rcParams["xtick.minor.visible"] = True
plt.rcParams["ytick.minor.visible"] = True
if plotstyle == 1:
fig, ax1 = plt.subplots(1, 1, figsize=(14,6))
ax1.semilogx(freqout,ampout,color = 'b', label = 'Freq Response')
ax1.semilogx(freqout2h,ampout2h,color = 'g', label = '2nd Harmonic')
ax1.semilogx(freqout3h,ampout3h,color = 'r', label = '3rd Harmonic')
ax1.set_ylabel("Amplitude (dB)")
ax1.set_xlabel("Frequency (Hz)")
ax1.legend(loc=4)
if plotstyle == 2:
fig, ax1 = plt.subplots(1, 1, figsize=(14,6))
ax2 = ax1.twinx()
if max(ampout) <7:
ax1.set_ylim(-25, 7)
if max(ampout) < 4:
ax1.set_ylim(-25,5)
if max(ampout) < 2:
ax1.set_ylim(-29,3)
if max(ampout) < 0.5:
ax1.set_ylim(-30,2)
if ovdylim == 1:
ax1.set_ylim(*ovdylimvalue)
ax1.semilogx(freqout,ampout,color = 'b', label = 'Freq Response')
ax2.semilogx(freqout2h,ampout2h,color = 'g', label = '2nd Harmonic')
ax2.semilogx(freqout3h,ampout3h,color = 'r', label = '3rd Harmonic')
align_yaxis(ax1, ax2)
ax1.set_ylabel("Amplitude (dB)")
ax2.set_ylabel("Distortion (dB)")
ax1.set_xlabel("Frequency (Hz)")
fig.legend(loc=4, bbox_to_anchor=(0.899, 0.11))
if plotstyle == 3:
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(14,6))
ax2.grid(True, which="major", axis="both", ls="-", color="black")
ax2.grid(True, which="minor", axis="both", ls="-", color="gainsboro")
if max(ampout) <= 1.75 and min(ampout) >= -1.75:
ax1.set_ylim(-2,2)
ax1.semilogx(freqout,ampout,color = 'b', label = 'Freq Response')
ax2.semilogx(freqout2h,ampout2h,color = 'g', label = '2nd Harmonic')
ax2.semilogx(freqout3h,ampout3h,color = 'r', label = '3rd Harmonic')
ax1.set_ylabel("Amplitude (dB)")
ax2.set_ylabel("Distortion (dB)")
ax2.set_xlabel("Frequency (Hz)")
ax1.legend(loc=4)
ax2.legend(loc=4)
ax1.grid(True, which="major", axis="both", ls="-", color="black")
ax1.grid(True, which="minor", axis="both", ls="-", color="gainsboro")
bbox_args = dict(boxstyle="round", color='b', fc='w', ec='b', alpha=1)
arrow_args = dict(arrowstyle="->")
ax1.annotate('+' + str(deltah) + ', ' + u"\u2212" + str(deltal) + ' dB', color='b', \
xy=(1000,0), xycoords='data', \
xytext=(-15, -20), textcoords='offset points', \
ha="left", va="center", \
bbox=bbox_args, \
#arrowprops=arrow_args \
)
ax1.set_xticks([1000,10000])
ax1.set_xticklabels(['1k','10k'])
ax1.set_title(infoline + "\n", fontsize=16)
mod_date = datetime.datetime.fromtimestamp(os.path.getmtime(_FILE))
plt.figtext(.17, .13, infoline + "\n" + _FILE + "\n" + \
mod_date.strftime("%b %d, %Y %H:%M:%S"), fontsize=8)
plt.savefig(infoline.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.75, dpi=96)
plt.show()
print('\nDone!')
#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.
HOME = '/Users/jjones/Documents/polar/'
_FILE = 'FM_WF_SONY.wav'
info_line = 'WF_SONY'
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
stereo_channel = int(stereo_channel)
filter_freq = float(filter_freq)
sec_per_rev = float(sec_per_rev)
Hz_per_tick = float(Hz_per_tick)
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(HOME + _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)]
mm = np.mean(if1)
mm2 = np.mean(if2)
for x in range(0,len(if1)):
if if1[x] > mm+10:
if1[x] = mm+10
if if1[x] < mm-10:
if1[x] = mm-10
for x in range(0,len(if2)):
if if2[x] > mm2+10:
if2[x] = mm2+10
if if2[x] < mm2-10:
if2[x] = mm2-10
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.
#ax.text(222.*dgr, 27, 'Mean Rev1 {:4.2f}Hz'.format(np.mean(if1)), fontsize=8)
#ax.text(223.*dgr, 27.5, 'Mean Rev2 {:4.2f}Hz'.format(np.mean(if2)), fontsize=8)
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.savefig(info_line.replace(' / ', '_') +'.png', bbox_inches='tight', pad_inches=.5, transparent=True)
plt.show()
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: # <-------- Change this to 5. f * 100 = Hz.
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])
I think I mentioned in a private conversation that changing that number doesn't work for me. I get the following:"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."
We can go full old-school chart-recorder and put guides on the plot for perfect response. I'd hope I can get a filter to correct it but not sure - I'm brute-forcing my way through most of this.
I haven't paid much attention to the difference in records - trying to get the scenarios debugged so I can consolidate the code. For instance, when I apply lower turnover to the STR-100 it breaks the interpolation routine. No idea why.
BTW, if you change the marked value to 5 you can plot 500Hz with the existing code. The downside is that the resolution is 100Hz.
View attachment 258194Code: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: # <-------- Change this to 5. f * 100 = Hz. 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])
norm = ampout[0]
ampout = ampout-norm #amplitude is in dB so normalize by subtraction at [0]
ampout2h = ampout2h-norm
ampout3h = ampout3h-norm
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
i = find_nearest(freqout, 1000)
norm = ampout[i]
ampout = ampout-norm #amplitude is in dB so normalize by subtraction at [i]
ampout2h = ampout2h-norm
ampout3h = ampout3h-norm
ax1.annotate('+' + str(deltah) + ', ' + u"\u2212" + str(deltal) + ' dB', color='b', \
xy=(1000,0), xycoords='data', \
xytext=(-15, -20), textcoords='offset points', \
ha="left", va="center", \
bbox=bbox_args, \
#arrowprops=arrow_args \
)
ax1.annotate('+' + str(deltah) + ', ' + u"\u2212" + str(deltal) + ' dB', color='b', \
xy=(freqout[0],(ampout[0]-1)), xycoords='data', \
xytext=(-15, -20), textcoords='offset points', \
ha="left", va="center", \
bbox=bbox_args, \
#arrowprops=arrow_args \
)
Would it be too much to ask for the same measurement without the KAB fluid damper? I’m considering one for myself but honestly not sure if it will provide any sonic benefit. I’m using a Shure V15V-MR so I’m interested to see the effect it has on your V15iii.Shure V15 III with Jico SAS/B stylus, Technics SL-1210GR with KAB fluid damper and stock headshell.
Latest script used, capacitance of arm and IC was measured with proper tool so it should be very close to real (+/- 10pF).
View attachment 264261