[tmtlib 0]

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.

[aser 0]

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

[SPACUM 0]

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

 

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

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

 

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

 

 

[tmtlib 0]

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

 

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 0]

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

Если синусоида идеальная, то любые три точки дадут абсолютно точное значение частоты и это никому не интересно.

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

[tmtlib 0]

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

 

Основная проблема - это как оптимально вытащить период из гистограммы.

Я представляю это так: для грубо определённого перида 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. Возможно можно как-то проще это сделать, если учесть все максимумы сразу, а не брать один пик. Но как это сделать не доходит.

 

 

[SPACUM 0]

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

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

[Kluwer 0]

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

 

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

[SPACUM 0]

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

У меня вопрос. По технике нет проблем и использую и работает. А в публикациях на графиках есть график 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)

[SPACUM 0]

А какая формула на самом деле? И в какой работе есть обоснование?

Вобщем вся информация о 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 больше. Этот график я и встретил.