Почему этот отфильтрованный шумовой сигнал не ведет себя так, как я ожидал?

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

Когда вы вычисляете БПФ, вы всегда должны явно указывать, по какой оси вы собираетесь это делать. По умолчанию 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 и ячейку на частоте Найквиста я установил на ноль, чтобы упростить вещи. Эти два должны быть реальными. Частота 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')

Не следует также применять масштабирование журнала к оси абсцисс, так как это будет смешно с частотным спектром симметрии c. Если вы действительно хотите сделать это, вы должны построить только половину спектра:

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