python - 为什么此滤波后的噪声信号的行为不像我预期的那样？

对于我的学士论文，我正在创建一个工具箱。我已经编写了许多功能，但是其中一个功能的运行不像我期望的那样。



首先，我创建了下面的函数。它会产生随机的白噪声，并假设它已经在该频域中（出于仿真目的），将其绘制在工频域中。然后，应用傅立叶逆变换获得模拟信号，然后应用常规傅立叶变换获得我自己创建的白噪声。最后一步是验证该功能是否按照我的预期进行，并且确实如此。

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)

最佳答案

在计算FFT时，应始终明确指定要沿着哪个轴进行计算。默认情况下，numpy.fft.fft沿最后一个轴进行计算，但是在您的代码中，您应该使用第一个。对于对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随机数！）。接下来，concatenate线生成频谱的另一半，这是前半部分的镜像复共轭版本。为简化起见，我将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,:]))

关于python - 为什么此滤波后的噪声信号的行为不像我预期的那样？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/59619261/