Имеем на выходе датчика синусоиду примерно 2кГц в течении 0.5 сек. , т.е. цуг из примерно 1000 периодов. В зависимости от физ. величины частота меняется на +-10% Считается, что во время цуга частота не меняется. Надо определить частоту с точностью хотя бы до 0.01Гц. На данный момент метод используетсяпрямое измерение периода каждого колебания с дальнейшей статистической обработкой. При существующем уровне помех, наводок и т.д. точность получается не лучше 0.1Гц. Вопрос: как методами ЦОС уточнить результат измерения? Это вообще возможно теоретически? Сразу должен оговориться, что я о цифровой обработке имею весьма общие представления и никакой практики. Мне в голову приходит только свертка с синусоидами от (f0 – 0.1Гц) до ( f0 + 0.1Гц) с шагом 0.01Гц (f0-приблицительно измеренная частота) и поиском "резонанса". Может быть можно как по другому, попроще?

У Вас есть все необходимое для такого измерения (Tизмх примерно 0.5 сек и вполне приемлемая частота Fx входного сигнала). Далее, предположим, у Вас есть опорная, точно известная частота Fo=1/To.
Если Tизмх=Тх*Nx ( Nx- отсчитываемое счетчиком1 периоды измеряемой частоты), Tизмo=To*nxo ( nxo- отсчитываемое счетчиком2 периоды опорной частоты за время прохождения Nx периодов) и, обеспечив условие Tизмх=Tизмo, получаем после несложного преобразования Fx=Fo*Nx/nxo.

Существует теоретический предел точности измерения частоты по максимуму правдоподобия - предел Крамера-Рао (CRLB).
var(W) = 6/(N*(N-1)*(N-1)*(Es/No))
W=2*pi*f - частота, pi- это пи, N- длина блока данных, Es/No-отношение сигнал/шум
В принципе, если сигнал/шум стремится к бесконечности точность может быть как угодно хорошей.
Существуют и алгоритмы приближающие эту оценку.
В частности, могут измеряться отсчёты DFT вблизи области интереса, находится максимум спектра и по трём точкам вблизи максимума строится парабола - аналитическое положение вершины параболы принимается за оценку частоты
Здесь есть матлабовские модели
http://home.comcast.net/~kootsoop/EricJ2/index.htm
и здесь
http://home.comcast.net/~kootsoop/freqalgs.htm

Судя по тому, что вы производите статистическую обработку каждого периода, они у вас разные получаются. Значит сигнал-то не чистая синусойда, а имеет фазовые шумы. Боюсь что 1000 периодов для требуемой точности будет недостаточно.

А каковы требования к длине выборки при оценке частоты сигнала путём оценки максимума спектра?
Как я понимаю чем выборка больше тем лучше?

Но вот скажем сигнал: синус частотой 40...60Гц. Требуемая точность =>0,01Гц.
Частота дискретизации 6400 (для примера примеру).

для fontp: какова должна быть длина выборки, что бы Macleods estimtor дал результат с требуемой точностью?

Чем больше выборка, естественно тем лучше. Спектральное разрешение всегда на сложных спектрах у нас
определяется критерием Релея df*T порядка 1. Но на простых сингулярных спектрах - когда спектральный пик 1, или даже несколько, но расположеных далеко, точность измерения частоты может быть значительно выше и определяется приведеным выше статистическим критерием Крамера-Рао. Он определяет теоретический предел для измерения параметров в присутствии шума. Для такого параметра как частота я его привёл, в нормированых единицах. Все стремятся добиться чего-то похожего на этот критерий.

При возрастании отношения сигнал/шум к бесконечности предельная точность может быть как угодно высокой при любой длине блока. Это по теории. А практические методы могут отличаться в трех отношениях

1. При очень низком SNR не работают никакие методы. Поэтому можно говорить о пороговом снизу значении SNR когда отказывает метод. Понятно, что для каждого метода этот параметр разный.
2. Энергетические потери в рабочем диапазоне. Хорошие методы в этой промежуточной области SNR все дают меньше 1 дб потерь. Или наоборот приличный метод по порядку с точностью до десятков процентов ложится на CRLB
Иначе их не публикуют в научных журналах.
3. Систематические ошибки при достаточно высоких SNR. При высоких SNR кривая CRLB уходит сама по себе, а точность реального метода становится постоянной и более от SNR не зависит.
Если не делать никакой интерполяции - то эта систематическая ошибка плохая и равна примерно половине размера бина преобразования. Если не добавлять нулей, то метод Маклеода - лучший среди известных. Умелым добавлением нулей можно добиться и более хороших результатов. Но это затратно вычислительно

