psilogic: (wpriz)
[personal profile] 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) получается ЭТО:



Есть идеи?

Date: 2015-12-26 01:16 pm (UTC)
From: [identity profile] metaclass.livejournal.com
https://en.wikipedia.org/wiki/Pitch_detection_algorithm
тут названия алгоритмов в frequency domain

Date: 2015-12-26 01:19 pm (UTC)
From: [identity profile] psilogic.livejournal.com
Ага, я в курсе, даже нашел один интересный метод, но к этой задачке он прямого отношения не имеет.

Date: 2015-12-26 01:44 pm (UTC)
From: [identity profile] metaclass.livejournal.com
Все pitch-tracking алгоритмы с точностью до fft bin что-ли определяют?

Сходу приходят в голову всякие идеи вроде найти первый и второй максимумы на fft и между ними искать частоту максимума или половинным делением (но тогда надо придумать функцию которая монотонно зависит от положения максимума между бинами) или попытаться определить ширину и количество лепестков (sin x/x) а из них - таблично искать положение максимума. Но это все надо экспериментировать.

Date: 2015-12-26 02:59 pm (UTC)
From: [identity profile] psilogic.livejournal.com
[ Все pitch-tracking алгоритмы с точностью до fft bin что-ли определяют? ]

насчет всех не скажу, а обычный fft находит с точностью до (sample rate/N)

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

Date: 2015-12-26 03:01 pm (UTC)
From: [identity profile] psilogic.livejournal.com
есть еще такой алгоритм
берут ДВА fft для соседних участков
потом смотрят, как изменилась фаза одной гармоники
и на основании изменения фазы сразу получают уточненную частоту

но это надо именно два последовательных fft иметь

Date: 2015-12-26 02:21 pm (UTC)
From: [identity profile] kray-zemli.livejournal.com
Преобразование Фурье -- лженаука. Смотри автокорреляционную функцию.

Date: 2015-12-26 02:51 pm (UTC)
From: [identity profile] psilogic.livejournal.com
насколько я понимаю, надо будет смотреть корреляцию сигнала с самим собой с задержкой
и задержку увеличивать-уменьшать, пока не найдешь максимум
то есть i штук итераций над функцией со сложностью N, чтобы найти ОДНУ частоту из тех, что есть в сигнале
для m частот понадобится сложность i*N*m

а ПФ за NlogN находит сразу все m
но "размазанные"

Date: 2015-12-26 03:10 pm (UTC)
From: [identity profile] kray-zemli.livejournal.com
Автокорреляционная функция вычисляется через БПФ за NlogN. Нужно все коэффициенты взять по модулю и возвести в квадрат, потом сделать обратное преобразование. Она самая и получится.

Там, конечно, есть нюансы и трюки, так как БПФ описывает только периодические сигналы. Например, можно после куска сигнала столько же нулей дописать, чтобы конец не наезжал на начало.

Date: 2015-12-26 03:21 pm (UTC)
From: [identity profile] psilogic.livejournal.com
я так понял, что ты вот об этом?

В результате мы получим N значений функции для разных задержек (периодов, частот). А если у нас искомая задержка лежит в промежутке между вычисленными значениями, то возвращаемся к той же проблеме - так?

Date: 2015-12-26 04:41 pm (UTC)
From: [identity profile] 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, то вторым слагаемым можно смело пренебречь, после чего задачка решается аналитически.

Date: 2015-12-26 05:15 pm (UTC)
From: [identity profile] psilogic.livejournal.com
Спасибо вам огромное, то что нужно!!!

Обязательно буду пробовать вариант 1.

Вариант 2 ясен, но для высокой точности требует вычислений огромных БПФ, что по затратам времени будет шило на мыло.

Вариант 3 - ну, собственно, приведенный выше мой алгоритм дает идеальную точность при одной гармонике в исходном сигнале, и при небольшом числе далеко разнесённых тоже всё OK. Так что я теперь заморачиваюсь вопросом производительности.

Вариант 4 - тоже обязательно попробую! И оформлю на сайте, сославшись на вас как на автора идеи - так что, если хотите увидеть там реальное ФИО вместо taras_, кидайте сюда :)

Date: 2015-12-26 05:22 pm (UTC)
From: [identity profile] psilogic.livejournal.com
[ то вторым слагаемым можно смело пренебречь ]

- а вот это вот не понял, почему пренебречь?

допустим k поcредине в районе N/4, m - там же, тогда exp(-2pi*j*(m+k)/N) ~ exp(-pi*j/2) ~ -j ... э... ну и?

Date: 2015-12-26 06:12 pm (UTC)
From: [identity profile] http://users.livejournal.com/taras_/
[допустим k поcредине в районе N/4, m - там же, тогда exp(-2pi*j*(m+k)/N) ~ exp(-pi*j/2) ~ -j ... э... ну и?]

В знаменателе второго слагаемого будет -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/

Date: 2015-12-26 06:49 pm (UTC)
From: [identity profile] psilogic.livejournal.com
а, блин, тормозил, не туда смотрел
ещё раз спасибо :)

Date: 2016-01-02 09:54 am (UTC)
From: [identity profile] pofig.livejournal.com
Я не математик, но мне тоже когда-то приходилось искать амплитуду и фазу дискретизованного синусоидального сигнала.

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

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

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

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

Не знаю, актуально ли это сейчас для вас, просто мне было приятно тряхнуть стариной :)
Edited Date: 2016-01-03 12:11 pm (UTC)

Date: 2016-01-03 06:58 pm (UTC)
From: [identity profile] psilogic.livejournal.com
Да, оконные функции эти размазанные пики сужают, так что смысл имеют, да.
Page generated Sep. 21st, 2017 03:26 am
Powered by Dreamwidth Studios