![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
Вопрос математикам
Я тут пытаюсь оптимизировать свой алгоритм анализа спектра. И походу опять вылез вопрос с решением одной хитрой системы уравнений. Я её сейчас решаю численно, методом последовательных приближений. Что медленно. Но, возможно есть более прямое решение? Если френды-математики смогут что-то посоветовать хотя бы в какую сторону копать, буду благодарен и напишу благодарность и ссылочку на вас в соответствующем материале на сайте.
Сами мы не местные, в смысле, я же не математик и понимаю, что мой вопрос может быть из разряда просьбы к френду-программисту "а напиши мне к завтрему с нуля аналог Microsoft Word". А, может, наоборот - решение на поверхности, а я его в упор не вижу, так как не математик.
Итак.
Дана вот такая функция Z(m,φ,k) =
где
N - целое, положительное, известно
k - целое от 0 до N/2, известно
m - действительное, превосходящее k не более, чем на 1, неизвестная
φ - действительное (фаза) от -пи до +пи, неизвестная
j - мнимая единица,
Z(m,φ,k) - комплексная функция.
Цель - найти m,φ
Известно отношение модулей |Z(m,φ,k)|/|Z(m,φ,k+1)|, известны комплексные фазы (аргументы) Z(m,φ,k) и Z(m,φ,k+1).
Сейчас я численно решаю систему
|Z(m,φ,k)|/|Z(m,φ,k+1)| = A
Arg(Z(m,φ,k+1)) = B
- путем попеременного варьирования m,φ
Физический смысл задачи - найти частоту и фазу гармоники, которая при разложении в дискретный спектр Фурье "размазалась" на множество гармоник по той причине, что не попала точно в одну из фурье-гармоник.
Музыкальный смысл задачи - уточнить высоту ноты, которую играют или поют в данное мгновение.
Иллюстрация: нота ля-бемоль в музыкальной записи с частотой дискретизации 44100 Гц. В исходном сигнале - единственная гармоника в районе 415 Гц и амплитудой 10000, а в результате преобразования Фурье (N = 1000) получается ЭТО:
Есть идеи?
Сами мы не местные, в смысле, я же не математик и понимаю, что мой вопрос может быть из разряда просьбы к френду-программисту "а напиши мне к завтрему с нуля аналог Microsoft Word". А, может, наоборот - решение на поверхности, а я его в упор не вижу, так как не математик.
Итак.
Дана вот такая функция Z(m,φ,k) =
где
N - целое, положительное, известно
k - целое от 0 до N/2, известно
m - действительное, превосходящее k не более, чем на 1, неизвестная
φ - действительное (фаза) от -пи до +пи, неизвестная
j - мнимая единица,
Z(m,φ,k) - комплексная функция.
Цель - найти m,φ
Известно отношение модулей |Z(m,φ,k)|/|Z(m,φ,k+1)|, известны комплексные фазы (аргументы) Z(m,φ,k) и Z(m,φ,k+1).
Сейчас я численно решаю систему
|Z(m,φ,k)|/|Z(m,φ,k+1)| = A
Arg(Z(m,φ,k+1)) = B
- путем попеременного варьирования m,φ
Физический смысл задачи - найти частоту и фазу гармоники, которая при разложении в дискретный спектр Фурье "размазалась" на множество гармоник по той причине, что не попала точно в одну из фурье-гармоник.
Музыкальный смысл задачи - уточнить высоту ноты, которую играют или поют в данное мгновение.
Иллюстрация: нота ля-бемоль в музыкальной записи с частотой дискретизации 44100 Гц. В исходном сигнале - единственная гармоника в районе 415 Гц и амплитудой 10000, а в результате преобразования Фурье (N = 1000) получается ЭТО:
Есть идеи?
no subject
тут названия алгоритмов в frequency domain
no subject
no subject
Сходу приходят в голову всякие идеи вроде найти первый и второй максимумы на fft и между ними искать частоту максимума или половинным делением (но тогда надо придумать функцию которая монотонно зависит от положения максимума между бинами) или попытаться определить ширину и количество лепестков (sin x/x) а из них - таблично искать положение максимума. Но это все надо экспериментировать.
no subject
насчет всех не скажу, а обычный fft находит с точностью до (sample rate/N)
ну я примерно так и делаю, как ты говоришь - половинным делением решаю систему.
но это же n итераций
а вдруг, думаю, эта система имеет решение, выражаемое в виде формулы, а я его просто не вижу, т.к. непрофессиональный глаз не заточен видеть такие вещи
no subject
берут ДВА fft для соседних участков
потом смотрят, как изменилась фаза одной гармоники
и на основании изменения фазы сразу получают уточненную частоту
но это надо именно два последовательных fft иметь
no subject
no subject
и задержку увеличивать-уменьшать, пока не найдешь максимум
то есть i штук итераций над функцией со сложностью N, чтобы найти ОДНУ частоту из тех, что есть в сигнале
для m частот понадобится сложность i*N*m
а ПФ за NlogN находит сразу все m
но "размазанные"
no subject
Там, конечно, есть нюансы и трюки, так как БПФ описывает только периодические сигналы. Например, можно после куска сигнала столько же нулей дописать, чтобы конец не наезжал на начало.
no subject
В результате мы получим N значений функции для разных задержек (периодов, частот). А если у нас искомая задержка лежит в промежутке между вычисленными значениями, то возвращаемся к той же проблеме - так?
no subject
no subject
Обязательно буду пробовать вариант 1.
Вариант 2 ясен, но для высокой точности требует вычислений огромных БПФ, что по затратам времени будет шило на мыло.
Вариант 3 - ну, собственно, приведенный выше мой алгоритм дает идеальную точность при одной гармонике в исходном сигнале, и при небольшом числе далеко разнесённых тоже всё OK. Так что я теперь заморачиваюсь вопросом производительности.
Вариант 4 - тоже обязательно попробую! И оформлю на сайте, сославшись на вас как на автора идеи - так что, если хотите увидеть там реальное ФИО вместо taras_, кидайте сюда :)
no subject
- а вот это вот не понял, почему пренебречь?
допустим k поcредине в районе N/4, m - там же, тогда exp(-2pi*j*(m+k)/N) ~ exp(-pi*j/2) ~ -j ... э... ну и?
no subject
В знаменателе второго слагаемого будет -j-1, то есть его модуль sqrt(2). Допустим, m-k = 0.5, N = 1000, тогда в модуль знаменателя первого слагаемого
|exp(2pi*j/2000) - 1| = 3.1416e-03
то есть примерно в 450 раз меньше.
Теперь посмотрим на числители. Поскольку exp(2*pi*j*k) = 1 для всех целых k, то
exp(2*pi*j*(m-k)) = exp(2*pi*j*(m+k)) = exp(2*pi*j*m)
из чего видно, что числители являются комплексно сопряженными. Следовательно, они имеют одинаковый модуль. Т.о. первое слагаемое по модулю примерно в 450 раз больше второго.
Если ссылаться, то можно, например, сюда: http://www.tvd-home.ru/
no subject
ещё раз спасибо :)
no subject
Насколько я тогда для себя понял (на глубокую научность я не претендую), все эти широкие крылья в спектре обусловлены тем, что синусоида как бы обрывается на краях диапазона. Поскольку DFT неявно предполагает, что диапазон дискретизации как бы повторяется периодически, то если период сигнала не укладывается целое число раз в этот диапазон, то отсюда и лезут эти разрывы. Фурье преобразование любые разрывы очень не любит (даже в производных, причём чем ниже степень производной, тем сильнее реакция) и отвечает на их добавлением вот этих крыльев. Фактически, это свёртка истинного спектра сигнала (который практически дельта функция) со спектром прямоугольного окна.
С этим принято бороться (в английской вики -- Window_function) перемножением сигнала во временной области с более хорошим окном. Причём в зависимости от того, что именно нужно извлечь из сигнала используют разные окна. Идея в том, чтобы так обнулить функцию на краях диапазона, чтобы она оставалась там гладкой. Потерю информации в таком случае можно восполнить, беря соседние диапазоны с нахлёстом друг на друга.
Я в свое время использовал традиционное окно Ганна (косинус с подставкой). Это фактически эквивалентно двукратному дискретному дифференцированию спектра в частотной области. Крылья дифференцированием практически обнуляются, остаются только два значения в центре, причём сумма их амплитуд как раз хорошо совпадала с амплитудой синусоиды. Мне этого было достаточно. Фазу (разность фаз) я из теоремы косинусов тогда искал, суммируя сигнал с другим сигналом, но получалось гораздо хуже.
В вашем случае, может быть имеет смысл взять узкий Гаусс (так чтобы на краях диапазона он был бы экспоненциально мал). Фурье преобразование от Гаусса тоже будет Гауссом (причём все его нелинейные параметры аналитически рассчитываются), так что подгонкой спектра можно было бы точно найти его середину.
Не знаю, актуально ли это сейчас для вас, просто мне было приятно тряхнуть стариной :)
no subject