Давайте посчитаем пример

Стандартное отклонение df/Fs = sqrt(6/((2*pi)**2*N*N*N*SNR))
0.01/6400 = примерно sqrt(1/6.6*N*N*N*SNR)

Если SNR=100
N**1.5 = 25000, N ~ 800

Спасибо.
Да, ссылки на оценщики больше не работают. Может выложите статью из журнала...?

так она выложена в другой, новой теме, недавно, вместе с другими
http://electronix.ru/forum/index.php?showt...mp;#entry695667

Ссылку на исчерпывающий систематический подход к квадратичной интерполяции спектра с дополнением сигнала нулями я тоже давал
https://ccrma.stanford.edu/STANM/stanm/node3.html
Документы от STANM-114 до STANM-118
Там есть всё - от использования окон и борьбы со смещением до влияния нестабильности параметров и модуляции
"... Это все, что я могу сказать про войну во Вьетнаме"©

А как насчёт способа когда определяем две опорные точки (максимумы или переходы через ноль), далее интерполируем в окрестностях этих точек (что бы уточнить положение) и по ним определяем частоту.
Этот способ может дать выигрыш в скорости? (по сравнению с алгоритмом по максимуму БПФ)

Это всё может работать только когда практически отсутствует шум. Это же дифференцирование сигнала.
Причём двойное - один раз сигнала для определения точки пересечения с нулём, а второй раз разность во времени делта-Т для определения периода.
Здесь на форуме приводились же работы китайцев по численному дифференцированию сетевого напряжения.
У них там очень хорошие результаты, когда шума нет.
В присутствии шума там приведена какая-то ерунда - чем больше шум, тем точнее (в одиночных измерениях).
Какая-то недостоверная статистически китайская ерунда.
С шумом китайцы не разобрались, про шум им, видимо, не читали, бакалаврам
При средних и сильных шумах обычно используется ДПФ

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

Можете пробовать разные методы, для разных задач какие-нибудь оригинальные могут оказаться почему-то предпочтительней.
Но всегда нужно стремиться в качестве ориентира к CRLB, поскольку эта предельная оценка максимального правдоподобия от метода не зависит. Хотя, конечно, если шума нет, то этот критерий не играет и всё определяется только систематической погрешностью метода - тогда ДПФ лучше не использовать.

Если частота меняется на +/- 10%, то разбейте ваш сигнал на куски по 20 периодов. Далее, для каждого такого куска сосчитайте скалярное произведение первой половины на вторую его половину. Возьмите аргумент получившегося комплексного числа и поделите его на количество сэмплов в половине куска. Вы получите набег фазы на сэмпл. Сосчитайте эту величину для каждого куска и усредните. Поделите на 2*pi и умножьте на частоту дискретизации. Вы получите значение частоты.

Это тоже приемлемый метод. Delay-Multiply-Add (DMA). Но он при заданной длине блока в области низких и средних уровней шума даёт немного меньшую (но соизмеримую) точность чем интерполяция спектра Фурье. При высоком отношении SNR может быть даже предпочтительней по точности, хотя никогда не выходит на предельный уровень CRLB. Вычислительно проще спектральных оценок, но требует вычисления арктангенса.
Подробно описан например в атачменте применительно к модемостроению

Существует модификация метода, позволяющая измерять таким способом частоту в любом широком диапазоне.
"Энергетические потери" по отношению к CRLB не превышают 0.5 дб при SNR>1, правда они на последнем этапе дополнили метод таки квадратичной интерполяцией спектра вблизи получившегося значения частоты (это так называемое ML-extension можно навесить в принципе на любой измеритель частоты).
Я приводил когда-то статью по этому Delay-Multiply-Add-Rotate-Add-Decimate(DMA-RAD). Смысл в том, что чем больше сдвиг в вычисляемых корреляциях, тем выше точность, но меньше диапазон измерения. Вот и проводятся вычисления итеративно начиная от малых задержек в корреляции к большим, убирая сначала большие смещения частоты, измеренные с меньшей точностью и последовательно уточняя оценку в меньшем диапазоне

http://electronix.ru/forum/index.php?act=A...st&id=20411

Я так понимаю, что это всё рекомендации для общего случая.

