Интерполяция в частотной плоскости дискретного Фурье

[MGM]

Дано дискретное разложение Фурье:

${x_n} = \frac{1}{N}\sum\limits_{k = 0}^{N - 1} {{X_k}} {e^{\frac{{2\pi i}}{N}kn}}$
Возьмем две соседние по частоте собственные функции разложения:
${e^{\frac{{2\pi i}}{N}kn}}$ и ${e^{\frac{{2\pi i}}{N}\left( {k + 1} \right)n}}$
Можно ли как выразить коэффициент разложения ${X_{k + dk}}$ ( $ 0 < dk < 1$) несобственной функции ${e^{\frac{{2\pi i}}{N}\left( {k + dk} \right)n}}$ через коэффициенты разложения двух соседних собственных
${X_{k + dk}} = f\left( {{X_k},{X_{k + 1}}} \right)$
Проще говоря, есть ли методика интерполяции в частотном пространстве, или проще брать «ближайшую» собственную частоту?
И что является ближайшей частотой, учитывая, что частота обратно пропорциональна длине волны?

[B@R5uk]

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

Заметьте, что для синуса, не укладывающегося целое число раз в период DFT, величина хвостов, а так же соотношение рассматриваемых вами амплитуд (и фаз) будет зависеть от того, как этот синус сфазирован относительно концов отрезка.

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

Так что, уточните, пожалуйста, ваши условия: что известно про сигнал и используете ли вы оконные функции и какие?

-- 29.07.2025, 15:19 --

Есть метод Прони ещё. Применимость ограничена, но экспоненциально затухающие сигналы он хорошо должен пережёвывать.

[MGM]

B@R5uk в сообщении #1695759 писал(а):

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

Заметьте, что для синуса, не укладывающегося целое число раз в период DFT, величина хвостов, а так же соотношение рассматриваемых вами амплитуд (и фаз) будет зависеть от того, как этот синус сфазирован относительно концов отрезка.

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

Так что, уточните, пожалуйста, ваши условия: что известно про сигнал и используете ли вы оконные функции и какие?

-- 29.07.2025, 15:19 --


Есть метод Прони ещё. Применимость ограничена, но экспоненциально затухающие сигналы он хорошо должен пережёвывать.

Задача проекции под некоторым углом. Похоже на интерференцию точечного сигнала с разделением двумя оттверстиями. То есть прективная волна перпендикулярна а второе отверстие под углом к плоскости. и фазовая длина волны второго источника увеличивается. То, что было высокой частотой - становится более низкой.

Но по сути, пытаюсь сделать быстрое дискретное преобразование Радона. А посему сложные коррекционные вычисления. Знаю, как минимум три быстрых алгоритма. Но мне они кажутся слишком неестесвенными. А самые лучшие алгоритмы - которые простые.

[B@R5uk]

Если модель известна, то, наверное, лучше будет подогнать данные под модель.

[MGM]

B@R5uk в сообщении #1695762 писал(а):

Если модель известна, то, наверное, лучше будет подогнать данные под модель.

Проблема в том, что алгоритмы обратного Радона пубикуют исключительно в интегральной форме. Моя попытка просто перевести формулу из Вики в дискрет, как это ни странно привела меня к удивительному успеху. Синтезированные синограммы обычных тестовых картинок - типа портретов, с точностью до контраста (var) реконструируются с точностью до 35dB - где-то близко к повороту дискретной картинки на произвольный угол и обратно. Подозреваю, комерческие алгоритмы так и устроены. Но это не быстрый алгоритм. $O(N^3)$
А все быстрые, как раз включают интерполяцию уже в частотной плоскости. DeepSeek посоветовал брать ближайшее. Ну или линейную комбинацию. Последнее мне не кажется правдоподобным.

[MGM]

В продолжение. Вот результат некой интерполяции на дискрете.

