python - 为什么这个经过过滤的噪声信号的表现不像我预期的那样？

Question

我正在为我的学士论文创建一个工具包。我编写了许多函数，但其​​中一个函数的行为与我预期的不同。

<小时/>

首先我创建了下面的函数。它会产生随机白噪声，将其绘制在电源频域中，“假设”它已经位于该域中(用于模拟目的)。然后，应用逆傅里叶变换来获得模拟信号，然后应用常规傅里叶变换来获得我自己创建的白噪声。最后一步是验证该函数的行为是否符合我的预期，而且确实如此。

def white_noise(n: int, N: int, slope: int = grad):
    x = np.linspace(1, 100, n)
    slope_loglog = (10 ** (slope * np.log10(x) + 1))

    whitenoise = rnd.rand(n, N)
    whitenoise_power = whitenoise ** 2  # quadratic of the white noise to retrieve the power spectrum

    whitenoise_filtered = (whitenoise_power.T * slope_loglog).T
    whitenoise_signal = fft.ifft(whitenoise_filtered)
    whitenoise_retransformed = fft.fft(whitenoise_signal)

    return whitenoise, whitenoise_filtered, whitenoise_signal, whitenoise_retransformed, slope_loglog

此后，我通过绘制结果来检查生成的白噪声和双重变换的白噪声是否相同。如下图所示，它们是相同的，从而验证了我的脚本是否有效。

<小时/>

现在，我的问题出现了。在上面脚本的修改版本中(见下文)。我生成的白噪声和双重变换的白噪声的行为不同。修改后的脚本添加了 modified_roll 的功能(一个小函数，将函数滚动到自身上以模拟随时间变化的干扰)。

def white_noise_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    whitenoise, whitenoise_filtered, whitenoise_signal = white_noise(n, N)[:3]

    # only showing the selected arrays
    arrays_to_select = random_arrays(N, num)
    selected_whitenoise = whitenoise[:, arrays_to_select].copy()
    selected_whitenoise_filtered = whitenoise_filtered[:, arrays_to_select].copy()
    selected_whitenoise_signal = whitenoise_signal[:, arrays_to_select].copy()

    # shifting the signal as a field of different refractive index would do
    if operations == 0:
        shifted_signal = selected_whitenoise_signal
    else:
        shifted_signal = modified_roll(selected_whitenoise_signal.copy(), shift, operations)

    # fourier transform back to the power frequency domain
    shifted_whitenoise = fft.fft(shifted_signal)

    return selected_whitenoise, selected_whitenoise_filtered, selected_whitenoise_signal, shifted_signal, \
           shifted_whitenoise

正如人们所见，填写 white_noise_signal_shift(n, N, N, 0, 0)这样应该等于 white_noise(n, N) (假设您使用的是 numpy.random.seed()) 。但是，对于大 N 则不然，如下图所示。在这两个脚本中，返回功率频域的步骤是 fft.fft("signal")可以看到图中的信号也是相同的。我做错了什么？

<小时/>

复制此内容以获取第二个数字的结果。

import matplotlib.pyplot as plt
import numpy as np
import numpy.fft as fft
import numpy.random as rnd

grad = -5/3.

def white_noise(n: int, N: int, slope: int = grad):
    x = np.linspace(1, 100, n)
    slope_loglog = (10 ** (slope * np.log10(x) + 1))

    whitenoise = rnd.rand(n, N)
    whitenoise_power = whitenoise ** 2  # quadratic of the white noise to retrieve the power spectrum

    whitenoise_filtered = (whitenoise_power.T * slope_loglog).T
    whitenoise_signal = fft.ifft(whitenoise_filtered)
    whitenoise_retransformed = fft.fft(whitenoise_signal)

    return whitenoise, whitenoise_filtered, whitenoise_signal, whitenoise_retransformed, slope_loglog

# random array selection
def random_arrays(N: int, num: int):
    res = np.asarray(range(N))
    rnd.shuffle(res)
    return res[:num]


def modified_roll(inp, shift: int, operations: int):
    count = 0
    array = inp[:]
    array_rolled = array.copy()
    for k in range(operations):
        count += shift
        array = np.roll(array, shift, axis=0)
        array[:count] = 0
        array_rolled += array

    out = array_rolled / operations
    return out

def white_noise_signal_shift(n: int, N: int, num: int, shift: int, operations: int):
    whitenoise, whitenoise_filtered, whitenoise_signal = white_noise(n, N)[:3]

    # only showing the selected arrays
    arrays_to_select = random_arrays(N, num)
    selected_whitenoise = whitenoise[:, arrays_to_select].copy()
    selected_whitenoise_filtered = whitenoise_filtered[:, arrays_to_select].copy()
    selected_whitenoise_signal = whitenoise_signal[:, arrays_to_select].copy()

    # shifting the signal as a field of different refractive index would do
    if operations == 0:
        shifted_signal = selected_whitenoise_signal
    else:
        shifted_signal = modified_roll(selected_whitenoise_signal.copy(), shift, operations)

    # fourier transform back to the power frequency domain
    shifted_whitenoise = fft.fft(shifted_signal)

    return selected_whitenoise, selected_whitenoise_filtered, selected_whitenoise_signal, shifted_signal, \
           shifted_whitenoise

