在Python中四舍五入非常小/接近零的复数值

Question

考虑以下这个序列/信号：

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])

我定义了一个函数来计算离散傅里叶变换（DFT）：

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)

在这里，我们的 X 对于非零频率的值非常小（这很正常，因为信号 x 没有变化）。然后，我尝试用 np.round(X, 6) 来四舍五入，结果是：

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])

但是，如果我尝试计算 X_round 的角度，使用 X_round_angle = np.angle(X_round)，结果并不是我预期的零（你可以试试使用 scipy 的 fft 函数）。这是 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.        ])

使用 np.angle(fft(x)) 可以返回所有频率的正确相位值为零。我该如何修正我的 DFT 函数呢？

编辑：补充一下答案，我找到这个链接，它也做了同样的事情：

当 Im(X) 和 Re(X) 都等于零时，相位值是未定义的。可以将其设为零，正如我们下面将要做的，这样相位谱会更容易阅读。

scipy 四舍五入信号处理频率分析数值稳定性离散傅里叶变换相位谱复数值

在Python中四舍五入非常小/接近零的复数值

1 个回答

撰写回答