Нулевое дополнение в дискретном ряду Фурье: индексация, DC-бин и корректная частотная шкала

Почему нулевое дополнение спектра в ДРФ даёт комплексный, смещённый сигнал. Разбираем индексацию, DC-бин и шкалу 2D, показываем корректное растяжение и 5 копий.

Нулевое дополнение спектра, чтобы растянуть дискретный сигнал, кажется делом простым — пока не подводит арифметика индексов. Ниже — компактное объяснение, почему вставка нулей между коэффициентами рядов Фурье привела к комплексной, смещённой форме волны вместо пяти аккуратных копий реального прямоугольного сигнала, и как это исправить, правильно настроив индексацию и частотную шкалу.

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

Начинаем с базового прямоугольного сигнала a[n] на диапазоне индексов [-1000, 1000], равного 1 при |n| < 100 и 0 вне этого диапазона. Далее вычисляются его коэффициенты дискретного ряда Фурье. Затем «растягиваем» спектр: после каждого коэффициента вставляем четыре нуля и масштабируем на 0.2, то есть 0.2·[a0, 0, 0, 0, 0, a1, 0, 0, 0, 0, a2, …]. Ожидалось, что во времени получится реальный сигнал из пяти одинаковых копий исходного, разнесённых по оси n.

Код, демонстрирующий проблему

Ниже приведён фрагмент, который воспроизводит ситуацию: обратное преобразование возвращает комплексный сигнал, а копии оказываются сдвинутыми. Имена переменных условные, но вычислительная логика соответствует описанному подходу.

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Прямоугольный сигнал: 1 при |n| < 100
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# Подрезаем микрошум с плавающей точкой и «почти целые» значения
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# Прямое/обратное преобразование дискретного ряда Фурье (ошибочная нормировка)
def fs_transform_bad(data, halfspan, inv=False):
    win_len = 2 * halfspan + 1  # используется и как длина массива, и как частотная шкала
    out = np.zeros(win_len, dtype=complex)
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / win_len)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / win_len) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / win_len)
    return out
# Спектр прямоугольника
spec_a = fs_transform_bad(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# Растягиваем по частоте вставкой нулей и масштабированием
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
# Восстанавливаем тем же преобразованием
guess_L = (len(spec_f) - 1) // 2
time_f = fs_transform_bad(spec_f, guess_L, inv=True)
time_f = snap_near_integers(time_f)

Что пошло не так и почему это происходит

Корень проблемы — несоответствие индексации между D и N, из‑за которого «бин» DC смещается, а вместе с ним уезжает и спектральное содержимое. Если исходный полуинтервал D=L и длина массива N=2D+1, то непосредственное умножение N на фактор может сделать новую длину чётной или, даже если она остаётся нечётной, увести «середину» от DC. В нашем случае с фактором 5 середина у массива формально есть, но после вставки с шагом нулей коэффициент нулевой частоты не попадает в настоящий центр. Нулевая частота, которая раньше находилась в индексе L, сдвигается на пару позиций — это и вызывает смещение восстановленного сигнала и появление мнимой составляющей.

Это видно напрямую при просмотре дополнённого спектра: центральный индекс растянутого массива — не то место, куда попадает DC после вставки нулей, поэтому обратное преобразование берёт неверно выровненный DC-бин. Этого сдвига достаточно, чтобы объяснить фазовый снос и неожиданную мнимую часть, хотя исходный a[n] — строго вещественный и чётный.

Есть и второе, тесно связанное, нарушение в реализации дискретного ряда Фурье. Длина массива равна 2D+1, однако нормирующий множитель в показателе экспоненты должен строиться по 2D, а не по 2D+1. Использование 2D+1 одновременно как длины массива и как знаменателя в аргументе экспоненты смещает частотную сетку, а при последующем растяжении усиливает ошибку. Размер массива остаётся 2D+1, но нормировка в показателе должна быть 2D.

Исправление: растягивайте D, а не N, и используйте верную частотную шкалу

Две точечные правки решают проблему и дают пять вещественных, корректно выровненных копий сигнала. Во‑первых, сохраняйте инвариант N = 2D + 1. При растяжении умножайте D на фактор и уже из нового D восстанавливайте N. Во‑вторых, в коде дискретного ряда Фурье используйте 2D как нормировку в показателе экспоненты, даже если в массиве по‑прежнему 2D+1 отсчёт. После этих поправок коэффициент DC оказывается ровно в центральном индексе, и обратное преобразование возвращает ожидаемый вещественный рисунок.

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Прямоугольный сигнал
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# Чистка численного шума
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# Корректный дискретный ряд Фурье: размер массива 2D+1, частотная шкала — 2D
def fs_transform_ok(data, halfspan, inv=False):
    freq_span = 2 * halfspan            # шкала для показателя экспоненты
    out = np.zeros(freq_span + 1, dtype=complex)  # длина массива по-прежнему 2D+1
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / freq_span)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / freq_span) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / freq_span)
    return out
# Прямое преобразование
spec_a = fs_transform_ok(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# Растягивание: вставляем нули между коэффициентами, масштабируем на фактор
factor = 5
L_new = factor * L
M_new = 2 * L_new + 1
spec_g = np.zeros(M_new, dtype=complex)
spec_g[::factor] = spec_a / factor
# Обратное преобразование
time_g = fs_transform_ok(spec_g, L_new, inv=True)
time_g = snap_near_integers(time_g)

Если хотите наглядно увидеть рассинхрон DC-бина в «сломленной» версии, быстрый тест всё покажет: центральная ячейка — не DC после вставки нулей; в примере с фактором 5 нужное значение оказывается на паре индексов в сторону.

# Демонстрация смещения DC в неверном подходе
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
mid = spec_f.shape[0] // 2
print(spec_f[mid])       # оказывается 0
print(spec_f[mid - 2])   # здесь лежит масштабированный DC

Почему это важно

Когда вы управляете рядами Фурье по индексам, связь между полуинтервалом D, длиной N = 2D + 1 и положением DC — это не формальность, а основа сетки, на которой выполняются частотные операции. Если растягивать, напрямую трогая N, и перепутать шкалу в экспоненте, весь спектр сдвигается; обратное преобразование вынуждено собирать сигнал с неправильной частотной сетки. Отсюда и неожиданная мнимая часть, и потеря чётких периодических копий во времени.

Итоги

Держите инвариант N = 2D + 1 и растягивайте D, а не N. В дискретном ряде Фурье используйте 2D в нормировке показателя экспоненты, пока сигнал всё ещё имеет 2D + 1 отсчёт. После вставки нулей между спектральными отсчётами убедитесь, что коэффициент DC попадает в центр нового спектра. Соблюдая эти три правила, на выходе получите ожидаемую вещественную, повторяющуюся форму без лишних фазовых артефактов. Если понадобится ускорение, можно векторизовать внутренние суммы, но корректность здесь целиком про индексацию и частотную шкалу.