А у меня сигнал известен - синусоида сетевой частоты (42.5... 57.5Гц). Синусоида может содержать гаромники кратные 5Гц. Амплитуда гармоник не превышает 10% от основной.
Погрешность вычисления частоты должна быть >0,01Гц.
Алгоритм тестируется с гауссовым шумом при ОСШ >40дБ.

А если отфильтровать сигнал полосовиком! Диапазон частоты ведь известен. Скажем полосовой БИХ-фильтр на частоту 40-70 Гц. Тогда же можно будет мерять частоту по двум опорным точкам (если нтерполировать вблизи них). Да! Частота должна ещё усредняться по нескольким десяткам периодов.

P.S. Просто здаётся мне: что БПФ, что автокореляция - это много процессорного времени.

Вариант такой есть. Умножить на комплексную экспоненту и взять фазу. Угол наклона фазы покажет разность частоты синусоиды и частоты комплексной экспоненты. соответсвенно df = dPhi/dT. Таким образом увеличивая интервал обработки вы получите сколь угодно высокую точность. Такую обработка лучше чем просто бпф.

Увеличивая интервал обработки можно получить как-угодную точность каким-угодно методом. В этих задачах
оптимальность означает - как быстро растёт точность при увеличении интервала. Оптимальность это не только научные понты для ученых. Мы не можем реально значительно увеличивать интервал измерения, только в пределах стабильности самих параметров сигнала (амплитуды, частоты)

Такая обработка ничем не лучше спектрального оценивания, хотя по определению должна давать точность по максимуму правдоподобия (CRLB), если честно подогнать минимум квадратов. Только вычислительно это даже сложнее, ибо Вы должны
1. Умножить на центральную комплексную экспоненту (снести сигнал в 0) и отфильтровать ФНЧ.
2. В каждой точке (возможно после какого то усреднения-сабсамплинга) взять арктангенс с поправкой на непрерывность фазы (2*pi туда-сюда)
3. Подогнать линейную регрессию Ф=at+B, a - частота
Вычислительно это будет скорее сложнее чем спектральное оценивание, а не проще, даже с применением CORDIC для арктангенса. (Не говоря уже о том, что это будет работать только при достаточно высоком SNR)

Сравните с оцениванием спектра и обратной квадратичной интерполяцией. Если интервал частот узкий, то сразу можно начинать с того, что называют ML-extension

1. Умножить на центральную комплексную экспоненту (снести сигнал в 0) и отфильтровать ФНЧ. Соответственно прорежение.
2. Взять 5-7 сумм типа ДПФ в интервале [-7.5, 7.5]. Оценка энергии
3. Найти максимум энергии и по трём точкам построить параболу. Аргумент максимума - частота

Точность будет примерно такая-же.

Хмм.. много или не много?.. Complex 16-bit FFT:

Интересно, почему тогда системы слежения выделяют фазу и производят подстройку частоты именно на основе фазового измерения, а не на основе БПФ?
На самом деле такая оценка дает точность соответсвующую методу максимального правдоподобия. Именно поэтому она лучше БПФ. При бпф после умножения на комплексную экспоненту все значения складываются (этакий фильтр) в результате вы получите некую среднюю фазу и не сможете понять как она меняется на интервале обработки. Здесь же у вас есть информация о фазе. Зная что у нас одна синусоида это линейная функция. Это можно использовать для определения df. Кстати в этом случае раскрыть периодичность арктангенса вполне не сложно читай например здесь

ЗЫ кстати данная реализация куда быстее бпф, так как арктангенс можно вообще не вычислять а задать таблично.

Она не лучше спектрального оценивания ("БПФ"), именно потому, что в очень широких пределах SNR (0-40дб) после интерполяции спектра это "БПФ" даёт именно точность максимального правдоподобия, соответствующую предельной оценке Крамера-Рао. Ваше утверждение про "лучше" - голословно, а я приводил ссылки на работы, где предельные оценки достигнуты. Спектральные оценки ДПФ это то же самое синхронное детектирование и линейную фазу гармоники ДПФ фильтруют оптимально, чего уж там "понимать". Нельзя быть святее Папы.

Возможно "такая оценка" с измерением фазы и будет лучше, но только при ещё более высоких SNR.
При SNR<10 дб "такая оценка" наоборот работать будет плохо, поскольку сама возможность восстанавливать линейную фазу по непрерывности по значению в пределах 2*pi исчезает - всё забито шумом. Это я говорю по собственному опыту. Я все эти способы пробовал.

