Balle Clorin
Major Contributor
- Joined
- Dec 26, 2017
- Messages
- 1,383
- Likes
- 1,256
Trying to run your frequency response plotting script, I something wrong with my file directory, file location here?Here's an updated script. Also includes some rude scaling logic and config to override and set a custom ylim. I haven't tested this version on Windows so there'll likely be a couple compatibility issues with the file attribute handling.
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 VERSION 6BClean """ 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 = 'v15v-str100.wav' infoline = 'V15VMR / 47k 350pF / STR-100' roundlvl = 2 plotstyle = 2 plotdataout = 0 ovdylim = 0 ovdylimvalue = [-25,5] #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 try: input = raw_input except NameError: pass y = read(_FILE) Fs = int(y[0]) input_sig = y[1] filelength = input_sig.shape[0] / Fs print(f"Length: {filelength}s") print('Sample Rate: ' + str("{:,}".format(Fs) + 'Hz')) 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!')
Moved my Wav file to the directory where Python is installed and now I get this error?
The sweep file is 50 seconds
Cut off the last millisec part of the file with low signal and Voila!
'
'
No Ju just need to find out how to plot the whole frequency range and a script for the Wow&Flutter polar plot
Comparing plots by overlay
Attachments
Last edited:
Well, that is going to be the best most of us will be able to do. Still, giving us an extra 500Hz is a generous gift.
