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.

Что-то тут не то. Периодограмма сама считается как усреднение БПФ. Как она может быть проще БПФ?

Материал из Википедии — свободной энциклопедии

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

И картинка непонятная. А новые методы это интересно. Нужно показать как это раньше измерялось, какая получалась точность и на сколько Ваш метод точнее или быстрее.

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

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 будет говорить о
сдвиге фаз квадратур в ту или иную сторону (вращение вектора по или против часовой стрелке).

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

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

Я сделал простую программку, но остаётся нерешенными ряд проблем.

Основная проблема - это как оптимально вытащить период из гистограммы.
Я представляю это так: для грубо определённого перида T=100 (в битах) имеется некий реальный период T=10001111011110000 (в битах), условно говоря, рассмотрение периодов на гистограмме для каждого увеличения периода в 2 раза даёт 1 бит точности

Гармоника T 100
Гармоника 2T 100[0]
Гармоника 3T ?
Гармоника 4T 100[01]
...
Гармоника 8T 100[011]
... 32T 100[01111] - в квадратных скобках вытащили 5 бит для 32-го максимума гистограммы.
Но так как пики - это некие статистичекские попадания, то там ещё есть некий "центр масс" распределения, т.е. можно проинтерполировать период. Т.е. если максимум для периода 32T находится по оси X на отметке 142 и T=142/32, то подсчитав центр масс T=142.45 мы получим более точное значение T=142.45/32. Возможно можно как-то проще это сделать, если учесть все максимумы сразу, а не брать один пик. Но как это сделать не доходит.

Да запросто, взять от них БПФ, лучше с окном Гаусса и определить "центр масс". Это разновидность кепстрального метода, который не применяется из-за низкой стойкости к шумам. Попробуйте, может Ваш будет лучше. Хотя бы для SNR = 1;

На самом деле, ваш метод, как и многие другие, имеет право на существование. По-сути, он является очередным подоптимальным методом оценки частоты гармоники. Как правильно написал SPACUM, эффективность вашего решения не известна. Что бы понять как оно будет работать, нужно, хотя бы численно (хотя, в вашем случае, вроде не вижу особых проблем посчитать и аналитически) сравнить ваш метод, эффективный метод, полученный методом МП и значением границы Рао-Крамера для оценки частоты гармоники в шуме.

Если окажется, что ваш метод в определенном диапазоне ОСШ близок к оптимуму, то - и прекрасно, будем иметь его в виду. А пока - это не более, чем эвристика.

У меня вопрос. По технике нет проблем и использую и работает. А в публикациях на графиках есть график CRLB. Если подобрать к нему формулу, то средний квадрат частотной ошибки
MSFE ~= 4 / (SNR * sqrt(N * N * N)), где SNR - отношение сигнала к шуму, а N - объем выборки. (формулу я придумал)
Однако даже в этом форуме fontp пишет:
Существует теоретический предел точности измерения частоты по максимуму правдоподобия - предел Крамера-Рао (CRLB).
var(W) = 6/(N*(N-1)*(N-1)*(Es/No)).
В статье "A novel frequency estimator" написано CRLB ~= 12 / (SNR * N * (N * N - 1)).
А на графике точная MSFE названная CRLF.
И в других публикациях формулы разные а графики одинаковые. В последнем случае как и у fontp формулы и графики отличаются более чем в 100 раз.
А какая формула на самом деле? И в какой работе есть обоснование?
(статья: http://www.google.ru/url?sa=t&source=w...eg&cad=rjt)

Вобщем вся информация о CRLB для оценки разных параметров есть в
(http://ebookbrowse.com/eece-522-notes-08-ch-3-crlb-examples-in-book-pdf-d62846321)
для определения чистого синуса в шумах это:
var(f)/Fs^2 = 12 / ((2 * PI)^2 * SNR * N * (N^2 - 1)).
Если извлечь из обеих частей корень, то получим:
RMS(df) / Fs = .551 / ((RMS(сигнал) / RMS(шум)) * sqrt(N * (N^2 - 1))).
А максимальная погрешность раза в 4 больше. Этот график я и встретил.