Как получить лог-нормальные сигналы с заданной PSD с помощью IAAFT

Как получить лог-нормальные сигналы с целевой PSD (синий и фиолетовый шум) без разрушения спектра: IAAFT, ранговое отображение и лог-лог интерполяция.

Генерация лог-нормальных сигналов с возрастающей плотностью спектральной мощности на словах кажется простой: сгенерировать гауссов шум, придать ему нужный спектр и возвести в степень. На практике последний шаг ломает спектральную структуру для высокочастотных вариантов вроде синего или фиолетового шума. В этом руководстве разберёмся, почему так происходит, и как получить сигналы, которые одновременно соответствуют целевому одномерному распределению (например, лог-нормальному) и целевой PSD с помощью IAAFT.

Постановка задачи

Цель — синтезировать случайные последовательности, чьи значения подчиняются длиннохвостому распределению на положительной полуоси (например, лог-нормальному), при этом плотность спектральной мощности мала на низких частотах и растёт с частотой (например, синий шум с наклоном +10 дБ/декаду или фиолетовый с +20 дБ/декаду). Взять выборку из лог-нормального распределения несложно. Формовать гауссов шум под заданную PSD через FFT — тоже стандартная практика. Трудность возникает, когда эти шаги сочетают, применяя нелинейное отображение (экспоненту) к уже сформованному процессу.

Демонстрация: сначала формование, затем экспонента

Следующий код генерирует цветной гауссов шум с несколькими каноническими PSD, затем применяет экспоненту, чтобы получить лог-нормальные выборки. После возведения в степень форма PSD сохраняется для белого и розового шума, но для синего и фиолетового спектр деградирует к белоподобному.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
# Спектральное формование через rFFT: умножаем на целевой амплитудный спектр
# и выполняем обратное преобразование.
def shape_spectrum(n_points, psd_fn=lambda f: 1.0):
    freq_side = np.fft.rfft(np.random.randn(n_points))
    spec = psd_fn(np.fft.rfftfreq(n_points))
    spec = spec / np.sqrt(np.mean(spec ** 2))
    colored = freq_side * spec
    return np.fft.irfft(colored)
# Вспомогательная функция: «замораживает» функцию PSD в генератор вида длина->сигнал.
def with_psd(psd_fn):
    return lambda n: shape_spectrum(n, psd_fn)
@with_psd
def make_white(f):
    return 1.0
@with_psd
def make_blue(f):
    return np.sqrt(f)
@with_psd
def make_violet(f):
    return f
@with_psd
def make_brown(f):
    return 1.0 / np.where(f == 0, float('inf'), f)
@with_psd
def make_pink(f):
    return 1.0 / np.where(f == 0, float('inf'), np.sqrt(f))
def synthesize_colored(n_points, kind="white"):
    if kind == "white":
        return make_white(n_points)
    elif kind == "pink":
        return make_pink(n_points)
    elif kind == "blue":
        return make_blue(n_points)
    elif kind == "violet":
        return make_violet(n_points)
    elif kind == "brown":
        return make_brown(n_points)
# Демонстрация: сравниваем гауссовские и лог-нормализованные версии при разных типах PSD.
n_points = 100_000
kinds = ["white", "pink", "blue", "violet"]
colors = ["lightgray", "magenta", "blue", "violet"]
fig, axes = plt.subplots(2, 2, sharex=False, sharey="row")
bins_gauss = np.linspace(-5, 5, 101)
for nm, col in zip(kinds, colors):
    x = synthesize_colored(n_points, nm)
    axes[0, 0].hist(x, bins=bins_gauss, color=col, label=nm, histtype="step", density=True)
    fx, px = welch(x)
    axes[1, 0].loglog(fx, px, color=col, label=nm)
bins_logn = np.linspace(0, 10, 101)
for nm, col in zip(kinds, colors):
    y = np.exp(synthesize_colored(n_points, nm))
    axes[0, 1].hist(y, bins=bins_logn, color=col, label=nm, histtype="step", density=True)
    fy, py = welch(y)
    axes[1, 1].loglog(fy, py, color=col, label=nm)
axes[0, 0].set_title("Gaussian samples")
axes[0, 1].set_title("Exp(Gaussian) samples")
axes[0, 0].set_ylabel("Density")
axes[0, 0].set_xlabel("Value")
axes[0, 1].set_xlabel("Value")
axes[1, 0].set_ylabel("Power")
axes[1, 0].set_xlabel("Frequency [Hz]")
axes[1, 1].set_xlabel("Frequency [Hz]")
axes[0, 0].legend(loc="upper left")
fig.tight_layout()
plt.show()

Почему спектральная форма рушится после возведения в степень

Экспонента — нелинейная операция. В отличие от линейной фильтрации, она не сохраняет линейную корреляционную структуру, которую кодирует PSD. Для высокочастотных спектров, таких как синий и фиолетовый шум, применение экспоненты сильно искажает спектр, и итоговая PSD стремится к белому шуму. Это согласуется с наблюдаемым поведением.

Отдельный, независимый момент в том, что нормировка формующей кривой для сохранения средней энергии, как часто делают при генерации цветного шума в частотной области, не компенсирует последующую нелинейную дисторсию. Эффект «отбеливания» вызван именно нелинейным отображением, а не нормировкой.

IAAFT: одновременно соблюдаем PSD и маргинальное распределение