По производительности. Это арктангенс двух переменных с восстановлением непрерывности фазы
Ph[m] = arctg2(SIm,SQm);
if (ABS(Ph[m-1]-Ph[m]) > M_PI) {// phase jump is detected
if (ABS(Ph[m-1]-Ph[m]-2*M_PI) <= M_PI) {
Ph[m] += 2*M_PI;
}else {
Ph[m] -= 2*M_PI;
}
}

1. 2 сравнения чтобы выяснить квадрант
2. Деление (x/y) или (y/x) после сравнения кто из них меньше, чтобы не делить на примерно 0 и не получать примерно бесконечность
3. Чтобы считать арктангенс уже от одной переменной в лоб по таблице с приличной точностью таблица должна быть уж очень большой. Поэтому обычно считают по быстрому алгоритму CORDIC по небольшой "фрактальной" таблице для одного квадранта.
4. Восстановление непрерывности фазы - еще несколько сравнений и вычитаний как написано выше

Всё про всё займет если не аппаратно тактов 50-100. Итого имеем пусть 50*N - N длина блока. Вычисление FFT требует примерно 4*N*log(N), DFT на современной архитектуре с двумя МАС-аккумуляторами - просто k*N, k - число вычисляемых гармоник в диапазоне. Плюс корреляция (Ph[i]*i), необходимая для подгонки линейной модели (регрессии, сумма(i*i) пусть расчитано предварительно, частота по линейной регрессии = сумма(Ph[i]*i)/сумма(i*i)). Арктангенсы выгодней считать получается только ассимптотически при очень больших блоках N и то только в некоторых случаях, когда 50 гармоник DFT (синхронных детекторов) недостаточно.

ЗЫ. Что касается "систем слежения" (видимо ФАПЧ c обратной связью) то можно совершенно определённо сказать, что всякий ФАПЧ в синхронизм входит медленней, чем это делают измерительные feed-forward системы, причем последние для области низких SNR работают как раз чаще всего на оценивании спектра через интерполяцию.
Если не ставить целью провести измерение частоты максимально быстро, то в принципе, задачу можно решать стандартно именно с помощью ФАПЧ. Вполне приемлемое решение, но оно никогда не будет оптимально по времени измерения. Есть даже такие готовые микросхемы квадратурной демодуляции с фазовым детектором и цифровым генератором. Код управляющий частотой и соответствует частоте входного сигнала. Только это всегда будет немного тормозное решение, уступающее схемам по прямому измерению без обратной связи по точности/времени измерения

Без всякого деления и арктангенса с одним сравнением сразу кордиком посчитаем аргумент комплексного числа.



Прежде чем "арктангенсы" считать ещё проредим простейшим скользящим средним CIC-ом раз в 10.

Кордик в любом случае даст десятки тактов. Скорее 50, чем 10, поскольку там тоже всё на сравнениях. А для DSP сравнения значительно затратнее МАС. Я просто пример привёл с СORDIC-ом на один квадрант на вскидку.
CIC может решать. Раз в десять только, если диапазон измеряемых частот в десять раз меньше частоты Найквиста.
Там даже не CIC, а простое усреднение нужно до вычисления арктангенса, если блочный алгоритм. Усреднять возможно только пока не подавляется синком полезный сигнал. Усреднение ещё и полезно, чтобы не позволить шуму делать сбои при восстановлении фазы по непрерывности, увеличивает SNR. Усреднение принято ))
Но похожую предварительную обработку можно проделать и для DFT, это не аргумент при сравнении.

Хорошо. Такой алгоритм будет оптимален по точности и вычислительно. Но не такой робастный как спектральное оценивание и грохнется раньше при снижении SNR. Из-за необходимости устранять скачки фазы

Но именно так в модемах и делают, работает при плохом сигнал/шум, грохнуться ведь в принципе любой метод может, в общем прям так сразу не очевидно что спектральное оценивание робастней.

В BURST-модемах делают все три метода (спектр, корреляция и линейная фаза). В общем-то все три уязвимы.
Но особенно тот, который строит модель линейной фазы по циклической.
Шумовой импульс по фазе в любой точке с амплитудой порядка ПИ погубит модель.
Вам известен простой алгоритм с интегральными свойствами? Чтобы строил линейную модель игнорируя одиночные выбросы. Так сказать алгоритм циклической регрессии, но чтобы без перебора. Мне нет.
(Хотя и можно, наверное, что то такое приладить как вычисление разности фаз между соседними точками с последующим ранговым усреднением медианой. Но это уже будет много более затратно)