# this plots white_noise_signal_shift
def plt_white_noise_signal_shift(n: int, N: int, num_ar: int, shift, operations, size=(10, 7.5)):

    whitenoise, whitenoise_filtered, whitenoise_signal, shifted_signal, shifted_whitenoise \
        = white_noise_signal_shift(n, N, num_ar, shift, operations)

    fig = plt.figure(figsize=size)

    ax1 = plt.subplot2grid((3, 2), (0, 0), rowspan=1, colspan=2)
    ax2 = plt.subplot2grid((3, 2), (1, 0), rowspan=1, colspan=1)
    ax3 = plt.subplot2grid((3, 2), (2, 0), rowspan=1, colspan=1)
    ax4 = plt.subplot2grid((3, 2), (1, 1), rowspan=1, colspan=1, sharey=ax2)
    ax5 = plt.subplot2grid((3, 2), (2, 1), rowspan=1, colspan=1, sharey=ax3)

    ax1.set_title('1) Original white noise')
    ax2.set_title('2) Filtered original white noise'), ax2.set_ylabel('Log(P)'), ax2.set_xlabel('Log(f)')
    ax3.set_title('3) Original signal')
    ax4.set_title('5) White noise from the shifted signal'), ax4.set_ylabel('Log(P)'), ax4.set_xlabel('Log(f)')
    ax5.set_title('4) Shifted signal')

    # plotting the whitenoise
    ax1.plot(whitenoise)

    # plotting white the original data
    ax2.loglog(whitenoise_filtered)
    ax3.plot(whitenoise_signal)

    # plotting the shifted data
    ax4.loglog(shifted_whitenoise)
    ax5.plot(shifted_signal)

    plt.tight_layout()
    plt.show()


rnd.seed(50)

# to run the script
plt_white_noise_signal_shift(100, 50, 50, 0, 0)

Answer 1

当您计算 FFT 时，您应该始终明确指定您打算沿哪个轴进行计算。默认情况下，numpy.fft.fft 沿last 轴进行计算，但在您的代码中，您应该使用第一个轴。对于对 fft 和 ifft 的所有调用，请添加 axis 参数:

whitenoise_signal = fft.ifft(whitenoise_filtered, axis=0)

<小时/>

代码中的另一个问题是，您没有考虑实值信号预期的频谱对称性，并且一开始就没有生成复值频谱。这会导致复值时间信号，我认为这不是您想要的。您可以生成白噪声并以这种方式过滤它:

def white_noise(n: int, N: int, slope: int = grad):
    x = np.linspace(1, 100, n//2)
    slope_loglog = (10 ** (slope * np.log10(x) + 1))

    whitenoise = rnd.randn(n//2, N) + 1j * rnd.randn(n//2, N)
    whitenoise[0, :] = 0  # zero-mean noise
    whitenoise_filtered = whitenoise * slope_loglog[:, np.newaxis]

    whitenoise = np.concatenate((whitenoise, whitenoise[0:1, :], np.conj(whitenoise[-1:0:-1, :])), axis=0)
    whitenoise_filtered = np.concatenate((whitenoise_filtered, whitenoise_filtered[0:1, :], np.conj(whitenoise_filtered[-1:0:-1, :])), axis=0)

    whitenoise_signal = fft.ifft(whitenoise_filtered, axis=0)
    whitenoise_signal = np.real_if_close(whitenoise_signal)
    if np.iscomplex(whitenoise_signal).any():
        print('Warning! whitenoise_signal is complex-valued!')
    whitenoise_retransformed = fft.fft(whitenoise_signal, axis=0)

    return whitenoise, whitenoise_filtered, whitenoise_signal, whitenoise_retransformed, slope_loglog

我首先生成了信号长度一半的复值正态分布噪声(这仍然是n个随机数!)。接下来，连接线生成频谱的另一半，它是前半部分的镜像和复共轭版本。为了简化事情，我将 0 频率和奈奎斯特频率处的 bin 设置为零。这两个应该是实值。如果预期信号的均值为零，则 0 频率应为 0。

请注意，返回的whitenoise_signal实际上不是白色的，因为它的频谱已被过滤。它是粉红噪声，它在较低频率下具有较高能量。

但还要注意，whitenoise_signal 是实值。 ifft 返回复数，但虚部实际上接近于零(由于 FFT 计算中的舍入误差，不完全为零)。 np.real_if_close 丢弃虚部，因为它在 0 的小容差范围内。

要绘制频谱，请使用np.abs:

ax2.plot(np.abs(whitenoise_filtered))
ax2.set_yscale('log')

您也不应该对 x 轴应用对数缩放，因为这对于对称频谱来说看起来很有趣。如果您确实想这样做，您应该只绘制频谱的一半:

ax2.loglog(np.abs(whitenoise_filtered[0:n//2,:]))