Надёжный способ выполнить оба ограничения — итеративно переходить от одного к другому. Алгоритм IAAFT (Iterative Amplitude Adjusted Fourier Transform) делает ровно это. Он чередует две проекции: на одном шаге навязывается требуемый модуль спектра; на другом восстанавливается целевое распределение значений через ранговое отображение. Такой подход описан для стохастических сигналов и хорошо задокументирован в профильной литературе.

Дополнительная деталь реализации, заметно улучшающая точность соответствия спектральной цели, — строить модуль ДПФ с использованием кусочно-линейной интерполяции в лог–лог шкале, когда PSD задана в лог–лог виде. Это соответствует тому, как обычно задают такие наклоны (например, «+10 дБ/декаду»), и даёт гораздо более плотный контроль над PSD.

import numpy as np
from scipy import signal as sig
from scipy.stats import rankdata
from numpy.fft import rfft, irfft
# Конфигурация
n_pts = 2 ** 16
fs_rate = 1.0
max_iter = 50
eps_tol = 1e-6  # порог относительной ошибки PSD
# Целевая маргиналь: лог-нормальное распределение через экспоненту от гауссовских выборок
target_seed = np.exp(np.random.randn(n_pts))
asc_template = np.sort(target_seed)
# Целевой модуль спектра: пример +10 дБ/декаду («синий шум»)
# Строим через кусочно-линейную интерполяцию на лог–лог осях.
brk_freqs = np.array([1e-5 * fs_rate, fs_rate / 2])
brk_psd_db = 10 * np.log10(brk_freqs)
freq_bins = np.fft.rfftfreq(n_pts, d=1 / fs_rate)
# Добавляем первую точку, чтобы избежать log10(0), и явно задаём DC.
brk_freqs = np.insert(brk_freqs, 0, fs_rate / n_pts)
brk_psd_db = np.insert(brk_psd_db, 0, brk_psd_db[0])
psd_interp_db = np.interp(
    np.log10(freq_bins[1:]),
    np.log10(brk_freqs[1:]),
    brk_psd_db[1:]
)
# Преобразуем в амплитудный спектр; DC-бин ставим в ноль для ВЧ-форм.
target_amp = np.empty_like(freq_bins)
target_amp[0] = 0.0
target_amp[1:] = 10 ** (psd_interp_db / 20)
# Ранговое отображение: заменяем значения на значения целевой маргинали с теми же рангами.
def remap_by_rank(values, sorted_ref):
    ranks = rankdata(values, method='ordinal') - 1
    return sorted_ref[ranks]
# Цикл IAAFT: поочерёдно навязываем спектр и восстанавливаем маргиналь.
synth = target_seed.copy()
prev_err = np.inf
for i in range(max_iter):
    # Навязываем целевой модуль спектра.
    X = rfft(synth)
    X = target_amp * np.exp(1j * np.angle(X))
    shaped = irfft(X, n=n_pts)
    # Восстанавливаем целевую маргиналь ранговым отображением.
    synth = remap_by_rank(shaped, asc_template)
    # Периодически проверяем сходимость.
    if (i % 5 == 0) or (i == max_iter - 1):
        cur_err = np.linalg.norm(np.abs(rfft(synth)) - target_amp) / np.linalg.norm(target_amp)
        print(f"iter {i:2d}: PSD rel-err = {cur_err:.3e}")
        if cur_err > prev_err - eps_tol:
            print("converged.")
            break
        prev_err = cur_err

Сравнение с подходами на основе фильтрации

Линейная фильтрация лог-нормальных выборок фильтрами с сине-шумовым формованием (FIR или IIR) может обеспечить нужные наклоны PSD. Результат остаётся строго положительным и с правой асимметрией, но уже не является лог-нормальным из-за линейного смешивания отсчётов. В приведённых результатах подход IAAFT дал наилучшее одновременное соответствие PSD и маргинали среди протестированных методов, плюс бонусом — лог–лог интерполяция для точного отслеживания PSD, заданных по декадам.

Тем не менее, тесты показали, что для лог-нормальной целевой маргінали и небелой PSD эмпирическое распределение может отклоняться и иметь моду в нуле, а не около единицы, хотя хвост остаётся длинным. Это важное предупреждение при валидации результатов.

Зачем это важно

Многие последующие алгоритмы предполагают определённые особенности PSD или маргинальные свойства. Если что-то не совпадает — например, PSD разваливается после нелинейного отображения или распределение уплывает после фильтрации — бенчмарки, детекторы или имитационные анализы могут вводить в заблуждение. Совместное соблюдение обеих сторон защищает от таких расхождений.

Практические выводы

Если вам нужны и заданная PSD, и конкретное распределение значений, например лог-нормальное, итеративная схема вроде IAAFT — эффективный инструмент. Она чередует навязывание спектра и восстановление амплитуд по рангам, а при построении PSD в лог–лог шкале позволяет точно отслеживать наклоны вроде +10 дБ/декаду. При валидации смотрите и на PSD, и на маргиналь. Помните, что для небелых целей эмпирическое распределение всё же может отклоняться; проверьте положение моды и поведение хвоста относительно ожиданий. Если важна только PSD, прямая фильтрация может быть проще; если критичны обе стороны, уместнее итеративный подход.

Вкратце: формуйте спектр и восстанавливайте амплитуды совместно, проверяйте оба представления данных и отдавайте предпочтение лог–лог обработке, когда PSD задана по декадам.