tmtlib
Sep 29 2011, 09:25
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.
анатолий
Sep 29 2011, 10:50
Что-то тут не то. Периодограмма сама считается как усреднение БПФ. Как она может быть проще БПФ?
SPACUM
Sep 29 2011, 11:03
Материал из Википедии — свободной энциклопедии
"Периодограмма — оценка спектральной плотности мощности (СПМ), основанная на вычислении квадрата модуля преобразования Фурье последовательности данных с использованием статистического усреднения:
Периодограмма не является состоятельной оценкой СПМ, поскольку дисперсия такой оценки сравнима с квадратом её математического ожидания. С ростом числа используемых отсчётов значения периодограммы начинают всё быстрее флуктуировать."
И картинка непонятная. А новые методы это интересно. Нужно показать как это раньше измерялось, какая получалась точность и на сколько Ваш метод точнее или быстрее.
tmtlib
Sep 29 2011, 15:38
В моём сообщении проблема с терминологией. Короче, это не периодограмма, а гистограмма.
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
Sep 29 2011, 17:12
Цитата(tmtlib @ Sep 29 2011, 19:38)

Короче, кто что думает по этому поводу? Есть смысл копаться?
Если синусоида идеальная, то любые три точки дадут абсолютно точное значение частоты и это никому не интересно.
Множество людей занималось и занимается определением параметров синусоиды в шуме или среди других сигналов. Возможность обнаружения синусоиды в шуме описывается неравенством Крамера-Рао. И имеется множество методов близких к этому пределу. Сравните Ваш метод с этим пределом например при SNR = 1. Как ни странно у Вас множество конкурентов. Если в поиске набрать "estimation"(оценка), то покажется, что все хотят оценить именно это. Значит пока оценивают недостаточно точно или недостаточно быстро. У Вас есть шанс.
tmtlib
Sep 30 2011, 03:39
Я сделал простую программку, но остаётся нерешенными ряд проблем.
Основная проблема - это как оптимально вытащить период из гистограммы.
Я представляю это так: для грубо определённого перида 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
Sep 30 2011, 06:04
Цитата(tmtlib @ Sep 30 2011, 07:39)

если учесть все максимумы сразу, а не брать один пик. Но как это сделать не доходит.
Да запросто, взять от них БПФ, лучше с окном Гаусса и определить "центр масс". Это разновидность кепстрального метода, который не применяется из-за низкой стойкости к шумам. Попробуйте, может Ваш будет лучше. Хотя бы для SNR = 1;
Kluwert
Oct 3 2011, 10:05
На самом деле, ваш метод, как и многие другие, имеет право на существование. По-сути, он является очередным подоптимальным методом оценки частоты гармоники. Как правильно написал SPACUM, эффективность вашего решения не известна. Что бы понять как оно будет работать, нужно, хотя бы численно (хотя, в вашем случае, вроде не вижу особых проблем посчитать и аналитически) сравнить ваш метод, эффективный метод, полученный методом МП и значением границы Рао-Крамера для оценки частоты гармоники в шуме.
Если окажется, что ваш метод в определенном диапазоне ОСШ близок к оптимуму, то - и прекрасно, будем иметь его в виду. А пока - это не более, чем эвристика.
Цитата(Kluwert @ Oct 3 2011, 14:05)

Если окажется, что ваш метод в определенном диапазоне ОСШ близок к оптимуму, то - и прекрасно, будем иметь его в виду. А пока - это не более, чем эвристика.
У меня вопрос. По технике нет проблем и использую и работает. А в публикациях на графиках есть график 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 @ Oct 3 2011, 19:53)

А какая формула на самом деле? И в какой работе есть обоснование?
Вобщем вся информация о 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 больше. Этот график я и встретил.
Для просмотра полной версии этой страницы, пожалуйста,
пройдите по ссылке.