Восстановление фазы по непрерывности я пробовал и мне сравнительно не понравилось. Очень быстро и резко кончается.
В спектре Фурье одиночные импульсы размазываются и никому не мешают. Хотя, конечно, наоборот гармониченская помеха в рабочем диапазоне критична.

Думаю, что поставленные требования к точности оценки частоты подразумевают, что отношение сигнал-шум велико (гораздо больше 10 дБ). Именно поэтому я и предложил такой вариант.



Как раз в этом случае высокая точность арктангенса и не нужна, поскольку мы измеряем не мгновенную фазу, а "медленное" изменение этой фазы. Так например если есть таблица из 360-и значений арктангенса, то средняя ошибка измерения одной фазы будет 0.5 градуса. 0.5 градуса на интервале T = 1 сек приведут к ошибке по частоте равной 0.5*pi/180 =0.008 Гц.

Для высокого отношения сигнал/шум метод может оказаться хорошим



1 сек это не по максимуму правдоподобия.
По максимуму правдоподобия это 1/8 сек (у TigerSHARCa), согласно критерию при 20 дб, а при большем отношении SNR ещё меньше.
У Вас получится примерно в среднем sqrt(2)*0.008*8 (sqrt(2) поскольку разность фаз на интервале по Вашей логике по двум точкам в конце и начале интервала), т.е. за 1/8 секунды по двум арктангенсам Вы так получите почти 0.1 гц ошибки. Мы ведь говорим о предельных оценках. Теория говорит, что
var(F) = 3*(Fs*Fs)/(2*pi*pi*N*(N-1)*(N-1)*(Es/No))
При этом нужная таблица получится в 10 раз больше для точности 0.01 гц если по двум точкам считать.
(По двум точкам на границах всего интервала считать нельзя, поскольку набег фазы должен быть меньше pi,
но для оценок, ладно сойдет. Мы будем на самом деле считать разность фаз между соседними арктангенсами, а потом усредним по всем разностям)

На самом деле строя линейную регрессию по минимуму квадратов по многим измерениям фазы точность повышается примерно как корень из числа измерений и по отношению к систематическим ошибкам в таблице. Здесь есть ещё небольшой резерв в несколько раз. В любом случае таблица нужна в несколько раз больше. И всё это ещё нужно очень хорошо сбалансировать, по равенству (<=) систематической ошибки примерно случайной, вызваной физическим шумом - посредством баланса по сколько отсчетов усреднять и сколько арктангенсов брать для регрессии=(полное число отсчетов)/(размер лаптя усреднения до вычисления фазы).

Всем здравствуйте!
Тема интересная, позвольте выложить матлабовский скрипт, в котором я набросал алгоритм оценки частоты комплексной экспоненты на основе планарной фильтрации.
Ничего сверхестественного, просто один из adhoc способов. Идея в том, чтобы умножая отсчеты комплексной экспоненты на свою задержанную и взятую с комплексным сопряжением копию, получить последовательность фазоров, угол которых пропорционален частоте исходной экспоненты. Затем нужно сложить полученные фазоры как комплексные числа, получив на комплексной плоскости длиный вектор, угол которого опять же пропорционален искомой частоте(эта операция и есть планарная фильтрация). Таким образом уменьшается воздействие шума на оценку. Далее вычисляем угол этого длинного вектора, умножаем на коэффициент Fs/(2*pi) и получаем оценку частоты.
Правда полученная таким образом оценка достаточно сильно подвержена влиянию шума.
Чтобы несколько улучшить оценку, я пошел немного дальше и используя вычисленную оценку частоты, создаю комплексный КИХ фильтр с ИХ в виде 1 периода комплексной экспоненты на полученной оценкой частоте. Получается такой комплексный ППФ со средней частотой совпадающей с нашей оценкой частоты.
Затем пропускаю входной сигнал через этот ППФ и снова нахожу оценку по вышеописанному методу.
Засчет второй итерации после фильтрации, погрешность оценки может уменьшиться на порядок по сравнению с первой оценкой.
Вот такой вот способ.
Правда, вторая итерация имеет смысл только в определенном диапазоне отношения сигнал шум (если ОСШ выше 45 дБ, то полезный эффект от второй итерации уже незначителен)
Планарная фильтрация используется в ряде прямых(feedforward) методов синхронизации, описанных например в
в часто упоминающейся на форуме книге Synchronization Techniques for Digital Receivers.

