1#General imports 2import subprocess, pylab, os, tempfile, shutil 3import numpy as np 4import scipy.signal 5import matplotlib.pylab as plt 6import itertools 7from itertools import count 8 9#Filters design imports 10from scipy import signal 11import math, numpy 12from matplotlib import pyplot 13 14#Curve fitting 15from scipy.optimize import curve_fit 16 17#Filters wrappers 18from scipy import signal 19import matplotlib.pyplot as plt 20 21#Text parsers 22import re 23from parse import * 24 25""" 26 OC-2 Octaver precomputation script. 27 28 The next basic steps have place: 29 - Convert .sch files to netlists. 30 - Run ngspice and extract required data from here. 31 - Parse data. 32 - Extract required data from ngspice. 33 - Write it to file, which will be included to faust .dsp script 34 35 Partially crafted from Guitarix/trunk/tools/ampsim/DK/spice.py 36 May be will be merged with ampsim in future. 37""" 38 39# Spice get characteristics 40def _add_ac_cmd(cmd, output, step, fmin, fmax, prefix = "ngspice"): 41 l = list(cmd) 42 l.append(".ac dec %g %g %g" % (step, fmin, fmax)) 43 l.append(".control") 44 l += ["run", 45 "wrdata %s %s" % (prefix, ",".join(["db(v(%s))" % output])), 46 ".endc", 47 ".end", 48 ] 49 return "\n".join(l) 50 51def _add_dc_cmd(cmd, output, source, step, fmin, fmax, prefix = "ngspice"): 52 l = list(cmd) 53 l.append(".dc %s %g %g %g" % (source, fmin, fmax, step)) 54 l.append(".control") 55 l += ["run", 56 "wrdata %s %s" % (prefix, ",".join(["v(%s)" % output])), 57 ".endc", 58 ".end", 59 ] 60 return "\n".join(l) 61 62# Read netlist from gschem .sch file 63def _mk_spice_netlist(fname): 64 try: 65 p = subprocess.Popen( 66 ("gnetlist","-g","spice-sdb","-o","/dev/fd/1","-q",fname), 67 stdout=subprocess.PIPE, stderr=subprocess.PIPE) 68 except OSError as e: 69 raise RuntimeError( 70 "can't start gnetlist -- please check gnetlist installation [%s]" % e) 71 out, err = p.communicate() 72 if p.returncode: 73 raise ValueError("error calling gnetlist: %s" % err) 74 if err: 75 print err 76 l = [] 77 for line in out.split("\n"): 78 line = line.rstrip() 79 if line == ".end": 80 break 81 l.append(line) 82 83 return l 84 85def _write_input_file(iname, lines): 86 with open(iname, "w") as f: 87 for i in lines: 88 f.write(i) 89 f.write("\n") 90 return len(lines) 91 92def _stream(script, out_file, prefix = "ngspice"): 93 try: 94 p = subprocess.Popen( 95 ("ngspice","-n","-b","/dev/fd/0"), cwd=out_file, stdin=subprocess.PIPE, 96 stdout=subprocess.PIPE, stderr=subprocess.PIPE) 97 except OSError as e: 98 raise RuntimeError( 99 "can't start ngspice -- please check ngspice installation [%s]" % e) 100 out, err = p.communicate(script) 101 # bug in ngspice: returncode always 1 102 def checkerr(buf): 103 #Find "error" key word in output stream 104 print buf 105 if "Error" in buf: 106 l = [line for line in buf.split("\n") if line.startswith("Error")] 107 raise RuntimeError("error in ngspice call:\n%s" % "\n".join(l)) 108 #checkerr(out) 109 #checkerr(err) 110 return 111 112def _read_result(rname): 113 with open(rname) as f: 114 lines = f.readlines() 115 result = None 116 for i, line in enumerate(lines): 117 parts = [float(part) for part in line.split()] 118 if result is None: 119 result = np.zeros((len(lines), len(parts))) 120 result[i] = parts 121 return result 122 123# Returns array with frequency and magnitude in db 124def get_ac(sch_fname, out_name, step, fmin, fmax): 125 outdir = tempfile.mkdtemp() 126 _stream(_add_ac_cmd(_mk_spice_netlist(sch_fname), out_name, step, fmin, fmax), outdir) 127 return _read_result(outdir+"/ngspice.data") 128 129# Returns array with output DC values 130def get_dc(sch_fname, out_name, source, step, vmin, vmax): 131 outdir = tempfile.mkdtemp() 132 _stream(_add_dc_cmd(_mk_spice_netlist(sch_fname), out_name, source, step, vmin, vmax), outdir) 133 return _read_result(outdir+"/ngspice.data") 134 135 136# Results processing 137def normalize(arr): 138 #Extract magnituce 139 mag = [] 140 for i, val in enumerate(arr): 141 mag.append(arr[i][1]) 142 #Normalize it 143 mag=mag - max(mag) 144 #Apply normalized 145 for i, val in enumerate(arr): 146 arr[i][1] = mag[i] 147 148 return arr 149 150def _find_first_more(arr, value): 151 for i, val in enumerate(arr): 152 if value < arr[i]: 153 return i 154 155def _find_first_less(arr, value): 156 for i, val in enumerate(arr): 157 if value > arr[i]: 158 return i 159 160def hp_freq(arr, val_db): 161 mag = [] 162 for i, val in enumerate(arr): 163 mag.append(arr[i][1]) 164 return arr[_find_first_more(mag, val_db)][0] 165 166def lp_freq(arr, val_db): 167 mag = [] 168 for i, val in enumerate(arr): 169 mag.append(arr[i][1]) 170 return arr[_find_first_less(mag, val_db)][0] 171 172# Custom plot functions 173def parse_data(data): 174 x = [] 175 y = [] 176 for i, val in enumerate(data): 177 x.append(data[i][0]) 178 y.append(data[i][1]) 179 180 return (x, y) 181 182def plot_ac_to_compare(filter_data_from_spice, coef_b, coef_a, FS, visualise=False): 183 if visualise: 184 #get and plot data from spice 185 w1, h1 = parse_data(filter_data_from_spice) 186 #get and plot data from coesfs 187 w, h = signal.freqz(coef_b, coef_a) 188 plt.figure() 189 plt.plot(w1, h1) 190 plt.semilogx(w/max(w)*FS/2, 20 * np.log10(abs(h))) 191 plt.xlabel('Frequency [ Hz ]') 192 plt.ylabel('Amplitude [dB]') 193 plt.margins(0, 0.1) 194 plt.grid(which='both', axis='both') 195 plt.show() 196 197 198def plot_dc_to_compare(data_to_fit, polyfunc, visualise=False): 199 if visualise: 200 vin, vout = parse_data(data_to_fit) 201 x_new = np.linspace(vin[0], vin[-1], 50) 202 y_new = polyfunc(x_new) 203 plt.figure() 204 plt.plot(vin,vout,'o', x_new, y_new) 205 plt.show() 206 207# Functions, which are needed for data exporting 208def save_filter_to_file(fileName, b, a, prefixStr, description=""): 209 try: 210 with open(fileName, "a+") as f: 211 l = list("\n") 212 l.append("//Description: %s" % (description)) 213 i = 0 214 for b_val, a_val in itertools.izip_longest(b, a, fillvalue=0): 215 l.append("b%g_%s = %s; a%g_%s = %s;" % (i, prefixStr, b_val, i, prefixStr, a_val)) 216 i += 1 217 f.write("\n".join(l)) 218 f.close() 219 except IOError as e: 220 print "I/O error({0}): {1}".format(e.errno, e.strerror) 221 222def save_dc_curve_to_file(fileName, poly, prefixStr, description=""): 223 coefs = poly.c 224 num_of_coefs = len(coefs) 225 try: 226 with open(fileName, "a+") as f: 227 l = list("\n") 228 l.append("//Description: %s" % (description)) 229 l.append("%s = " % (prefixStr)) 230 for i, coef in enumerate(coefs): 231 l.append("%g * pow(x, %d) + " % (coef, num_of_coefs -1 - i)) 232 if i == num_of_coefs - 2: 233 break 234 235 l.append("%g;" % (coefs[num_of_coefs -1])) 236 f.write("\n".join(l)) 237 f.close() 238 except IOError as e: 239 print "I/O error({0}): {1}".format(e.errno, e.strerror) 240 241# Read components values from netlist 242def _getValByRefDes(lines, refDes): 243 #Use reversed() to read main schematic from 244 #.end and avoid included models 245 for line in reversed(lines): 246 if line == "*============== Begin SPICE netlist of main design ============": 247 break 248 249 split = parse("%s {} {} {}" % (refDes), line) 250 if split != None: 251 return split[2] 252 253def _metricConvert(val): 254 num = re.findall(r'\d+', val) 255 num = float(''.join(map(str,num))) 256 letter = re.findall("[a-zA-Z]+", val) 257 return { 258 'M': num*1e6, 259 'k': num*1e3, 260 'm': num*1e-3, 261 'u': num*1e-6, 262 'n': num*1e-9, 263 'p': num*1e-12, 264 }[letter[0]] 265 266def getValByRefDes(lines, refDes): 267 #Find value, parse and apply metric conversion 268 return _metricConvert(_getValByRefDes(lines, refDes)) 269 270# Some info 271descriptionHeader = """ 272//BOSS OC-2 partial emulation header file. 273//Generated automatically using oc_2.py script. 274//So, do not modify it. 275//Also, see description in the .odg file 276""" 277 278# Set to True, if you want to learn figures with plots 279showPlots = True 280 281# Predefinitions 282FS=48000.0 283filtersOrder = 3 284polyOrder = 12 285 286# Filters emulation 287f1_mag = get_ac("oc_2_partial_emulation_0.sch", "out_f1", 20, 20, FS/2) 288f2_mag = get_ac("oc_2_partial_emulation_0.sch", "out_f2", 20, 20, FS/2) 289f3_mag = get_ac("oc_2_partial_emulation_0.sch", "out_f3", 20, 20, FS/2) 290 291# Emulate cached filters using butterworth function, 292# *_cutoff values picked up visually from plots 293f1_cutoff = lp_freq(f1_mag, -3) 294b1, a1 = signal.butter(filtersOrder, 2*f1_cutoff/FS, 'lowpass') 295plot_ac_to_compare(f1_mag, b1, a1, FS, showPlots) 296 297f2_cutoff = lp_freq(f2_mag, -3) 298b2, a2 = signal.butter(filtersOrder, 2*f2_cutoff/FS, 'lowpass') 299plot_ac_to_compare(f2_mag, b2, a2, FS, showPlots) 300 301f3_cutoff = lp_freq(f3_mag, -3) 302b3, a3 = signal.butter(filtersOrder, 2*f3_cutoff/FS, 'lowpass') 303plot_ac_to_compare(f3_mag, b3, a3, FS, showPlots) 304 305# Emulate trigger filters using bilinear transform 306#The values are cached from gschem file "oc_2_partial_emulation_1.sch" 307netlist = _mk_spice_netlist("oc_2_partial_emulation_1.sch") 308 309R1=getValByRefDes(netlist, "R1") 310C1=getValByRefDes(netlist, "C1") 311R2=getValByRefDes(netlist, "R2") 312C2=getValByRefDes(netlist, "C2") 313R9=getValByRefDes(netlist, "R9") 314 315# Some bilinear transform usage: see oc_2.dsp and *_1.sch to understand, what is going on. 316b41_1, a41_1 = scipy.signal.bilinear ([1.0], [R1*C1, 1.0], FS) 317b41_2, a41_2 = scipy.signal.bilinear ([C2*R9,1.0], [C2*(R2+R9),1.0], FS) 318b42, a42 = scipy.signal.bilinear ([C2*(R2+R9), 1.0], [C1*C2*(R2+R9)*R1, C2*(R2+R9)+R1*(C1+C2), 1.0], FS) 319 320R6=getValByRefDes(netlist, "R6") 321C4=getValByRefDes(netlist, "C4") 322b_int, a_int = scipy.signal.bilinear ([1.0], [R6*C4, 1.0], FS) 323 324# Switches emulation 325opened_state = get_dc("oc_2_partial_emulation_0.sch", "sw_out_opened", "V3", 0.1, -5, 5) 326vin1, vout1 = parse_data(opened_state) 327f_opened = np.poly1d(np.polyfit(vin1, vout1, polyOrder)) 328plot_dc_to_compare(opened_state, f_opened, showPlots) 329 330closed_state = get_dc("oc_2_partial_emulation_0.sch", "sw_out_closed", "V4", 0.1, -5, 5) 331vin2, vout2 = parse_data(closed_state) 332f_closed = np.poly1d(np.polyfit(vin2, vout2, polyOrder)) 333plot_dc_to_compare(closed_state, f_closed, showPlots) 334 335# All data ready, start file exporting 336outputFilePath = os.path.dirname(os.path.realpath(__file__)) 337outputFilePath += "/oc_2.lib" 338 339f = open(outputFilePath, "w") 340f.write(descriptionHeader) 341f.close() 342 343# Save filters 344save_filter_to_file(outputFilePath, b1, a1, "f1", "Filter 1") 345save_filter_to_file(outputFilePath, b2, a2, "f2", "Filter 2") 346save_filter_to_file(outputFilePath, b3, a3, "f3", "Filter 3") 347 348save_filter_to_file(outputFilePath, b41_1, a41_1, "f41_1", "Filter 4") 349save_filter_to_file(outputFilePath, b41_2, a41_2, "f41_2", "Filter 4") 350save_filter_to_file(outputFilePath, b42, a42, "f42", "Filter 4") 351save_filter_to_file(outputFilePath, b_int, a_int, "fint", "Integration filter") 352 353# Save switches curves 354save_dc_curve_to_file(outputFilePath, f_opened, 355 "sw_opened(x)", "Switch opened state curve") 356 357save_dc_curve_to_file(outputFilePath, f_closed, 358 "sw_closed(x)", "Switch closed state curve") 359 360 361