Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: m-labs/artiq
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: d34084be0f37
Choose a base ref
...
head repository: m-labs/artiq
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 3eef6229cca8
Choose a head ref
  • 4 commits
  • 3 files changed
  • 1 contributor

Commits on Dec 7, 2016

  1. Copy the full SHA
    7e0f3ed View commit details

Commits on Dec 8, 2016

  1. Copy the full SHA
    d303225 View commit details
  2. fir: add ParallelHBFCascade

    jordens committed Dec 8, 2016
    Copy the full SHA
    a629eb1 View commit details
  3. Copy the full SHA
    3eef622 View commit details
Showing with 270 additions and 10 deletions.
  1. +150 −0 artiq/gateware/dsp/fir.py
  2. +17 −10 artiq/gateware/dsp/sawg.py
  3. +103 −0 artiq/test/gateware/test_fir.py
150 changes: 150 additions & 0 deletions artiq/gateware/dsp/fir.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
from operator import add
from functools import reduce
import numpy as np
from migen import *


def halfgen4(width, n):
"""
http://recycle.lbl.gov/~ldoolitt/halfband
params:
* `up` is the passband/stopband width, as a fraction of
input sampling rate
* `n is the order of half-band filter to generate
returns:
* `a` is the full set of FIR coefficients, `4*n-1` long.
implement wisely.
"""

npt = n*40
wmax = 2*np.pi*width
wfit = (1 - np.linspace(0, 1, npt)[:, None]**2)*wmax

target = .5*np.ones_like(wfit)
basis = np.cos(wfit*np.arange(1, 2*n, 2))
l = np.linalg.pinv(basis)@target

weight = np.ones_like(wfit)
for i in range(40):
err = np.fabs(basis@l - .5)
weight[err > .99*np.max(err)] *= 1 + 1.5/(i + 11)
l = np.linalg.pinv(basis*weight)@(target*weight)
a = np.c_[l, np.zeros_like(l)].ravel()[:-1]
a = np.r_[a[::-1], 1, a]/2
return a


class FIR(Module):
"""Full-rate finite impulse response filter.
:param coefficients: integer taps.
:param width: bit width of input and output.
:param shift: scale factor (as power of two).
"""
def __init__(self, coefficients, width=16, shift=None):
self.width = width
self.i = Signal((width, True))
self.o = Signal((width, True))
n = len(coefficients)
self.latency = (n + 1)//2 + 1

###

# Delay line: increasing delay
x = [Signal((width, True)) for _ in range(n)]
self.sync += [xi.eq(xj) for xi, xj in zip(x, [self.i] + x)]

# Wire up output
o = []
for i, c in enumerate(coefficients):
# simplify for halfband and symmetric filters
if c == 0 or c in coefficients[i + 1:]:
continue
o.append(c*reduce(add, [
xj for xj, cj in zip(x[::-1], coefficients) if cj == c
]))

if shift is None:
shift = width - 1
self.sync += self.o.eq(reduce(add, o) >> shift)


class ParallelFIR(Module):
"""Full-rate parallelized finite impulse response filter.
:param coefficients: integer taps.
:param parallelism: number of samples per cycle.
:param width: bit width of input and output.
:param shift: scale factor (as power of two).
"""
def __init__(self, coefficients, parallelism, width=16, shift=None):
self.width = width
self.parallelism = p = parallelism
n = len(coefficients)
# input and output: old to young, decreasing delay
self.i = [Signal((width, True)) for i in range(p)]
self.o = [Signal((width, True)) for i in range(p)]
self.latency = (n + 1)//2//parallelism + 2 # minus .5

###

# Delay line: young to old, increasing delay
x = [Signal((width, True)) for _ in range(n + p - 1)]
self.sync += [xi.eq(xj) for xi, xj in zip(x, self.i[::-1] + x)]

if shift is None:
shift = width - 1

