Performance concerns

[bradleyharden]

Using this library was pretty straight forward. I was able to do what I wanted, but the performance left a lot to be desired.

Here's what I was trying to do. Maybe there's a faster way to implement it? If not, I think I have a proposed solution.

I am using an FPGA and high-speed DAC to synthesize some sinusoidal tones. We are aiming for a SFDR of 70+ dB. The DAC update frequency is 30 MHz. Right now, I'm using a linear interpolation scheme to reduce the size of my sine look-up tables. I have a 28-bit phase accumulator, and I have two LUTS containing 18-bit, signed, fixed point slopes and offsets. I use the upper bits of the phase accumulator to address the LUTs, and I use 17 of the lower bits as an unsigned "distance". I then multiple the distance by the slope and add the offset. I send the resulting sample stream to the DAC.

We're seeing some spurious peaks in the FFT of the DAC output. I wanted to do some analysis to see what level of spurious peaks I should expect from my linear interpolation scheme. That's where the fixedpoint library comes in. I dumped my slope and offset LUTs from my HDL simulator and loaded them into Python as lists of FixedPoint values. Then, I calculated my fixed point distances and performed the same multiplication and addition.

The problem is that I want my simulated FFT to match the real FFT. The real FFT uses 10 ms of data, so I need to generate 300,000 FixedPoint results.

Right now, this is the best I've come up with. dist and addr are numpy arrays of integers. I couldn't find a way to initialize the raw fractional bits of dist_fp directly from the integers in dist, so I had to convert from int to str and back. However, creating dist_fp "only" takes 5-6 seconds. Computing the samples takes something like 30+ seconds.

  dist_fp = [FixedPoint(hex(d), True, 1, FRAC_BITS) for d in dist]
  samples = []
  for a, d in zip(addr, dist_fp):
      offset = offsets[a]
      slope = slopes[a]
      samples.append(offset + slope * d)

Is there a better way to do this? My gut tells me I might be able to improve it slightly, but I think the better solution would be to offer something like a FixedPointArray that is backed by numpy arrays underneath. I think it is common to have arrays of fixed point numbers with identical types, so I think it would be useful. You would probably have to limit the total number of bits to 64, and use dtype="u8" under the hood, but that doesn't seem like a major limitation to me.

Thoughts?

[ZackSheffieldSEL]

I couldn't find a way to initialize the raw fractional bits of dist_fp directly from the integers in dist, so I had to convert from int to str and back.

That is correct, this is the fastest way to do it.

Is there a better way to do this?

When you do offset + slope * d, it actually creates 2 FixedPoint objects, which is expensive. one of which doesn't get used. First, the product of slope * d is generated, then the sum of offset + <product> is created. By minimizing the number of objects created, you'll save some time.

I'm able to trim off a few seconds off the loop by doing the following:

for a, d in zip(addr, dist_fp):
    val = FixedPoint(slopes[a])  # simple/quick copy, no error checking
    # val = slopes[a]  # even quicker if you don't mind changing the value of slopes[a]
    val *= d
    val += offsets[a]
    samples.append(val)

The reason this is quicker is because

• val is simply copied, and no error checking is performed (see FixedPoint.__init__
• Augmented operators += and *= do not create new FixedPoint objects; the result of the arithmetic operation is simply assigned to the existing val object.

So when I use tqdm to measure iteration time I get the following for your example above:
300000it [00:26, 11442.99it/s]

and this for my augmented operator example using val = FixedPoint(slopes[a]):
300000it [00:19, 15722.85it/s]

and this is for my augmented operator example using val = slopes[a]:
300000it [00:15, 19874.71it/s]

Here's my full test script:

from fixedpoint import FixedPoint
import random
from tqdm import tqdm

N = 300000
FRAC_BITS = 17
# I have a 28-bit phase accumulator, and I have two LUTS containing 18-bit,
# signed, fixed point slopes and offsets.
# I use the upper bits of the phase accumulator to address the LUTs
addr = range(N)
# and I use 17 of the lower bits as an unsigned "distance".
dist = [random.getrandbits(17) for _ in range(N)]
# I dumped my slope and offset LUTs from my HDL simulator and loaded them into
# Python as lists of FixedPoint values.
# I have two LUTS containing 18-bit, signed, fixed point slopes and offsets.
offsets = [FixedPoint(hex(random.getrandbits(18)), 1, 1, FRAC_BITS) for _ in addr]
slopes = [FixedPoint(hex(random.getrandbits(18)), 1, 1, FRAC_BITS) for _ in addr]

# Code sample
dist_fp = [FixedPoint(hex(d), True, 1, FRAC_BITS) for d in dist]
samples = []
for a, d in tqdm(zip(addr, dist_fp)):
    offset = offsets[a]
    slope = slopes[a]
    samples.append(offset + slope * d)

for a, d in tqdm(zip(addr, dist_fp)):
    val = FixedPoint(slopes[a])  # simple/quick copy, no error checking
    val *= d
    val += offsets[a]
    samples.append(val)

for a, d in tqdm(zip(addr, dist_fp)):
    val = slopes[a]  # quicker if you don't mind changing the value of slopes[a]
    val *= d
    val += offsets[a]
    samples.append(val)

I think the better solution would be to offer something like a FixedPointArray that is backed by numpy arrays underneath. I think it is common to have arrays of fixed point numbers with identical types, so I think it would be useful. You would probably have to limit the total number of bits to 64, and use dtype="u8" under the hood, but that doesn't seem like a major limitation to me.

This seems like a good idea; I am not as familiar with numpy as I should be but this is definitely worth looking into.

[bradleyharden]

Thanks for the tips. I think the problem is that I want to avoid having to micro-optimize my code. I think that would quickly become tedious for more complicated numerical applications. It also requires deeper knowledge of the FixedPoint class. You have to consider all possible ways to implement your task, and you have to know the cost of each operation within that set before you can determine the approach with the best performance. I think all of that would drop away with a FixedPointArray.