Собственно, скрипт:

%test_carr_freq_ff_est
script
clear
clc
close all

Fs=6400;
Tau=1.00;%
Nlen=Tau*Fs;
t=[1/Fs:1/Fs:Nlen/Fs].';
SNR=6;
Fc=47.333333333333333333333;%
disp(['SNR = ' num2str(SNR)])
disp(['Fc = ' num2str(Fc)])
ini_phase = (rand(1)*2-1)*pi
y = 1.*exp(j*(2*pi*Fc*t+ini_phase));
Sinp = awgn(y,SNR,'measured');
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d)
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(11)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est = phase_pf_pi*Fs/2;
delta_F = Fc_est-Fc;
disp(['Оценка частоты: Fc_est = ' num2str(Fc_est)])
disp(['Погрешность оценки частоты: delta_F = ' num2str(delta_F)])
figure(12)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ сигнала после фильтра');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(13)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр сигнала после фильтра'),hold on;
% пофильтруем
h_len = round(Fs/Fc_est)*1;
hbpf = exp(j*2*pi*(Fc_est/Fs)*[1:h_len]);
Sinp_f = filter(hbpf,1,Sinp);

% еще разок
Sinp = Sinp_f;
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d);
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(21)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est2 = phase_pf_pi*Fs/2;
delta_F2 = Fc_est2-Fc;
disp(['Оценка частоты №2: Fc_est2 = ' num2str(Fc_est2)])
disp(['Погрешность оценки частоты: delta_F2 = ' num2str(delta_F2)])
figure(22)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ входного сигнала');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(23)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр вх сигнала'),hold on;
% % характеристики фильтра
% figure(31)
% plot(real(hbpf),'b- .'),grid on,hold on
% plot(imag(hbpf),'r- .'),grid on
% figure(32)
% freqz(hbpf,1)

для fontp:
хочу всё таки запустить метод маклеода.

Сигнал - это синусоида частотой 53.33Гц.
Дискретизируем его с частоой 960 Гц.
Берём 128 тоочек. (соответственно частотное разрешене DFT = 7.5Гц)
Делаем DFT - получаем размазаный спектр.

Теперь узнаю максимум - 8 бин DFT (т.е 60Гц)

В качестве вектора из трёх элементов отправляю на Маклеода массив из трёх значений соответсвтующих 7-му, 8-му и 9-му бинам DFT. (как положено - максимум и две точки окрестности)
Маклеод выдаёт -0.02
Стало быть смещение относительно максимума -0.02 бина. Это выходит 60 - 0.02*7,5 = 59,85 Гц

Объясните где я неправ. Может нужно разрешение увеличить?

Во всех спектральных оценивателях (Маклеоде, Квине) речь идёт о комплексной экспоненте. Непосредственно с действительной синусоидой не будет точно работать, поскольку синусоида - это две комплексных экспоненты. Поэтому для того, чтобы добраться до Маклеода (и всех других) нужно ещё сделать из реальной синусоиды комплексную, ну хотя бы тем же квадратурным детектором. Или подавить влияние второго пика широкими спектральными,например, гауссовыми окнами как описано по ссылке Стенфордского университета. Потом делать комплексное преобразование по типу DFT в пределах [-7.5 7.5]. В любом случае это комплексное DFT, а не действительное.

Средний бин будет соответствовать частотному нулю.
первый бин на стороне положительных частот соответствует нулю, восьмой 52.5

