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