Вопрос математикам

[psilogic]

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

Сами мы не местные, в смысле, я же не математик и понимаю, что мой вопрос может быть из разряда просьбы к френду-программисту "а напиши мне к завтрему с нуля аналог 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) получается ЭТО:



Есть идеи?

Comments

[metaclass.livejournal.com]

https://en.wikipedia.org/wiki/Pitch_detection_algorithm
тут названия алгоритмов в frequency domain

[kray-zemli.livejournal.com]

Преобразование Фурье -- лженаука. Смотри автокорреляционную функцию.

[http://users.livejournal.com/taras_/]

1. Если не нужна какая-то супер точность, то можно использовать простенькую интерполяцию: http://dspguru.com/dsp/howtos/how-to-interpolate-fft-peak 2. Идеальная интерполяция ДПФ получится, если увеличивать длину сигнала, дополняя его нулями (http://www.dsprelated.com/freebooks/sasp/Spectral_Interpolation.html). Но придется делать fft для более длинных сигналов. 3. Если нужна именно супер точность, то соответствующие методы есть, но я в них не силен. 4. Глядя на функцию Z можно заметить, что (m-k)/N - величина маленькая, значит exp(2pi*j*(m-k)/N) примерно равна 1. То есть знаменатель первого слагаемого очень мал. Если мы не находимся на краю диапазона, когда m+k близко к 0 или N, то вторым слагаемым можно смело пренебречь, после чего задачка решается аналитически.

[pofig.livejournal.com]

Я не математик, но мне тоже когда-то приходилось искать амплитуду и фазу дискретизованного синусоидального сигнала.

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

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

Я в свое время использовал традиционное окно Ганна (косинус с подставкой). Это фактически эквивалентно двукратному дискретному дифференцированию спектра в частотной области. Крылья дифференцированием практически обнуляются, остаются только два значения в центре, причём сумма их амплитуд как раз хорошо совпадала с амплитудой синусоиды. Мне этого было достаточно. Фазу (разность фаз) я из теоремы косинусов тогда искал, суммируя сигнал с другим сигналом, но получалось гораздо хуже.

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

Не знаю, актуально ли это сейчас для вас, просто мне было приятно тряхнуть стариной :)