Кстати все три рассмотреных метода с действительной синусоидой работать не будут, не только спектральный. Здесь в теме рассмотрены в разных сообщениях всего 3 работающих почти оптимально метода - интерполяция спектра Фурье вблизи максимума, вычисление набега фазы через arg(сумма по T(S()*Sсопряженное())/T (автокорреляционный) и просто в лоб линейное моделирование набега фазы после квадратичного детектирования. Видите, эти формулы даже не напишешь в форуме, читайте лучше статьи....

Когда-то загружал статью, с вариантом перехода от реальной синусоиды к комплексной
http://electronix.ru/forum/index.php?act=A...st&id=37367

Большое спасибо!

вот ещё нашелся тот матлабовский тест по ссылке, "которого больше нет уже в Интернете"
Начните с него с комплексных экспонент, он должен сразу работать
Там три работающих оценщика параметров комплексной экспоненты Macleod, Quinn1, Quinn2

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

Но без квадратурного детектирования - это нептимально вычислительно. Всегда вычислительно можно сэкономить снося частоту в ноль, фильтруя полосу измерения и проводя сабсэмплинг - так чтобы сигнал был отцифрован уже оптимально, а не как получилось

Мне вот, непонятно, зачем "умножить на центральную комплексную экспоненту"? При таком сносе отклонение вправо/влево измеряемого вещественного сигнала от частоты 50Гц станут неотличимы друг от друга.

Мне кажется, что проще "умножить на комплексную экспоненту" с частотой 42.5Гц и т.о., снести в 0 левый участок спектра 42.5Гц после чего уже сделать FFT на интервале [0, 15].

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

А я делаю так. Беру выборку размером в два периода минимальной частоты (1/42,5)*2 с.
заетем ищу выборки возле переходов через ноль. С поможью интерполяции уточнаю точки перехода чарез ноль. Нахожу период - считаю частоту.

Fs = 6400;
Ts = 1/Fs;
f = 42.5;
t = 0:Ts:(1/42.5)*2-Ts;
A = 10;
Y1 = A*sin(2*pi*f*t); % модель сигнала
snr = 40; % отношение сигнал/шум
Y = awgn(Y1, snr);
z1 = find(diff(Y<0));%возвращаем номера элементов массива, предшествующих смене знака

T1 = [ Y(z1(2)),Y((z1(2)+1))]; %создаём массив из двух элементов окружающих смену знака
e1 = [Ts*z1(2)*1000, Ts*(z1(2)+1)*1000]; %массив временных значений
e1i = e1(1):0.008:e1(2)-0.008; %дробим временные значения на 25 частей
U1 = interp1(e1,T1,e1i,'line'); % линейная интерполяция в 25 раз вблизи нуля
T01 = (find(diff(U1<0))*0.008)+e1(1)-0.008; %первая опорная точка

%найдём вторую точку

T2 = [ Y(z1(4)),Y((z1(4)+1))]; %создаём массив из двух элементов окружающих смену знака
e2 = [Ts*z1(4)*1000, Ts*(z1(4)+1)*1000]; %массив временных значений
e2i = e2(1):0.008:e2(2)-0.008; %дробим временные значения на 25 частей
U2 = interp1(e2,T2,e2i,'line'); % линейная интерполяция в 25 раз вблизи нуля
T02 = (find(diff(U2<0))*0.008)+e2(1)-0.008; %вторая опорная точка
>> f0 = 1/((T02 - T01)/1000);

здесь и шум смоделирван. Как бы всё получается. Такие результаты в диссертации на мою тему.
Вот отрывок:

Судя по всему измеряете частоту прецессии (измерения магнитного поля земли). Какой точности в гаммах хотите добиться? Если ловите десятые доли гамма - очень проблематично и требует периода измерений в секундах, с\ш не меньше 40 дБ. Не зря ведь оверхаузера изобрели (для увеличения с\ш). Поищите по патентам GEM, решений много.

Нет. Я измеряю частоту основной гармоники в электрической сети.
по идее на АЦП поступает сигнал с трансформатора тока (как в отрывке диссертации выше).
Основная гармоника фильтруется полосовиком, затем применяется выше изложенный алгоритм.

Берутся две опорные точки (переходы через "0", а лучше - максимумы) затем интерполируется дабы уточнить момент. Так по точкам находят период, а затем и частоту. Моделирую в матлаб - всё нормально при С/Ш >40.
Почитайте отрывок из дисертации и мой пост.
Так в диссертации достигнута точность чуть выше чем 0,01Гц, чего более чем достаточно.

Хочется услышать коментарии спецов.

Я не про Вас, а про автора поста (куда то он пропал). А диссертацию почитаю, задача актуальная.

Скажите мне, ну причём тут разговоры про сигнал/шум и прочее... если намереваюсь профильтроовать сигнал полосовым дискретным фильтром, дабы выделить основную гармонику, диапазон изменения которой известен (42,5...57,5Гц)..... какой там шум на выходе????
Тогда измеряю просто по переходам через ноль...не выйдет?

Именно шум ограничивает в конечном итоге предельно достижимую точность, которую можно получить в принципе. Лучше этой границы не получится никакими методами, хуже может получиться в зависимости от особенностей примененного метода обработки.

Шум есть всегда и именно он ограничивает точность измерения для всех методов измерения сразу, причем науке известно как - это оценка максимального правдоподобия.
Вы почему-то пытаетесь упорно игнорировать этот факт.
Предлагаемый Вами метод очень неустойчив к шумам. Он дважды дифференцирует сигнал. А дифференцирование - операция
по отношению к шумам недопустимая, поскольку их усиливает.
Ваш метод будет работать только при низком уровне шума.
Что значит "будет работать" достойное с научной точки зрения?

Это значит что точность будет соизмерима с предельно возможной при данном уровне сигнал/шум и длине выборки -
а это оценка максимального правдоподобия Крамера-Рао. Если Вы на модельном сигнале сможете обосновать, что метод
стремится к этой точности хотя бы при больших SNR - то метод имеет научную ценность. В противном случае - увы...

Хотя с технической точки зрения может быть оправданы методы и не дающие предельную точность. Может быть простотой или ещё чем... Например так - "В принципе при возрастании длины выборки N ошибка всякого метода стремится к нулю, что не обязательно верно для SNR стремящегося к бесконечности. Предлагается метод, ошибка которого стремится к 0 при возрастании SNR к бесконечности, хотя он и не выходит на оценку CRLB.... Известны другие методы, достигающие этой оценки, но их сложность - не очень..." )) В конце концов, если не писать диссертации, а решать конкретную задачу, то можно просто сказать - заказчика точность устраивает. Однако нужно помнить, что чем ближе Вы можете подобраться к критерию CRLB тем метод практически всё равно предпочтительней, поскольку это гарантирует нужную точность при меньшем кол-ве отсчетов (за меньшее время), а временнОе разрешения важно уже тем, что параметры гармонического сигнала могут изменяться за время измерения, если это время велико. Т.е. на другом конце нестационарность сигнала.