# wire up each output
for j in range(p):
o = []
for i, c in enumerate(coefficients):
# simplify for halfband and symmetric filters
if c == 0 or c in coefficients[i + 1:]:
continue
o.append(c*reduce(add, [
xj for xj, cj in zip(x[-1 - j::-1], coefficients) if cj == c
]))
self.sync += self.o[j].eq(reduce(add, o) >> shift)


def halfgen4_cascade(rate, width, order=None):
"""Generate coefficients for cascaded half-band filters.
:param rate: upsampling rate. power of two
:param width: passband/stopband width in units of input sampling rate.
:param order: highest order, defaults to :param:`rate`"""
if order is None:
order = rate
coeff = []
p = 1
while p < rate:
p *= 2
coeff.append(halfgen4(width*p/rate/2, order*p//rate))
return coeff


class ParallelHBFUpsampler(Module):
"""Parallel, power-of-two, half-band, cascading upsampler.
Coefficients should be normalized to overall gain of 2
(highest/center coefficient being 1)."""
def __init__(self, coefficients, width=16, **kwargs):
self.parallelism = 1
self.latency = 0
self.width = width
self.i = Signal((width, True))

###

i = [self.i]
for coeff in coefficients:
self.parallelism *= 2
# assert coeff[len(coeff)//2 + 1] == 1
hbf = ParallelFIR(coeff, self.parallelism, width, **kwargs)
self.submodules += hbf
self.comb += [a.eq(b) for a, b in zip(hbf.i[::2], i)]
i = hbf.o
self.latency += hbf.latency
self.o = i
27 changes: 17 additions & 10 deletions artiq/gateware/dsp/sawg.py
Original file line number Diff line number Diff line change
@@ -4,9 +4,10 @@
from misoc.interconnect.stream import Endpoint
from misoc.cores.cordic import Cordic

from .accu import PhasedAccu, Accu
from .accu import PhasedAccu
from .tools import eqh, Delay, SatAddMixin
from .spline import Spline
from .fir import ParallelHBFUpsampler, halfgen4_cascade


_Widths = namedtuple("_Widths", "t a p f")
@@ -145,13 +146,17 @@ def __init__(self, width=16, parallelism=4, widths=None, orders=None):

self.submodules.a1 = a1 = SplineParallelDDS(widths, orders)
self.submodules.a2 = a2 = SplineParallelDDS(widths, orders)
coeff = [[int(round((1 << 26) * ci)) for ci in c]
for c in halfgen4_cascade(parallelism, width=.4, order=8)]
hbf = [ParallelHBFUpsampler(coeff, width=width, shift=25)
for i in range(2)]
self.submodules.b = b = SplineParallelDUC(
widths._replace(a=len(a1.xo[0]), f=widths.f - width), orders,
parallelism=parallelism, a_delay=-a1.latency)
parallelism=parallelism, a_delay=-a1.latency-hbf[0].latency)
cfg = Config(widths.a)
u = Spline(width=widths.a, order=orders.a)
du = Delay(width, a1.latency + b.latency - u.latency)
self.submodules += cfg, u, du
du = Delay(width, a1.latency + hbf[0].latency + b.latency - u.latency)
self.submodules += cfg, u, du, hbf
self.u = u.tri(widths.t)
self.i = [cfg.i, self.u, a1.a, a1.f, a1.p, a2.a, a2.f, a2.p, b.f, b.p]
self.i_names = "cfg u a1 f1 p1 a2 f2 p2 f0 p0".split()
@@ -174,12 +179,14 @@ def __init__(self, width=16, parallelism=4, widths=None, orders=None):
Cat(a1.clr, a2.clr, b.clr).eq(cfg.clr),
]
self.sync += [
b.i.x.eq(self.sat_add(a1.xo[0], a2.xo[0],
limits=cfg.limits[0],
clipped=cfg.clipped[0])),
b.i.y.eq(self.sat_add(a1.yo[0], a2.yo[0],
limits=cfg.limits[1],
clipped=cfg.clipped[1])),
hbf[0].i.eq(self.sat_add(a1.xo[0], a2.xo[0],
limits=cfg.limits[0],
clipped=cfg.clipped[0])),
hbf[1].i.eq(self.sat_add(a1.yo[0], a2.yo[0],
limits=cfg.limits[1],
clipped=cfg.clipped[1])),
b.i.x.eq(hbf[0].o[0]), # FIXME: rip up
b.i.y.eq(hbf[1].o[0]),
eqh(du.i, u.o.a0),
]
# wire up outputs and q_{i,o} exchange
103 changes: 103 additions & 0 deletions artiq/test/gateware/test_fir.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
import numpy as np
import matplotlib.pyplot as plt

