19

I would like to create numpy.ndarray objects that hold complex integer values in them. NumPy does have complex support built-in, but for floating-point formats (float and double) only; I can create an ndarray with dtype='cfloat', for example, but there is no analogous dtype='cint16'. I would like to be able to create arrays that hold complex values represented using either 8- or 16-bit integers.

I found this mailing list post from 2007 where someone inquired about such support. The only workaround they recommended involved defining a new dtype that holds pairs of integers. This seems to represent each array element as a tuple of two values, but it's not clear what other work would need to be done in order to make the resulting data type work seamlessly with arithmetic functions.

I also considered another approach based on registration of user-defined types with NumPy. I don't have a problem with going to the C API to set this up if it will work well. However, the documentation for the type descriptor strucure seems to suggest that the type's kind field only supports signed/unsigned integer, floating-point, and complex floating-point numeric types. It's not clear that I would be able to get anywhere trying to define a complex integer type.

What are some recommendations for an approach that may work?

Whatever scheme I select, it must be amenable to wrapping of existing complex integer buffers without performing a copy. That is, I would like to be able to use PyArray_SimpleNewFromData() to expose the buffer to Python without having to make a copy of the buffer first. The buffer would be in interleaved real/imaginary format already, and would either be an array of int8_t or int16_t.

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
Jason R
  • 11,159
  • 6
  • 50
  • 81
  • 2
    This is pretty - non-standard. How do you define division for this type? for example, what do you expect if you do `(2+1j)/(3+0j)`? Do you expect it to give you a complex result or `(0+0j)`? – mgilson Dec 13 '12 at 16:02
  • @mgilson: Understood. That particular case would have to be handled specially. I wouldn't say that it's nonstandard; complex integers are used very frequently in signal processing. Complex division is not a common operation in such contexts; I could live with that operation not being supported. – Jason R Dec 13 '12 at 17:27
  • I am not comfortable with the C API, so I would subclass `ndarray`, force a dtype holding pairs of integers, and overload all the arithmetic operations for the new subclass. – Jaime Dec 13 '12 at 17:56
  • 2
    Out of curiosity, when in signal processing are complex integers used? I can't think of examples off the top of my head. – acjay Dec 13 '12 at 17:57
  • Which operations do you need? I simple `ndarray` subclass might be good enough. – mgilson Dec 13 '12 at 18:18
  • 1
    @acjohnson55: One application is in digital communications. Communications signals are often represented in what's called *complex baseband* format (this is also known as forming an *analytic signal*), which makes many common manipulations easier to model and implement. – Jason R Dec 13 '12 at 18:21
  • 1
    @mgilson: Likely addition/subtraction and multiplication. I also just thought of another constraint that I'll edit into the post. – Jason R Dec 13 '12 at 18:23
  • 1
    @JasonR -- I unfortunately don't have time to work on something of this magnitude right now, but it seems that you could accomplish something like this via a subclass using [view-casting](http://docs.scipy.org/doc/numpy/user/basics.subclassing.html#view-casting). You could add properties `imag` and `real` which would return views into the appropriate portions of the arrays. And you could override `__mul__`,`__add__`,`__sub__` accordingly, but I don't know exactly how numpy-like you actually *need* this to be, so I can't say for sure. – mgilson Dec 13 '12 at 18:38
  • @JasonR: I'm familiar with complex baseband format, but I haven't seen a system that's restricted to complex integers, let alone where the abstraction would fail if you did floating point arithmetic. But no biggie, I was just curious. – acjay Dec 13 '12 at 18:40
  • @acjohnson55: Floating-point format is the more typical case, but I will occasionally work with data that is stored as integers. While most processing is done in floating-point, for some common operations, I would like to be able to avoid a type conversion. – Jason R Dec 13 '12 at 18:53
  • 1
    @EMS: I'm aware that most work is done in floating-point. However, I have a requirement to interface with libraries that use integer formats, so I would like to have the capability to, at a minimum, view and do some basic manipulations on complex integers without necessarily having to do the conversion to floating point. My desire is driven strictly by performance (i.e. speed). – Jason R Jan 28 '13 at 20:03
  • 5
    @EMS: I fully understand the disadvantages of using integers for complex arithmetic. However, your argument isn't very constructive to this problem; suffice it to say that I have a *requirement* to sometimes interface with complex data that is formatted as integers. – Jason R Jan 29 '13 at 13:12
  • 4
    Just to add to the discussion, I too have a requirement for complex integers. It's basically for modelling fixed point implementations of algorithms. This most certainly is _not_ rare. That said, in my case it can be worked around using floating point complex values with suitable rounding (since I'm only dealing with multiplication). – Henry Gomersall Jul 16 '14 at 17:27
  • 2
    I have a requirement for complex integers too. My data size is around 1TBso the conversion time and memory usage is an issue. – J. Martinot-Lagarde Jul 21 '14 at 09:25
  • All three links are broken. The first times out and the other two are *"404. Not Found"*. – Peter Mortensen Apr 07 '22 at 16:57
  • Related: *[Is there a way to make a complex number in NumPy with the accuracy of the 'Decimal' type?](https://stackoverflow.com/questions/21750012/)* – Peter Mortensen Apr 07 '22 at 17:00

3 Answers3

18

I also deal with lots of complex integer data, generally basebanded data.

I use

dtype = np.dtype([('re', np.int16), ('im', np.int16)])

It's not perfect, but it adequately describes the data. I use it for loading into memory without doubling the size of the data. It also has the advantage of being able to load and store transparently with HDF5.

DATATYPE  H5T_COMPOUND {
    H5T_STD_I16LE "re";
    H5T_STD_I16LE "im";
}

Using it is straightforward and is just different.

x = np.zeros((3,3),dtype)
x[0,0]['re'] = 1
x[0,0]['im'] = 2
x
>> array([[(1, 2), (0, 0), (0, 0)],
>>        [(0, 0), (0, 0), (0, 0)],
>>        [(0, 0), (0, 0), (0, 0)]],
>>  dtype=[('re', '<i2'), ('im', '<i2')])

To do math with it, I convert to a native complex float type. The obvious approach doesn't work, but it's also not that hard.

y = x.astype(np.complex64) # doesn't work, only gets the real part
y = x['re'] + 1.j*x['im']  # works, but slow and big
y = x.view(np.int16).astype(np.float32).view(np.complex64)
y
>> array([[ 1.+2.j,  0.+0.j,  0.+0.j],
>>        [ 0.+0.j,  0.+0.j,  0.+0.j],
>>        [ 0.+0.j,  0.+0.j,  0.+0.j]], dtype=complex64)

This last conversion approach was inspired by an answer to What's the fastest way to convert an interleaved NumPy integer array to complex64?

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
Greg Allen
  • 575
  • 5
  • 14
  • I convert my array `f` (which looks like that: `array([(127, -128), (127, 127)], dtype=[('I', 'i1'), ('Q', 'i1')])`) like that: `f_complex = [np.complex(*x) for x in f]`. Its still not a complex integer, but I think a more elegant solution to convert from int to complex float. – Stefan D. Dec 01 '15 at 18:48
  • @StefanD. I expected your method to be slower, but I was shocked to measure >10000x slower on 1M points. It also only works for 1D (instead of any array), and `f_complex` is a list, not an array. – Greg Allen Jan 08 '16 at 23:22
  • really creative. i wonder when you do this... is there a copy operation that happens when you use `y = x.view(np.int16).astype(np.float32).view(np.complex64)`, yes or no? if no. that's amazing IMO. – Trevor Boyd Smith Mar 12 '20 at 21:51
  • astype() does a type conversion (and copy), but view() does not. So that's one copy-and-convert operation. – Greg Allen Mar 13 '20 at 18:42
1

Consider using matrices of the form [[a,-b],[b,a]] as a stand-in for the complex numbers.

Ordinary multiplication and addition of matrices corresponds to addition an multiplication of complex numbers (this subring of the collection of 2x2 matrices is isomorphic to the complex numbers).

I think Python can handle integer matrix algebra.

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
  • 1
    Thanks. It's possible that this could work *in principle*, but the main motivating factor was for efficiency, to allow creation of a view onto existing complex integer data. I think this method would require preprocessing the complex values into the matrix form that you suggested. – Jason R Oct 23 '20 at 14:14
  • 1
    @JasonR The nice thing is that this would allow, for example, exact multiplication of $(23232+54546i)*(25436+123132132i)$ without any inexact arithmetic coming from the use of floats. –  Oct 23 '20 at 15:46
-1

Python, and hence NumPy, does support complex numbers. If you want complex integers, just use np.round or ignore the decimal part.

For example,

import numpy as np

# Create 100 complex numbers in a 1D array
a = 100*np.random.sample(100) + (100*np.random.sample(100)*1j)

# Reshape to a 2D array
np.round(a)
a.reshape(10, 10)

# Get the real and imaginary parts of a couple of x/y points as integers
print int(a[1:2].real)
print int(a[3:4].imag)
Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
Jay M
  • 3,736
  • 1
  • 24
  • 33