В любом случае метод нужно моделировать и сравнивать с точностью максимального правдоподобия CRLB, поскольку это эталон
к которому нужно стремиться. Собственно не понятно, что Вы хотите ещё обсуждать...

случайно наткнулся на открытый проект по измерению магнитного поля , там применен интересный метод по измерению частоты http://www.jandspromotions.com/philips2005...ners/AR1746.htm

А не пробовали моделировать фильтр Герцеля? В смысле зависимость его характеристики от времени накопления при заданной частоте дискретизации? Там получаются весьма интересные колебания в зависимости от кратности частот дискретизации и измеряемой частоты. В указаном случае частота заранее определена достаточно точно и по поведению нескольких фильтров Герцеля (чтобы покрыть весь диапазон частот) можно оценить реальную частоту. Преимущество такого метода будут следующие:
1. Фильтрация исходного сигнала т.е. увеличение отношения сигнал/шум
2. Простота реализации БИХ фильтра
3. Маленькие вычислительные затраты
Думаю, вполне может получиться неплохой результат.

читаю и понимаю, что люди вопросы задают и ждут определанных ответов. Но если ответы не совпадают с их ожиданиями то они их просто игнорируют. Теперь по делу. В технике есть предельные оценки, например граница Шеннона скорости передачи данных. Также есть оценка Крамера-Рао точности измерения параметров сигнала. Так вот все что вы можете померить вы можете померить хуже чем эта оценка. Если вы хотите получить наилучшую из возможных оценок, то вы должны делать специальное устройство обработки, которое называеся оптимальным приемником. Надо сказать что когда мы говорим о наилучшей оценки, то надо понимать что мы задаем критерий, по которому данная оценка наилучшая. В данном случае - минимум ошибки измерения частоты при фиксированной выборке. Так вот в случае измерения частоты гармоники люди давно придумали оптимальные следящие контуры, которые подстраиваются под ваш сигнал с точностью до фазы и в режиме реального времени вы можете снимать с них данные о частоты вашего сигнала. Никакое fft не сравниться с контуром фапч по точности, поскольку контур фапч можно сделать на любой сетке частот все упирается в разрядность синтезатора частоты, и при частоте гармоники 50 Гц вытащить 0.01 Гц как нефиг нафиг, только надо обеспечить непрерывный поток данных с ацп, что при частоте 50 Гц тоже как нефиг нафиг. Короче разговоров про эти несчачтные 50 Гц на разных форумах пруд пруди уже не первый год, но что-то как-то все периливают из пустого в порожнее. За это время можно столько книжек прочитать и экспериментов поставить.