Rounding really small/near zero complex values in Python

Question

Consider the following sequence/signal:

import numpy as np
from scipy.fft import fft

# %% Discrete data

x = np.array([0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875,
       0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875, 0.7966875])

I have defined the following function to calculate a Discrete Fourier Transform:

def DFT(x):
    N = len(x)
    n = np.arange(0, N)
    k = n.reshape((N, 1))
    omega = 2*np.pi*k/N
    e = np.exp(-1j*omega*n)
    X = np.dot(x, e)
    X_round = np.round(X, 6)
    return X, X_round, N, n

X, X_round, N, n = DFT(x)

Our X here have really small values for non-zero frequencies (as expected as the signal x does not change). Then, I have tried to round it using np.round(X, 6), that returns:

array([38.241+0.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, -0.   -0.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,  0.   -0.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, -0.   +0.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,  0.   -0.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,
        0.   +0.j,  0.   -0.j,  0.   -0.j])

But, if I try to take the angle of X_round, using X_round_angle = np.angle(X_round) it is not zero as I would expect (you can try using the fft function from scipy). Here is the output of X_round_angle:

array([ 0.        ,  3.14159265, -3.14159265, -3.14159265, -0.        ,
        0.        , -3.14159265, -0.        ,  3.14159265, -3.14159265,
       -3.14159265, -3.14159265,  3.14159265, -0.        , -3.14159265,
       -3.14159265, -3.14159265, -0.        , -0.        , -0.        ,
        0.        ,  0.        ,  3.14159265,  3.14159265,  0.        ,
       -3.14159265,  0.        ,  3.14159265, -3.14159265, -3.14159265,
       -3.14159265, -0.        , -0.        , -0.        , -3.14159265,
       -3.14159265, -0.        , -0.        , -0.        ,  0.        ,
        3.14159265, -3.14159265, -0.        ,  3.14159265, -3.14159265,
        0.        , -0.        , -0.        ])

Using np.angle(fft(x)) returns the correct phase value of zero for all frequencies. How can i correct my DFT function?

EDIT: Complementing the answer, I have found this link, that does the same:

When Im(X) and Re(X) are both equal to zero, the value of phase is undefined is undefined. It can be set equal to 0, as we will do below, to make the phase spectrum easier to read.

Answer 1

The behavior you observe is not totally unreasonable, since there is no way to define a phase angle for the complex number 0. np.angle() is giving you a mix of 0, pi, and -pi, all of which are equally correct.

The documentation for np.angle() says this:

Although the angle of the complex number 0 is undefined, numpy.angle(0) returns the value 0.

If we test this exactly as it is stated, we get the expected result of 0. However, the behavior that you're seeing happens when we tweak the parameter z from 0 (integer) to 0.0 (float). That's technically an undocumented edge case, so we can't expect any consistency there. We can test different combinations of positive and negative 0.0, with and without a 0.0j, and results are mixed.

That's all reasonable mathematically speaking, but there are certainly reasons why you might want the phase to be consistent for all kinds of zero (for example, for plotting). In that case, I'd simply use a comprehension to replace the offending entries with 0.

phase = np.angle(X_round)
phase = np.array([0.0 if np.abs(a)<0.000001 else theta for a, theta in zip(X_round, phase)])