from migen import *
from migen.fhdl import verilog
from artiq.gateware.dsp import fir


class Transfer(Module):
def __init__(self, dut):
self.submodules.dut = dut

def drive(self, x):
for xi in x:
yield self.dut.i.eq(int(xi))
yield

def record(self, y):
for i in range(self.dut.latency):
yield
for i in range(len(y)):
yield
y[i] = (yield self.dut.o)

def run(self, samples, amplitude=1.):
w = 2**(self.dut.width - 1) - 1
x = np.round(np.random.uniform(
-amplitude*w, amplitude*w, samples))
y = np.empty_like(x)
run_simulation(self, [self.drive(x), self.record(y)],
vcd_name="fir.vcd")
x /= w
y /= w
return x, y

def analyze(self, x, y):
fig, ax = plt.subplots(3)
ax[0].plot(x, "c-.", label="input")
ax[0].plot(y, "r-", label="output")
ax[0].legend(loc="right")
ax[0].set_xlabel("time (1/fs)")
ax[0].set_ylabel("signal")
n = len(x)
w = np.hanning(n)
x = (x.reshape(-1, n)*w).sum(0)
y = (y.reshape(-1, n)*w).sum(0)
t = (np.fft.rfft(y)/np.fft.rfft(x))
f = np.fft.rfftfreq(n)*2
fmin = f[1]
ax[1].plot(f, 20*np.log10(np.abs(t)), "r-")
ax[1].set_ylim(-70, 3)
ax[1].set_xlim(fmin, 1.)
# ax[1].set_xscale("log")
ax[1].set_xlabel("frequency (fs/2)")
ax[1].set_ylabel("magnitude (dB)")
ax[1].grid(True)
ax[2].plot(f, np.rad2deg(np.angle(t)), "r-")
ax[2].set_xlim(fmin, 1.)
# ax[2].set_xscale("log")
ax[2].set_xlabel("frequency (fs/2)")
ax[2].set_ylabel("phase (deg)")
ax[2].grid(True)
return fig


class ParallelTransfer(Transfer):
def drive(self, x):
for xi in x.reshape(-1, self.dut.parallelism):
yield [ij.eq(int(xj)) for ij, xj in zip(self.dut.i, xi)]
yield

def record(self, y):
for i in range(self.dut.latency):
yield
for yi in y.reshape(-1, self.dut.parallelism):
yield
yi[:] = (yield from [(yield o) for o in self.dut.o])


def _main():
coeff = fir.halfgen4(.4/2, 8)
coeff_int = [int(round(c * (1 << 16 - 1))) for c in coeff]
if False:
coeff = [[int(round((1 << 26) * ci)) for ci in c]
for c in fir.halfgen4_cascade(8, width=.4, order=8)]
dut = fir.ParallelHBFUpsampler(coeff, width=16, shift=25)
print(verilog.convert(dut, ios=set([dut.i] + dut.o)))
elif True:
dut = fir.ParallelFIR(coeff_int, parallelism=4, width=16)
# print(verilog.convert(dut, ios=set(dut.i + dut.o)))
tb = ParallelTransfer(dut)
else:
dut = fir.FIR(coeff_int, width=16)
# print(verilog.convert(dut, ios={dut.i, dut.o}))
tb = Transfer(dut)

x, y = tb.run(samples=1 << 10, amplitude=.8)
tb.analyze(x, y)
plt.show()


if __name__ == "__main__":
_main()