JP
Major Contributor
This is pretty granular.
As a quick hack you can set ylim. For example, put this
Code:
plt.ylim(-15,15)
above the line
Code:
plt.semilogx(freqout,ampout,color = 'b')
This is pretty granular.
plt.ylim(-15,15)
plt.semilogx(freqout,ampout,color = 'b')
I don't know how to help answer your question as I've never seen a UP-4 in person and know almost nothing about it.
Maybe write to Nagaoka?
I have a Shure V15VxMR that I really like and it still sounds good but I know the stylus is approaching the end of its life. I've considered the Jico SAS replacement but others who've used it in the V15 indicate it raises the high frequencies and makes it sound very bright. I believe its because the boron (and now sapphire) cantilever has a resonance in the audio band whereas the original in the Shure is beryllium. Would changing the loading and capacitance reduce this effect? On the other hand the Jico stylus price is pretty steep and I wonder if sonically the Nagaoka MP500 is the best replacement for the Shure V15VxMR?
Is it preferable to separate these in to two plots so that they can scale individually?
View attachment 111569
The regular ripples on the top plot could be at the rotational frequency (LP eccentricity). This is an interesting artifact, you could verify this by taking the seconds per decade of the sweep and dividing by 1.8, I don't have a TRS1007.Is it preferable to separate these in to two plots so that they can scale individually?
View attachment 111569
The regular ripples on the top plot could be at the rotational frequency (LP eccentricity). This is an interesting artifact, you could verify this by taking the seconds per decade of the sweep and dividing by 1.8, I don't have a TRS1007.
BTW thanks for helping out here, I've lost my furnace and am distracted trying to keep any pipes from freezing.
The regular ripples on the top plot could be at the rotational frequency (LP eccentricity). This is an interesting artifact, you could verify this by taking the seconds per decade of the sweep and dividing by 1.8, I don't have a TRS1007.
BTW thanks for helping out here, I've lost my furnace and am distracted trying to keep any pipes from freezing.
@scott wurcer @JPJ
Is there a way to set the x-axis to something less than 100kHz?
I tried playing with:
x = minf*100
but it didn't seem to do anything.
"""
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.
Python does not read 24bit packed .wav's 16 bit is more than enough here.
The infoLine 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
"""
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
#edit here to add a HOME directory, etc.
#HOME = '/Users/jjones/Documents/polar/' #THIS IS NOT NEEDED IF RAN FROM A FILE
_FILE = 'V15VMR350pF_TRS1007.wav'
infoLine = 'V15-VMR / 47K / 350pF / TRS-1007'
plotStyle = 3
#end Edit
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
try:
input = raw_input
except NameError:
pass
#_FILE = input('Enter a .wav file to process ')
#File needs trimming with no leading or trailing silence (undefined f)
y = read(_FILE)
Fs = float(y[0])
input_sig = y[1]
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
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 min(ampout) >= -5:
ax1.set_ylim(-10,)
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")
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()
Here's a bit of a 'pick your poison' version. Ex:
View attachment 112043
View attachment 112044
View attachment 112045
Code:""" 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. Python does not read 24bit packed .wav's 16 bit is more than enough here. 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 """ 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 #edit here to add a HOME directory, etc. #HOME = '/Users/jjones/Documents/polar/' #THIS IS NOT NEEDED IF RAN FROM A FILE _FILE = 'V15VMR350pF_TRS1007.wav' infoLine = 'V15-VMR / 47K / 350pF / TRS-1007' plotStyle = 3 #end Edit 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 try: input = raw_input except NameError: pass #_FILE = input('Enter a .wav file to process ') #File needs trimming with no leading or trailing silence (undefined f) y = read(_FILE) Fs = float(y[0]) input_sig = y[1] 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 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.grid(True, which="major", axis="both", ls="-", color="black") ax1.grid(True, which="minor", axis="both", ls="-", color="gainsboro") 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() ax1.grid(True, which="major", axis="both", ls="-", color="black") ax1.grid(True, which="minor", axis="both", ls="-", color="gainsboro") if min(ampout) >= -5: ax1.set_ylim(-10,) 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)) ax1.grid(True, which="major", axis="both", ls="-", color="black") ax1.grid(True, which="minor", axis="both", ls="-", color="gainsboro") 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.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()
Where can I get the scripts and instructions? I looked back at the thread but only found bits and pieces and no link for the scripts... probably missed it.