为什么这个过滤后的噪声信号不像我预期的那样表现？
Why doesn't this filtered noise signal behave like I expect it should?
我正在为我的学士论文创建一个工具包。我编写了一些函数，但其中一个函数的行为与我期望的不同。
首先我创建了下面的函数。它会产生随机白噪声，将其绘制在电源频率域 "assuming" 中，它已经在该域中（用于模拟目的）。然后，应用逆傅立叶变换以获得模拟信号，然后应用常规傅立叶变换以获得我自己创建的白噪声。最后一步是验证函数的行为是否如我所料，确实如此。
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 沿 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 随机数！）。接下来，concatenate 行生成频谱的另一半，这是前半部分的镜像和复共轭版本。我将 0 频率和奈奎斯特频率的 bin 设置为零以简化操作。这两个应该是实值的。如果预期信号具有零均值，则 0 频率应为 0。
请注意，返回的 whitenoise_signal 实际上不是白色的，因为它的光谱已被过滤。它是粉红噪声，它在较低的频率具有较高的能量。
但还要注意 whitenoise_signal 是实数。 ifft returns 复数，但虚部确实接近于零（由于 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,:]))
我正在为我的学士论文创建一个工具包。我编写了一些函数，但其中一个函数的行为与我期望的不同。
首先我创建了下面的函数。它会产生随机白噪声，将其绘制在电源频率域 "assuming" 中，它已经在该域中（用于模拟目的）。然后，应用逆傅立叶变换以获得模拟信号，然后应用常规傅立叶变换以获得我自己创建的白噪声。最后一步是验证函数的行为是否如我所料，确实如此。
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 沿 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 随机数！）。接下来，concatenate 行生成频谱的另一半，这是前半部分的镜像和复共轭版本。我将 0 频率和奈奎斯特频率的 bin 设置为零以简化操作。这两个应该是实值的。如果预期信号具有零均值，则 0 频率应为 0。
请注意，返回的 whitenoise_signal 实际上不是白色的，因为它的光谱已被过滤。它是粉红噪声，它在较低的频率具有较高的能量。
但还要注意 whitenoise_signal 是实数。 ifft returns 复数，但虚部确实接近于零（由于 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,:]))