I tried a little experiment, implementing some code in numpy (usually I
build modules in c++ to interface to python). Since these operations are
all large vectors, I hoped it would be reasonably efficient.
The code in question is simple. It is a model of an amplifier, modeled by
it's AM/AM and AM/PM characteristics.
The function in question is the __call__ operator. The test program plots a
spectrum, calling this operator 1024 times each time with a vector of 4096.
Any ideas? The code is not too big, so I'll try to attach it.
from linear_interp import linear_interp
from numpy import *
from math import pi
from uvector import vector_Complex
def db_to_pwr (db):
return 10**(0.1*db)
def db_to_volt (db):
return 10**(0.05*db)
def db_to_gain (db_in, db_out):
return 10**(0.05*(db_out-db_in))
def from_gain_and_phase (gain, phase):
out = empty (len (gain), dtype=complex)
out.real = cos (phase) * gain
out.imag = sin (phase) * gain
return out
class ampl (object):
def __init__ (self, pin, pout, deg, max_pin, delta_v, ibo):
"pin, pout in dB normalized to sat, ibo in dB"
ampl_interp = linear_interp (vectorize (db_to_volt) (pin), db_to_volt (pout))
phase_interp = linear_interp (vectorize (db_to_volt)(pin), array(deg)*pi/180)
## These would be used if input was linear instead of db
## ampl_interp = linear_interp (sqrt (pin), sqrt (pout))
## phase_interp = linear_interp (sqrt (pin), array(deg)*pi/180)
eps = 1e-6
max_vin = sqrt (max_pin)
vin = arange (0, max_vin+eps, delta_v)
gain = vectorize (ampl_interp) (vin)/vin
##gain = vectorize (ampl_interp) (vin**2) / vin
#gain[0] = 1
gain[0] = gain[1] # avoid singularity
phase = vectorize (phase_interp) (vin**2)
self.cmplx_gain = from_gain_and_phase (gain, phase)
self.delta_v = delta_v
self.ibo = 10**(-0.05*ibo)
self._phase_comp = from_gain_and_phase (ones (len (phase)), -phase)
def __call__ (self, zin):
zin_ibo = zin * self.ibo
vin = abs (array (zin_ibo, dtype=complex))
index = vin / self.delta_v
lower = floor (index).astype (int)
upper = lower + 1
assert (alltrue (upper < len (self.cmplx_gain)))
delta = (index - lower)
cmplx_gain = self.cmplx_gain[lower]*(1-delta) + self.cmplx_gain[upper]*delta
return vector_Complex (cmplx_gain * zin_ibo)
def phase_comp (self, zin):
zin_ibo = zin * self.ibo
vin = abs (array (zin_ibo, dtype=complex))
index = vin / self.delta_v
lower = floor (index).astype (int)
upper = lower + 1
assert (alltrue (upper < len (self._phase_comp)))
delta = (index - lower)
return vector_Complex (self._phase_comp[lower]*(1-delta) + self._phase_comp[upper]*delta)
if __name__ == "__main__":
x = array ((-20, -10, 0, 10), dtype=float)
y = array ((-20, -10, 0, 0), dtype=float)
p = array ((-20, -10, 0, 0), dtype=float)
t = ampl (x, y, p, 10, .01, 0)
def pwr_to_cmplx (db):
return complex (10**(0.05*db))
db = (-20, -10, 0, 1)
res = array (t (vectorize (pwr_to_cmplx) (db)), dtype=complex)
#print res
def mag_sqr (z):
return z.real*z.real + z.imag*z.imag
db_out = log10 (vectorize(mag_sqr) (res)) * 10
phase_out = arctan2 (res.imag, res.real)*180/pi
print db_out, phase_out
class linear_interp (object):
def __init__ (self, x, y):
assert (len (x) == len (y))
self.n = len (x)
self.x = tuple (x)
self.y = tuple (y)
def __call__ (self, _x):
klo = 1
khi = self.n
while (khi - klo > 1):
k = (khi + klo) >> 1
if (self.x[k - 1] > _x):
khi = k
else:
klo = k
h = float (self.x[khi - 1]-self.x[klo - 1])
if (h == 0.0):
raise ValueError, "Bad x input to routine linear_interp"
a = (self.x[khi - 1]-_x) / h
b = (_x - self.x[klo - 1]) / h
return a * self.y[klo - 1] + b * self.y[khi - 1]
if (__name__ == "__main__"):
x = (xrange (10))
y = [e+1 for e in x]
l = linear_interp (x, y)
print l (0)
print l (1)
print l (0.5)
import sys
sys.path.append ('../mod')
from constellation import *
from boost_rand import rng, pn, normal_c, normal, uniform_real, uniform_int, gold_code_generator, beta
from nyquist import *
from fir import coef_from_func_double, FIR_Complex_double, InterpFIR_Complex_double, DecimFIR_Complex_double
from ampl import ampl
#from polar_ampl import ampl
from fft3 import *
from uvector import vector_double, vector_Complex, mag_sqr, log10
from stats import stat_Complex
from limit import Limit
import math
from optparse import OptionParser
parser = OptionParser()
parser.add_option ("--sps", type=int, default=4)
parser.add_option ("-n", "--pts", type=int, default=1024)
parser.add_option ("-m", "--sets", type=int, default=1024)
parser.add_option ("-a", "--alpha", type=float, default=0.25)
parser.add_option ("--si", type=int, default=16)
parser.add_option ("--ibo", type=float, default=0)
parser.add_option ("--linear", action="store_true")
(opt, args) = parser.parse_args()
scale = 1
xmit_coef = coef_from_func_double (RNyquistPulse (opt.alpha, scale), 1./opt.sps, opt.sps*opt.si)
mf_coef = coef_from_func_double (RNyquistPulse (opt.alpha, scale), 1./opt.sps, opt.sps*opt.si)/opt.sps
rcv_filt = DecimFIR_Complex_double (mf_coef, opt.sps, 0)
CONSTELLATION_POINTS = ( complex ( 0.9239, 0.3827 ), complex ( 0.9239, -0.3827 ),
complex ( -0.9239, 0.3827 ), complex ( -0.9239, -0.3827 ),
complex ( 0.3827, 0.9239 ), complex ( 0.3827, -0.9239 ),
complex ( -0.3827, 0.9239 ), complex ( -0.3827, -0.9239 ) )
const = QAMconstellation (CONSTELLATION_POINTS)
rng1 = rng()
pn3 = pn (rng1, 3)
the_fft = fft (opt.pts, FORWARD)
xmit_filt = InterpFIR_Complex_double (xmit_coef, opt.sps)
import numpy as np
Id= np.array ([50.0*10**(0.05*-100), 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0, 1050.0])
pout= np.array ([-100, 0.0, 7.4000000000000004, 11.4, 14.199999999999999, 16.300000000000001, 18.100000000000001, 19.300000000000001, 20.600000000000001, 21.600000000000001, 22.600000000000001, 23.399999999999999, 24.100000000000001, 24.899999999999999, 25.399999999999999, 26.100000000000001, 26.5, 26.899999999999999, 27.100000000000001, 27.300000000000001, 27.399999999999999, 27.399999999999999])
phase= np.array ([105.0, 105.0, 101.0, 97.0, 94.0, 91.0, 88.0, 85.0, 81.0, 79.0, 76.0, 74.0, 73.0, 70.0, 68.0, 63.0, 58.0, 52.0, 47.0, 39.0, 31.0, 23.0])
#phase= np.array ([105.0, 101.0, 97.0, 94.0, 91.0, 88.0, 85.0, 81.0, 79.0, 76.0, 74.0, 73.0, 70.0, 68.0, 63.0, 58.0, 52.0, 47.0, 39.0, 31.0, 23.0])
phase = np.zeros (len (pout))
#Normalize to saturation point
pout -= np.amax (pout)
Id /= Id[np.argmax (pout)]
phase -= phase[0]
pin = 20*np.log10(Id)
xaxis = np.array (xrange (the_fft.size),dtype=float)*opt.sps/the_fft.size
#print xaxis
for opt.ibo in 0,1,2,3:
the_sum = vector_double (opt.pts, 0)
pwr_stat = stat_Complex()
rcv_pwr_stat = stat_Complex()
#for c++ ampl
#amplifier = ampl (vector_double (pin), vector_double (pout), vector_double (phase), pin[-1]+10, 0.01, opt.ibo)
#for numpy python ampl
amplifier = ampl (pin, pout, phase, pin[-1]+10, 0.01, opt.ibo)
for block in xrange (opt.sets):
symbols = opt.pts / opt.sps
xsyms = pn3 (symbols)
xconst = const (xsyms)
mod_out = xmit_filt (xconst)
if (not opt.linear):
mod_out = amplifier (amplifier.phase_comp (mod_out) * mod_out)
the_sum += mag_sqr (the_fft (mod_out)/the_fft.size)
pwr_stat += mod_out
rcv_pwr_stat += rcv_filt (mod_out)
the_sum /= opt.sets
print >> sys.stderr, 'ibo:', opt.ibo, 'obo:', math.log10 (pwr_stat.Var())*10, 'rcv_pwr:', math.log10 (rcv_pwr_stat.Var()) * 10
from matplotlib import pylab
pylab.grid (True)
pylab.plot (xaxis, np.array (log10 (the_sum)*10))
the_sum = vector_double (opt.pts, 0)
pwr_stat = stat_Complex()
rcv_pwr_stat = stat_Complex()
for block in xrange (opt.sets):
symbols = opt.pts / opt.sps
xsyms = pn3 (symbols)
xconst = const (xsyms)
mod_out = xmit_filt (xconst)
if (not opt.linear):
mod_out = Limit (mod_out)
the_sum += mag_sqr (the_fft (mod_out)/the_fft.size)
pwr_stat += mod_out
rcv_pwr_stat += rcv_filt (mod_out)
the_sum /= opt.sets
print >> sys.stderr, 'obo:', math.log10 (pwr_stat.Var())*10, 'rcv_pwr:', math.log10 (rcv_pwr_stat.Var()) * 10
from matplotlib import pylab
pylab.grid (True)
pylab.plot (xaxis, np.array (log10 (the_sum)*10))
pylab.show()
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
