Зацените метод уточнения максимума спектра FFT по периодограмме. Опции

[tmtlib]

1) Берём FFT без перекрытия, например 1024 точки. Находим максимум F
2) Считаем периодограмму (гистограмму) скользящим окном 512 точек. Получаем гистограмму с максимумами в периодах T, 2*T, 3*T, 4*T,...,
3) По формуле T=1/F находим период Tфурье. Помножив этот период на N=2,3,4,5 получаем грубое попадание в максимумы "гармоник" гистограммы периодов. Корректируем значение максимума по реальному максимуму гистограммы. Например, это была 20-тая гармоника (20*T) с индексом J, тогда мы значем, что период T=J/20, F=1/T=20/J

Вот пример: синим - спектр от FFT (максимум белой линией), зелёным - гармоники периодограммы (белым 30-тая гармоника). Как видите, она "скачет", но реальный максимум можно найти по картинке.

Плюсы: периодограмма считается очень быстро, не требуется скользящего STFT.

Эскизы прикрепленных изображений

 

[tmtlib]

В моём сообщении проблема с терминологией. Короче, это не периодограмма, а гистограмма.

1. По аналогии с FFT, где спектр - это спектр в осях частота-амплитуда,
гистограмма - это спектр в осях период-количество отсчётов.

2. Гистограмма периодов получается следующим образом:
рассмотрим N-тый и M-тый отсчёты исходного сигнала. Для каждой такой пары отсчётов мы можем определить период T=M-N, измеряемый в отсчётах.
Например, для 50-того и 101-ого отсчётов период T=101-50=51.

3. Далее необходимо ввести критерий, по которому будет определяться, прошел ли полный период между отсчётами N и M. В качестве такого критерия можно использовать скалярное произведение между векторами-касательными к функции в точках N и M.
Пусть WAV - это звуковой поток (или I,Q - комплексный звуковой поток)

Для отсчёта N вектор будет иметь координаты (x,y) X=1, Y=WAV[N]-WAV[N-1]
Для отсчёта M вектор будет иметь координаты (x,y) X=1, Y=WAV[M]-WAV[M-1]
При этом возникает некоторая неоднозначность, так вектор всегда направлен "вправо".
Отметим, что при использовании комплексного (квадратурного) сигнала в качестве координат X можно брать непосредственно сам сигнал I:
Отсчёт N: (x,y) X=I[N]-I[N-1], Y=Q[N]-Q[N-1]
Отсчёт M: (x,y) X=I[M]-I[M-1], Y=Q[M]-Q[M-1]
Но так как в общем случае входной сигнал может быть некомплексным, то мы можем искуственно ввести квадратурную составляющую:
I[j]=WAV[j]
Q[j]=WAV[j+1] (по сути это тот самый embedding, что делал Dmitry Terez в своих Pitch Estimation)

Тогда для вещественного сигнала - две касательных до и после отсчёта (в какую сторону выпуклость, чтобы не спутать период с полупериодом и т.п.)
Отсчёт N: (x,y) X=WAV[N]-WAV[N-1], Y=WAV[N]-WAV[N-1]
Отсчёт M: (x,y) X=WAV[M+1]-WAV[M], Y=WAV[M+1]-WAV[M]

Для комплексного - две касательных к окружности
Отсчёт N: (x,y) X=I[N]-I[N-1], Y=Q[N]-Q[N-1]
Отсчёт M: (x,y) X=I[M]-I[M-1], Y=Q[M]-Q[M-1]

Скалярное произведение данных двух векторов - это косинус угла между ними. Критерием прохождения периода между отсчётами M и N является равенство этого угла 0 или близкого к нему (единицы градусов),
что соответствует скалярному произведению V1(x,y) * V2(x,y) = 0.85 ... 1.

4. На практике построение гистограммы периодов выглядит следующим образом:
для всех пар отсчётов N и M считается скалярное произведение касательных векторов. Если это произведение >=0.85, то для
периода T=M-N на гистограмме периодов увеличивается значение на 1. Зелёные линии показывают распределение значений гистограммы по времени.
У метода есть свои недостатки, например он не позволяет найти только основной тон (и то, есть исключения)

5. По указанной выше причине можно сделать так:

1) мы берём "дешёвое" FFT на 256 или 512 точек без перекрытия, определяем интересующий нас максимум F.
2) переводим F в период T (T=1/F)
3) корректируем значение T на основе пика K-того гармоники гистограммы периодов
4) переводим T обратно в частоту.
5) Весь смысл в том, чтобы уловить мелкие изменения частоты: на рисунке выше синяя линяя - это спектрограмма FFT для изменяющегося сигнала,
а зелёные - это спектрограмма от гистограммы периодов. Видно, что изменению частоты FFT на 1 бин соответствует изгибающаяся зелёная линия, более
точно описывающая изменение сигнала по времени.

В методе присутствует скалярное умножение для всех возможных пар отсчётов. Много это или мало?


Также имеется возможность разделять квадратурный сигнал по фазе в гистограмме:
для этого достаточно взять векторное произведение от векторов
(x,y) X=I[N], Y=Q[N] и X=I[N+1], Y=Q[N+1]. Через эти два вектора около комплексного отсчтёта N мы можем провести треугольник,
а векторное произведение - это нормаль Z в трёхмерном пространстве. Соответственно, нормаль Z=+/-1 будет говорить о
сдвиге фаз квадратур в ту или иную сторону (вращение вектора по или против часовой стрелке).

Короче, кто что думает по этому поводу? Есть смысл копаться?

[SPACUM]

Короче, кто что думает по этому поводу? Есть смысл копаться?

Если синусоида идеальная, то любые три точки дадут абсолютно точное значение частоты и это никому не интересно.
Множество людей занималось и занимается определением параметров синусоиды в шуме или среди других сигналов. Возможность обнаружения синусоиды в шуме описывается неравенством Крамера-Рао. И имеется множество методов близких к этому пределу. Сравните Ваш метод с этим пределом например при SNR = 1. Как ни странно у Вас множество конкурентов. Если в поиске набрать "estimation"(оценка), то покажется, что все хотят оценить именно это. Значит пока оценивают недостаточно точно или недостаточно быстро. У Вас есть шанс.