Справа, то что должно получиться, слева - что получилось при обратном преобразовании Фурье (цвета - красный - самый яркий, темно синий = 0, сигнал норсализован по мкс-мин, так что информативны только границы). Ясно, что есть некие пробемы. Когда-то (80-90е) были люди, умевшие поставить диагноз. Понятно, что что то со знаками коэффициентов разложения в ряд Фурье. Видимо, некоторая нессиметрия. Никто не в курсе, как лечить?
Насколько помню, если $F\left( {u,v} \right) = \int\limits_{ - R}^R {\int\limits_{ - R}^R {f\left( {x,y} \right){e^{ - i\left( {ux + vy} \right)}}dxdy} } $, то ${\mathop{\rm Re}\nolimits} \left\{ {F\left( {u,v} \right)} \right\} = {\mathop{\rm Re}\nolimits} \left\{ {F\left( { - u, - v} \right)} \right\}$ и
${\mathop{\rm Im}\nolimits} \left\{ {F\left( {u,v} \right)} \right\} = - {\mathop{\rm Im}\nolimits} \left\{ {F\left( { - u, - v} \right)} \right\}$. Или это не так?
Да и вообще, может ошибка и не в этом.

[MGM]

Ошибка, если кому интересно, в периодичноской модели любого дискретного преобразования Фурье. То есть если интегрирование по собственной двумерной функции идет по направлениям условно не x или y, то нужно достраивать переодическое дополнение. Что невозможно. Но есть выход. Удвоить границы интегрирования и заполнить нулями. Наверное.

[B@R5uk]

MGM в сообщении #1704779 писал(а):

Справа, то что должно получиться, слева - что получилось при обратном преобразовании Фурье

Очевидно, что вы что-то напутали и обратное преобразование не соответствует прямому (один оператор не является обратным ко второму). Попробуйте построить матрицу прямого преобразования (БПФ — это линейное преобразование, его можно записать в матричном виде), обратить её в лоб (для какого-нибудь разумно малого размера данных) и сравнить результат с тем, что вы делаете при обращении по формуле (записав это обращение тоже в матричной форме).

MGM в сообщении #1695753 писал(а):

Можно ли как выразить коэффициент разложения ${X_{k + dk}}$ ( $ 0 < dk < 1$) несобственной функции ${e^{\frac{{2\pi i}}{N}\left( {k + dk} \right)n}}$ через коэффициенты разложения двух соседних собственных
${X_{k + dk}} = f\left( {{X_k},{X_{k + 1}}} \right)$

Если я правильно понял вопрос, то вот эта статья обсуждает его в одномерном случае: Eric Jacobsen and Peter Kootsookos - Fast, Accurate Frequency Estimators. Однако, это всё от лукавого. В реальности шумы и/или другие частотные компоненты, которые не попадают на гребёнку БПФ могут заметно сместить оценку частоты. Чем ближе к друг другу частотные компоненты, тем сильнее смещение. Даже для чистого синусоидального сигнала результат может зависеть от того, как он сфазирован относительно прямоугольного окна. Поэтому лучше использовать окна типа окно Ганна или окно flap-top, у которых на краях и сама оконная функция и её производная нулевые. Далее, реальные сигналы как правило имеют амплитудную и частотную модуляцию. В этом случае уточнение частоты вообще смысл теряет, потому что в разные моменты времени она может быть разная. Надеюсь, мой ответ поможет хоть как-то.

[realeugene]

MGM в сообщении #1704779 писал(а):

сигнал норсализован по мкс-мин, так что информативны только границы)

Э... Фурье не обязан работать на границах. Проверьте сначала на дифференцируемой функции.

[B@R5uk]

MGM, кажется, до меня допёрло, что вы пытались сделать, поправьте, если я не прав. Ключевое слово у вас "интерполяция", то есть вы пытаетесь увеличить набор данных путём добавления новых точек в него на основе имеющихся, но делаете это в спектральной области. Не совсем понятно, правда, для чего: может быть для уменьшения трудоёмкости задачи, когда спектр сигнала/плотности двумерного объекта достаточно компактен? К сожалению, это с преобразованием Фурье это работать не будет (или будет, но не так как вам хочется).

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

[Gleb1964]

Попробуйте использовать метод "Zero padding " , добавить область нулей вокруг массива данных, от которого вы делаете Фурье преобразование. Это работает в двухмерных преобразованиях, должно работать и в одномерном случае. Это sinc-интерполяция в частотном пространстве.
Вбейте Zero padding в поисковик, найдёте подробности метода.