Измерение частоты основной гармоники (50 Гц) с точностью 0.01 Гц Опции

[Pridnya]

Здравствуйте!

Существует ли программный метод измерения частоты основной гармоники с точностью 0,01 Гц? Предполагается, что в системе есть антиалисинговый фильтр 0-1600 Гц, АЦП с частотой выборки 3200 Гц и микроконтроллер.

Сообщение отредактировал Pridnya - Sep 8 2015, 16:25

[Fat Robot]

Да. См методы MUSIC, ESPIRIT

http://mathworks.com/help/signal/ref/pmusic.html

Здравствуйте!

Существует ли программный метод измерения частоты основной гармоники с точностью 0,01 Гц? Предполагается, что в системе есть антиалисинговый фильтр 0-1600 Гц, АЦП с частотой выборки 3200 Гц и микроконтроллер.

[serjj]

MUSIC или ESPRIT могут и не залезть в микроконтроллер. Если искомая частота всегда около 50 Гц, то как вариант сделать сначала узкополосную фильтрацию (для MCU возможно лучше подойдёт IIR фильтр, но это детали), а потом частоту гармоники оценивать с помощью MLE (Maximum likelihood estimator), т.к. задача сведётся к оцениванию синусоиды на фоне аддитивного шума. Эта задача хорошо изучена, есть множество статей, в которых даны различные варианты реализации. Время на оценку (размер выборки) определяется с помощью границы Крамера-Рао (CRB) для минимального ожидаемого сигнал-шума. Формула CRB для MLE приводится в тех же статьях. Ресурс на реализацию намного меньше чем MUSIC/ESPRIT. И сам алгоритм проще.

Сообщение отредактировал serjj - Sep 8 2015, 17:44

[Fat Robot]

Для единственной синусоиды music сведется mle. А насчет 'не залезть' принцип простой ведь: без денег в кабаке делать нечего.

MUSIC или ESPRIT могут и не залезть в микроконтроллер. Если искомая частота всегда около 50 Гц, то как вариант сделать сначала узкополосную фильтрацию (для MCU возможно лучше подойдёт IIR фильтр, но это детали), а потом частоту гармоники оценивать с помощью MLE (Maximum likelihood estimator), т.к. задача сведётся к оцениванию синусоиды на фоне аддитивного шума. Эта задача хорошо изучена, есть множество статей, в которых даны различные варианты реализации. Время на оценку (размер выборки) определяется с помощью границы Крамера-Рао (CRB) для минимального ожидаемого сигнал-шума. Формула CRB для MLE приводится в тех же статьях. Ресурс на реализацию намного меньше чем MUSIC/ESPRIT. И сам алгоритм проще.

[serjj]

Для единственной синусоиды music сведется mle. А насчет 'не залезть' принцип простой ведь: без денег в кабаке делать нечего.

Статистически. MUSIC для одной частоты - это тот же MUSIC, что и для N частот. От svd не уйти. Только сигнальное подпространство будет всегда единичной размерности. В MLE же нужно только фильтр по входу, коррелятор, одно деление и арктангенс.

[Pridnya]

Спасибо ответившим! Много нового узнал.

[_pv]

Да. См методы MUSIC, ESPIRIT

может кто-нибудь в двух словах объяснить: для одной частоты чем это будет лучше чем просто посчитать Фурье, и потом около максимума наименьшими квадратами найти коэффициенты у функции Гаусса или параболы и положение максимума взять оттуда?
или сразу во временной области наименьшими квадратами искать в сигнале A*sin(wt+p).

[Guest_TSerg_*]

и потом около максимума наименьшими квадратами найти коэффициенты у или параболы и положение максимума взять оттуда?

Именно так и делал на 50 Гц: ПФ + поиск по параболе.
Вполне устраивало.

[blackfin]

Если искомая частота всегда около 50 Гц, то как вариант сделать сначала узкополосную фильтрацию, а потом частоту гармоники оценивать с помощью MLE..

В MLE же нужно только фильтр по входу, коррелятор, одно деление и арктангенс.

А зачем "фильтр по входу (..сделать сначала узкополосную фильтрацию)", если сам коррелятор уже и есть узкополосный фильтр?

[Fat Robot]

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

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

А так вы всё правильно рассказываете: наименьших, параболы и т.п. Только эту всю канитель надо как-то делать и проверять. Но может быть и заниматься этим не нужно из-за того, что всё равно ничего путного не выйдет.

Такие дела.

может кто-нибудь в двух словах объяснить: для одной частоты чем это будет лучше чем просто посчитать Фурье, и потом около максимума наименьшими квадратами найти коэффициенты у функции Гаусса или параболы и положение максимума взять оттуда?
или сразу во временной области наименьшими квадратами искать в сигнале A*sin(wt+p).

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

А зачем "фильтр по входу (..сделать сначала узкополосную фильтрацию)", если сам коррелятор уже и есть узкополосный фильтр?

[blackfin]

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

А кто сказал, что на входе коррелятора не аддитивная смесь сигнала и белого шума?

PS. И потом, коррелятор в моём понимании это Линейная Стационарная Система (с точностью до эффектов квантования).
То же, очевидно, касается и предполагаемого узкополосного фильтра на его входе.
А раз так, то обе системы можно поменять местами.
Но полоса пропускания коррелятора ~1/Tизм будет на порядки уже самого узкополосного фильтра, КМК..
Так зачем тогда делать двойную работу?

[Fat Robot]

"We can neither confirm nor deny the existence of the information requested..."

А кто сказал, что на входе коррелятора не аддитивная смесь сигнала и белого шума?

Все так. Только вот подавление у коррелятора вне его полосы пропускания слабенькое.

PS. И потом, коррелятор в моём понимании это Линейная Стационарная Система (с точностью до эффектов квантования).
То же, очевидно, касается и предполагаемого узкополосного фильтра на его входе.
А раз так, то обе системы можно поменять местами.
Но полоса пропускания коррелятора ~1/Tизм будет на порядки уже самого узкополосного фильтра, КМК..
Так зачем тогда делать двойную работу?

[serjj]

если сам коррелятор уже и есть узкополосный фильтр?

А что с чем мы будем коррелировать? В MLE обычно оценивается фрагмент АКФ случайного процесса, т.е. корреляция сигнала с задерженными копиями самого себя. Полоса такого фильтра определяется спектральным составом случайного процесса. Это не тоже самое, что корреляция входного процесса с некоторым гармоническим сигналом. С другой стороны, говорится, что MLE оптимальный метод для оценки комплексной экспоненты на фоне аддитвиного белого шума, т.е. фильтр нужен для подавления гармоник за пределами диапазона, в котором происходит оценка. Как сказал Fat Robot, это обеливание процесса на всякий случай. Если в системе может быть несколько гармоник или искомая частота варьируется в широком диапазоне, то да, MUSIC будет более разумным решением. Что же касается Фурье + интерполяция - почему бы и нет, только MLE по ресурсам выигрывает, иначе бы все ставили бы на feedforward синхронизацию скоростных модемов Фурье и радовались

Сообщение отредактировал serjj - Sep 9 2015, 09:13

[blackfin]

А что с чем мы будем коррелировать? В MLE обычно оценивается фрагмент АКФ случайного процесса, т.е. корреляция сигнала с задержанными копиями самого себя.

Если заранее известно, что сигнал гармонический, то на мой взгляд "корреляция с его задержанной копией" и есть просто корреляция с гармонической функцией, то есть с sin/cos.

[serjj]

Если заранее известно, что сигнал гармонический, то на мой взгляд и "корреляция с его задержанной копией" есть просто корреляция с гармонической функцией, то есть с sin/cos.

Заранее известно, что мы ищем сигнал с частотой близкой к 50 Гц в сигнале с символьной 3,2 кГц. Почему бы в таком сигнале не быть гармоникам на частотах 200, 300 Гц или каких либо ещё. Фильтр переведёт неизвестный случайный процесс к процессу вида "комплексная экспонента + белый шум", который является допустимым для работы MLE.

Сообщение отредактировал serjj - Sep 9 2015, 09:27

[blackfin]

Почему бы в таком сигнале не быть гармоникам на частотах 200, 300 Гц или каких либо ещё.

По условию задачи, гармоника одна. Всё, точка..

Понятно, что если это не так, то "простые" методы не работают. Но мы ведь не пытаемся "натянуть сову на глобус"?

[serjj]

Существует ли программный метод измерения частоты основной гармоники с точностью 0,01 Гц? Предполагается, что в системе есть антиалисинговый фильтр 0-1600 Гц, АЦП с частотой выборки 3200 Гц и микроконтроллер.

Измерение частоты основной гармоники, основной Карл (с)
Никто не говорит, что в диапазоне 0...1600 Гц она единственная.
А по поводу простоты - покажите, что для случая гармонический сигнал + шум есть более простой метод чем MLE и что эта оценка будет несмещённой и оптимальной.
Фильтр - эта та мелочь, которая позволит перейти от произвольного процесса к гармоника + шум, при условии, что мы априори знаем где лежит интересующая нас частота (это то, что нам дано по условию задачи).

[Fat Robot]

Ну вот хочет человек, чтобы во входном сигнале была единственная гармоника на фоне БГШ. Не нужно запрещать человеку хотеть. Дело за малым: приказать сигналу хорошо себя вести.

Никто не говорит, что в диапазоне 0...1600 Гц она единственная.

[blackfin]

Измерение частоты основной гармоники, основной Карл (с)

А, ну да.. это я проморгал..

А по поводу простоты - покажите, что для случая гармонический сигнал + шум есть более простой метод чем MLE и что эта оценка будет несмещённой и оптимальной.

Собсно, мне казалось, что ход мыслей примерно такой же.. "только в профиль"..

Поясню..

Вот, допустим, у нас есть "гармонический сигнал + шум".

Шум мы сразу с негодованием отбрасываем..

Тогда разыскиваемый сигнал можно записать в виде:

S(t) == S(n*∆t) = A*cos(ω*n*∆t+φ), где n = 0..N-1

Далее вычисляем сумму:

C_N = Σ{S(n*∆t)*e^{-j*ω₁*n*∆t}} = Σ{A*cos(ω*n*∆t+φ)*e^{-j*ω₁*n*∆t}}, где суммирование по всем: n = 0..N-1

В предположении, что частота ω₁ удовлетворяет соотношению: ω₁*N*∆t/2 = 2*pi*M, где M - целое,

находим отношение мнимой и действительной частей полученной суммы:

Im[C_N]/Re[C_N] = tg[ω*(N-1)*∆t/2 - ω₁*∆t/2 + φ] * tg[ω₁*∆t/2] / tg[ω*∆t/2].

В этом уравнении две неизвестных: ω и φ. Чтобы их исключить, сумму C_N нужно вычислить как минимум дважды для разных значений частоты ω₁.

Затем, с помощью MLE можно вычислить оба искомых параметра: ω и φ..

Понятно, что шум, гармоники сигнала и прочие эффекты могут сместить оценку, но сделав несколько измерений для разных значений частоты ω₁ их влияние можно сильно уменьшить, КМК..

Как-то так..

[attachment=95186:EstimateFrequency.doc]

[Fat Robot]

Вот кстати, наглядный пример того, что шумовой процесс должен быть белым для MLE. Если для этих "разных значений частоты ω₁" узкополосные процессы будут коррелированны, то это будет смещать оценку.

Это так.. общее замечание.

Понятно, что шум, гармоники сигнала и прочие эффекты могут сместить оценку, но сделав несколько измерений для разных значений частоты ω₁ их влияние можно сильно уменьшить, КМК..

[petrov]

но сделав несколько измерений для разных значений частоты ω₁ их влияние можно сильно уменьшить, КМК..

В несколько синк фильтров пролазит куча неизвестных гармоник, не очень понятно как можно сильно уменьшить их влияние.

[blackfin]

Если для этих "разных значений частоты ω₁" узкополосные процессы будут коррелированны, то это будет смещать оценку.

"Дьявол скрывается в деталях." В приложении к ЦОС можно даже сказать, что "Дьявол скрывается в значениях цифровых величин"..

Когда узкополосный фильтр вырезает из случайного процесса все лишние частоты, он точно так же смещает оценку,

если, конечно, этот фильтр не обладает как по волшебству идеально плоской частотной характеристикой..

Любая АЧХ имеющая наклон в точке измерения приведет к смещению оценки..

И таки да, для тех кто ещё не понял, C_N это обычные коэффициенты разложения сигнала S(t) в ряд Фурье..

[serjj]

В этом уравнении две неизвестных: ω и φ. Чтобы их исключить, сумму CN нужно вычислить как минимум дважды для разных значений частоты ω1.
...
Затем, с помощью MLE можно вычислить оба искомых параметра: ω и φ..
...
И таки да, для тех кто ещё не понял, CN это обычные коэффициенты разложения сигнала S(t) в ряд Фурье..

Вай! Сложно у вас всё!
Зачем считать Фурье на одной частоте, на другой частоте, получать систему с нелинейными уравнениями?
Почему бы не сделать вот так?
 Fast_carrier_offset_estimation.pdf ( 869.05 килобайт ) Кол-во скачиваний: 421

Вы пытаетесь оценить частоту используя сравнение входного сигнала с некоторым паттерном - делаете ли вы корреляцию или расчитываете Фурье для некоторой частоты - вы задаетесь своим значением ω1 и пытаетесь получить ω.
В основе спектрального анализа лежит дуальность АКФ и СПМ. MLE в статье (вариант реализации) основан на получении частоты из АКФ. Не нужно подбирать частоты ω1, для рассчёта потребуется только входной сигнал. Идея, как я её понял, проста: с одной стороны у нас есть оценки АКФ, полученные по данным, с другой - аналитическое значение АКФ для комплексной экспоненты. Сопоставляя оценку с аналитическим выражением, можно получить частоту сигнала. Увеличивая размер выборки и порядок эстиматора увеличиваем точность и (при росте порядка эстиматора) уменьшаем полосу захвата.

Тут же кстати и объяснение тому, что гармоники отличные от искомой вносят искажение в оценку. При оценке АКФ они вносят свой вклад, а для аналитической формы АКФ мы считаем, что в сигнале только одна гармоника. При рассчёте это несоответствие выливается в смещение оценки.

MUSIC, Capon и прочие методы, кстати, также опираются на АКФ, только в форме автоковариационной матрицы, в отличие от Фурье, в котором в сущности происходит сравнение сигнала с набором комплексных экспонент, частоты которых заданы и определяются размерностью.

Сообщение отредактировал serjj - Sep 9 2015, 11:17

[blackfin]

В несколько синк фильтров пролазит куча неизвестных гармоник, не очень понятно как можно сильно уменьшить их влияние.

Смысл в том, что если усреднить отклик от четного числа последовательно расположенных по частоте синк-фильтров,

то суммарный вклад "неизвестной гармоники" будет уменьшен, так как АЧХ синк фильтров меняет знак при переходе к следующему бину..

Но это так.. Гипотеза..

Зачем считать Фурье на одной частоте, на другой частоте, получать систему с нелинейными уравнениями?

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

Отличие состоит лишь в способе вычисления частоты.. Я просто предложил посчитать частоту "в лоб", минуя вычисление точки экстремума параболы.

И я отнюдь не призываю использовать предложенный способ в боевых условиях. Мне просто было интересно узнать мнение коллектива..

За сим, и спасибо!..

[petrov]

Смысл в том, что если усреднить отклик от четного числа последовательно расположенных по частоте синк-фильтров,

то суммарный вклад "неизвестной гармоники" будет уменьшен, так как АЧХ синк фильтров меняет знак при переходе к следующему бину

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

[blackfin]

Почему бы просто не проектировать сразу нужный FIR фильтр?

ОК, уговорили.. Пойду проектировать нужный фильтр..

[petrov]

MUSIC, Capon и прочие методы, кстати, также опираются на АКФ, только в форме автоковариационной матрицы, в отличие от Фурье, в котором в сущности происходит сравнение сигнала с набором комплексных экспонент, частоты которых заданы и определяются размерностью.

Не смущает, что в одном случае мы умножаем зашумлённый сигнал на зашумлённый сигнал, а в другом случае на чистенькую комплексную экспоненту?

[serjj]

Не смущает, что в одном случае мы умножаем зашумлённый сигнал на зашумлённый сигнал, а в другом случае на чистенькую комплексную экспоненту?

Я знаю два способа оценки АКФ: по данным и параметрический. В первом случае мы имеем данные, во втором - модель. Второй может дать нам АКФ, которая будет максимально близкой к аналитической, т.е. к истинной, но она требует знание а) количества сигналов, входящих в смесь и б) их частот. Если количество может быть дано априорно, то частоты мы не знаем наверняка, т.к. именно их мы и оцениваем. Так что при оценке частот АКФ мы вынуждены оценивать по данным с ограниченной размерностью и в случае, как например с AR экстраполировать. Параметрический подход применим когда например нужно посчитать оптимальный фильтр, если даны все параметры модели, полученные каким-либо другим способом.

Что-то упускаю? )

Сообщение отредактировал serjj - Sep 9 2015, 12:39

[petrov]

Что-то упускаю? )

Что же касается Фурье + интерполяция - почему бы и нет, только MLE по ресурсам выигрывает, иначе бы все ставили бы на feedforward синхронизацию скоростных модемов Фурье и радовались

Ну вот я так делаю и радуюсь, конечно всё БПФ смысла делать нет, только скользящие несколько бинов.

А вы говорите:

В MLE обычно оценивается фрагмент АКФ случайного процесса, т.е. корреляция сигнала с задерженными копиями самого себя.

Типа что-то сигнал у вас слишком хорош, давайте-ка мы его ещё на шум умножим.

Что-то упускаю? )

[serjj]

Типа что-то сигнал у вас слишком хорош, давайте-ка мы его ещё на шум умножим.

Автокорреляционная матрица сигнала, состоящего из детерменированной составляющей (например комплексная синусоида) и аддитивного белого шума, x(t) = s(t) + n(t), где n(t) - нормально распредёленный шум -
R_xx = R_ss + N * I, где N - мощность шума, а I - единичная диагональная матрица. Мы можем рассматривать первую строку этой матрицы как правую половину АКФ. Т.к. белый шум не коррелирует сам с собой, то он вносит своё искажение только на диагональные элементы матрицы, иначе говоря для нулевого лага. Оценка АКМ тем ближе к тому, что написано выше чем больше отчётов мы имеем. При этом она не зависит от уровня шума, при условии что шум белый.
Шум при этом разумеется искажает оценку параметра, но смещение, которое он вносит не увеличивается от того, что мы оцениваем АКФ только по данным, т.е. по смещённым копиям самого сигнала. Почему это так, я привёл выше.

Сообщение отредактировал serjj - Sep 9 2015, 13:40

[Pridnya]

Измерение частоты основной гармоники, основной Карл (с)
Никто не говорит, что в диапазоне 0...1600 Гц она единственная.
А по поводу простоты - покажите, что для случая гармонический сигнал + шум есть более простой метод чем MLE и что эта оценка будет несмещённой и оптимальной.
Фильтр - эта та мелочь, которая позволит перейти от произвольного процесса к гармоника + шум, при условии, что мы априори знаем где лежит интересующая нас частота (это то, что нам дано по условию задачи).

Товарищи ученые, доценты с кандидатами... Пока мало что понятно, но я стараюсь.
Реь идет о электрораспределительныой сети.
В сети всегда есть основная гармоника и гармоники до 13-й, но основная всегда около 50 Гц, например 45-55 Гц.
Все гармоники в сумме всегда значительно меньше основной.
50 Гц - 90%
100 Гц 0%
150 Гц 3%
200 Гц 0%
250 Гц 2%
...
Хочется выделить основную гармонику (цифровым фильтром) и измерить программно её частоту с точностью 0,01 Гц.

Что я могу:
С помощью ДПФ посчитать амплитуду основной гармоники.
С помощью БПФ посчитать амплитуды всех гармоник.
С помощью цифрового фильтра выделить полосу около 50 Гц и получить отфильтрованную выборку из 64-х точек.

[Tiro]

Реь идет о электрораспределительныой сети.
В сети всегда есть основная гармоника и гармоники до 13-й, но основная всегда около 50 Гц, например 45-55 Гц.
...
Хочется выделить основную гармонику (цифровым фильтром) и измерить программно её частоту с точностью 0,01 Гц.

Почему бы Вам не поинтересоваться, как это реализовано в регистраторах показателей качества электроэнергии?
Парма, к примеру.

Надо еще учесть, что в сети всегда присутствует АМ и ЧМ, фликкер и проч. Откройте ГОСТ 13109 и забудьте про цифровой фильтр. Нужен образцовый сигнал.

[petrov]

Pridnya

Выделяйте комплексным полосовым КИХ фильтром интересующий диапазон, АЧХ в полосе пропускания должна быть как можно более равномерная, гармоники(в том числе на отрицательных частотах) должны быть задавлены как можно сильнее, затем вычисляете несколько бинов ДПФ, далее по статье считаете частоту:

http://electronix.ru/forum/index.php?s=&am...t&p=1141831

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

[Милливольт]

Удалено, не понял задачи и сморозил чушь.

Сообщение отредактировал Милливольт - Sep 10 2015, 16:02

[Pridnya]

Почему бы Вам не поинтересоваться, как это реализовано в регистраторах показателей качества электроэнергии?
Парма, к примеру.

Надо еще учесть, что в сети всегда присутствует АМ и ЧМ, фликкер и проч. Откройте ГОСТ 13109 и забудьте про цифровой фильтр. Нужен образцовый сигнал.

Обязательно поинтересуюсь, спасибо за ссылку на прибор, который вы считаете образцом.
В сети много чего есть (в неизвестно какой электросети), они же все разные.
Упоминаний о наличии АМ и ЧМ фликкера сетях 6-10 кВ нигде не встречал, даже слово такое впервые слышу,
никогда не нужно было бороться с АМ и ЧМ фликкером в электросетях.
Исходя из используемой вами терминологии, ссылки на ГОСТ (электромагнитная совместимость) и предложением отказаться от цифровой обработки сигналов, предположу, что вы являетесь специалистом в области аналоговых систем связи.
Картинки из ГОСТ 13109 "колебания напряжения", "временное перенапряжение" (в вашей терминологии АМ-фликкер) и "перепады частоты сигнала"
(такой картинки нет, но в вашей терминологии ЧМ-фликкер), это как раз то, что мне нужно и чем я занимаюсь.
Перепады напряжений могут быть в диапазоне от 0 до 2-х кратного номинального, а отклонение частоты не более 10%.
Чтобы все это считать, используется аппаратная часть (в первом моем сообщении указано), микроконтроллер с ядром ARM и модулем FPU осуществляет обработку данных (это о доступной вычислительной мощности, ну, то есть, не компьютер под Windows с MatLab-ом это все обсчитывает).
ГОСТы на ЭМС необходимы только для сертификации готового изделия на ЭМС (как работает изделие в условиях помех в питающей сети ~220 вольт, помех в линиях дискретных входов, линий связи...), анализируются и регистрируются (всякие там аварийные журналы с параметрами аварийных режимов, векторные диаграммы токов и напряжений основной гармоники, ток нулевой, прямой, обратной последовательности основной частоты, напряжение нулевой последовательности основной частоты, углы между всеми этими векторами, среднеквадратичное значение суммы верхних гармоник тока нулевой последовательности...осциллограмма аварийного режима для аналоговых и цифровых входов в формате COMTRADE). Поэтому никак без цифровой фильтрации нельзя, везде фильтрация, свертка...И кстати, качество электроэнергии - это не основная функция, которая меня интересует.

Pridnya

Выделяйте комплексным полосовым КИХ фильтром интересующий диапазон, АЧХ в полосе пропускания должна быть как можно более равномерная, гармоники(в том числе на отрицательных частотах) должны быть задавлены как можно сильнее, затем вычисляете несколько бинов ДПФ, далее по статье считаете частоту:

http://electronix.ru/forum/index.php?s=&am...t&p=1141831

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

Спасибо за информацию! Чувствуется, что человек понимает о чем пишет. Приятно читать.
Вот симулинком не пользовался никогда. Наверное, придется попробовать. И за симулинк спасибо!

Если позволите...
Задача, поставленная автором не так проста... Во-первых, гармоники сетевой помехи далеко не всегда малы, а в ряде случаев имеют амплитуду больше, чем амплитуда основной частоты (sic!). Во-вторых: измерение с требуемой точностью требует и длинной реализации, и неиспользования любых аппроксимационных решений, т.к. в условиях шумов (априори неизвестных уровней) потенциальное смещение оценки может быть обусловлено именно ими.
А почему бы не использовать в этих условиях вариант синхронного фильтра (извиняюсь за самоцитирование ( http://www.tredex-company.com/ru/kvaziopti...i-i-ee-garmonik ), но не для режекции, а для выделения сигнала. Если делать ресэмплинг входных данных со смещающейся частотой сдвига (как пример ГКЧ), то отклик синхронного фильтра в какой-то момент пройдет максимум. И именно этот максимум будет соответствовать искомой частоте. Дополнительным «плюсом» такого решения будет использование энергии как четных, так и нечетных гармоник.
С вычислительной точки зрения это будет самое экономичное решение.

В сетях 6-10 киловольт в основном присутствует основная гармоника (50 Гц), поэтому по ней и работают. Это в сетях 0,4 киловольт (внутридомовые...) гармоник полно...
Вариант синхронного фильтра посмотрю. Тоже не встречал и не пользовался таким фильтром.

Сообщение отредактировал Pridnya - Sep 10 2015, 07:13

[Tiro]

Прошу не приписывать мне то, что не было сказано.

Обязательно поинтересуюсь, спасибо за ссылку на прибор, который вы считаете образцом.

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

Упоминаний о наличии АМ и ЧМ фликкера сетях 6-10 кВ нигде не встречал, даже слово такое впервые слышу,
никогда не нужно было бороться с АМ и ЧМ фликкером в электросетях.

Понятия "амплитудная модуляция", "частотная модуляция", "доза фликкера" применяются раздельно. Речь о борьбе с ними не шла.

Исходя из используемой вами терминологии, ссылки на ГОСТ (электромагнитная совместимость) и предложением отказаться от цифровой обработки сигналов, предположу, что вы являетесь специалистом в области аналоговых систем связи.

ГОСТ 13109-97 вообще-то называется ЭЛЕКТРИЧЕСКАЯ ЭНЕРГИЯ. НОРМЫ КАЧЕСТВА ЭЛЕКТРИЧЕСКОЙ ЭНЕРГИИ В СИСТЕМАХ ЭЛЕКТРОСНАБЖЕНИЯ ОБЩЕГО НАЗНАЧЕНИЯ, то есть речь никак не об ЭМС. Терминология взята в том числе из него.

Предложение отказаться от ЦОС не вносилось, а было предложение, альтернативное использованию цифрового фильтра. И с областью моих интересов Вы промахнулись.

(такой картинки нет, но в вашей терминологии ЧМ-фликкер), это как раз то, что мне нужно и чем я занимаюсь.

Это Ваша терминология.

Верно определяйте слова, и вы освободите мир от половины недоразумений. Рене Декарт

[Pridnya]

Прошу не приписывать мне то, что не было сказано.

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

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

http://www.parma.spb.ru/catalog/equipment/...tvo/parma_rk301
Что-то как-то не ярко, похоже очень старая разработка:
угла сдвига фаз между каналами напряжения (только РК03.01ПТ)
Оснащены интерфейсом Centronics для подключения матричного принтера и интерфейсом RS-232 для подключения ПК
Частота f Гц от 45 до 55 Гц ±0,02 Гц (абсолютная, при интервале усреднения 20 секунд)
4.4.2.1 Регистратор обеспечивает измерение текущих значений ПКЭ:
− частоты входного сигнала f от 40 до 70 Гц;
...
Таблица 1 – Нормируемые метрологические характеристики регистраторов
Частота f Гц от 45 до 55
...
4.4.2.1 Время установления рабочего режима – не более 60 с.
4.4.2.2 Потребляемая мощность регистратора не более 15 ВА.

Я считаю фазы всех векторов, направление мощности...
Использую изолированный интерфейс USB, Ethernet, т.е. к любому ноутбуку без всяких переходников могу подключиться.
А мне нужно круче чем ±0,02 Гц, чтобы 50,00 Гц "показометром" показывать (без усреднения за минуту), ну, т.е., чтобы быстро среагировать на изменение частоты.
Измеряем в диапазоне 40-70, а погрешность нормиркется в диапазоне 45-55 Гц.
И куда ж там 15 Вт на энергопотребление уходит, микрокомпьютер что-ли засунули с крякнутым MatLab-ом.

Не годится.

Сообщение отредактировал Pridnya - Sep 10 2015, 09:23

[rudy_b]

Парма - это, как раз, образец безграмотной разработки.

Стандартный способ точного определения частоты основной гармоники сети такой (примерно).
1. Выборка с частотой порядка 25 кГц точек (24 бита) за, примерно, 300 мс (это параметры реального прибора класса 0.05 с частотным диапазоном порядка 8 кГц).
2. На выборку накладывается окно Гаусса с аккуратно подобранными параметрами. Делается FFT (само FFT - float, обработка - double). При этом провал бин между соседними гармониками оказывается порядка 80 дБ.
3. По МНК (парабола или, лучше, гаусс, по 5 - 7 соседним бинам с весами) определяется частота максимума. Если она сильно смещена относительно центра выбранных бин - операция повторяется с коррекцией центра.

При этом точность определения частоты основной гармоники не хуже 0.005 Гц даже при 10% гармониках.

[Pridnya]

Стандартный способ точного определения частоты основной гармоники сети такой (примерно)...

Спасибо!!! Буду пробовать сам метод, не смотря на то, что у меня АЦП 16 бит и частота выборки 3200 Гц (могу поднять до 6400).

Сообщение отредактировал Pridnya - Sep 10 2015, 11:35

[Tiro]

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

Кажется, кажется.. Еще у меня мультиметр Fluke, а осциллоскоп простенький LeCroy. Давайте определяйте по приборам, какой из этих фирм я сотрудник?

Парма - это, как раз, образец безграмотной разработки.

Стандартный способ точного определения частоты основной гармоники сети такой (примерно).
1. Выборка .... за, примерно, 300 мс

При этом точность определения частоты основной гармоники не хуже 0.005 Гц даже при 10% гармониках.

Вот и выросло поколение воинствующих инженеров, не отличающих точность от разрешающей способности, мгновенную частоту от частоты основной гармоники

[Pridnya]

Кажется, кажется.. Еще у меня мультиметр Fluke, а осциллоскоп простенький LeCroy. Давайте определяйте по приборам, какой из этих фирм я сотрудник?


Вот и выросло поколение воинствующих инженеров, не отличающих точность от разрешающей способности, мгновенную частоту от частоты основной гармоники

Тепрь-то можно написать хоть про Снеговика с ведром на голове, на содеянное это не повлияет. ))

На данном этапе можно любую терминологию использовать, лишь бы Src и Dst понимали друг друга, а для других - хоть белый шум. Кстати, единого языка с единой терминологией для нучного мира так и не создали.

Сообщение отредактировал Pridnya - Sep 10 2015, 19:58

[rudy_b]

Вот и выросло поколение воинствующих инженеров, не отличающих точность от разрешающей способности, мгновенную частоту от частоты основной гармоники

Точность и разрешающая способность - это совершенно разные вещи, указанная цифра дает именно значение погрешности, т.е. среднеквадратичное отклонение. Разрешающая способность на 2 порядка выше.

Понятия "мгновенная частота" не существует.

При периодическом сигнале частота основной гармоники совпадает с измеренной за любые несколько периодов. Иначе сигнал периодическим ( в строгом смысле этого слова) не является и говорить о частоте основной гармоники бессмысленно.

[Fat Robot]

Дарю. Пользуйтесь.

Понятия "мгновенная частота" не существует.

[rudy_b]

Спасибо, конечно, но это не метрологическое, а чисто теоретическое понятие, которое применяется только на чистом синусе (ну, в крайнем случае, при узкополосном синале, не содержащем гармоник). Непосредственному измерению не подлежит, поскольку определяется через другой, опять же, непосредственно не измеряемый параметр - скорость изменения фазы, которую опять же невозможно определить по отсчетам АЦП, поскольку, даже при моночастотном сигнале, его амплитуда на момент измерения неизвестна. Для понимания этого полезно почитать литературу, указанную в списке литературы под номером [3] в вашей же ссылке (Финк Л. М. Сигналы, помехи, ошибки…).

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

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

[serjj]

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

f(t) = d/dt(arg(y(t))), простейший частотный дискриминатор, известный со времён старого доброго ЧМ, нет? )

А ещё оказывается в современном ЦОС есть целый класс алгоритмов, т.н. DIFM, например
 GE_47.pdf ( 358.18 килобайт ) Кол-во скачиваний: 386

Вполне себе практические вещи делают.

[blackfin]

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

А ещё был такой чудак Гейзенберг, так тот вообще всем доказывал, что ошибка "∆F" при оценке измерении мгновенной частоты "F" равна бесконечности..

И даже формулу в свое оправдание привел:

∆E*∆t ≥ ћ/2,

и, с учетом того, что:

∆E == ћ*∆ω == ћ*2*pi*∆F,

такую же формулу для принципа неопределенности:

∆F*∆t ≥ 1/(4*pi).

Так что да, оценить мгновенную частоту сложно, а вот в какую-нить формулу запихнуть - легко!

[rudy_b]

f(t) = d/dt(arg(y(t))), простейший частотный дискриминатор, известный со времён старого доброго ЧМ, нет? )

Нет, вы не отличите изменение фазы от изменения амплитуды.

А ещё оказывается в современном ЦОС есть целый класс алгоритмов, т.н. DIFM, например
 GE_47.pdf ( 358.18 килобайт ) Кол-во скачиваний: 386

Вполне себе практические вещи делают.

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

[serjj]

Нет, вы не отличите изменение фазы от изменения амплитуды.

Да вы что!
Ну расскажите мне, где вы в формуле f(t) = d/dt(arg(y(t))) увидили изменение амплитуды? И как же это у нас FM радио тогда ловит?..

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

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

Сообщение отредактировал serjj - Sep 11 2015, 14:46

[rudy_b]

Да вы что!
Ну расскажите мне, где вы в формуле f(t) = d/dt(arg(y(t))) увидили изменение амплитуды? И как же это у нас FM радио тогда ловит?..

Ну тогда уж сначала вы расскажите мне, как вы измеряете arg(y(t))?

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

Как я уже говорил ранее, понятие "мгновенная частота" в метрологии отсутствует. Оно есть только в головах пишущих формулы теоретиков.

Попробуйте почитать книжку, указанную в посте

...Финк Л. М. Сигналы, помехи, ошибки…

[serjj]

Ну тогда уж сначала вы расскажите мне, как вы измеряете arg(y(t))?

Арктангенсом, а вы?

И, только тогда, можно оценить "мгновенную" частоту.
...
Как я уже говорил ранее, понятие "мгновенная частота" в метрологии отсутствует. Оно есть только в головах пишущих формулы теоретиков.

То у вас её можно оценить, то она только в "головах у теоретиков". Вам приводят конкретные статьи и примеры, а вы говорите какие-то общие фразы и не хотите их количественно подтвердить. Речь же о точном предмете )

[rudy_b]

Арктангенсом, а вы?

Что-то я не знаю приборов измеряющих арктангенс произвольного аналогового сигнала, может просветите?

Кроме того, арктангенс - разрывная функция и брать производную от него...

Несерьезно как-то, лучше книжку почитайте.

[анатолий]

Хочется выделить основную гармонику (цифровым фильтром) и измерить программно её частоту с точностью 0,01 Гц.

Что я могу:
С помощью ДПФ посчитать амплитуду основной гармоники.
С помощью БПФ посчитать амплитуды всех гармоник.
С помощью цифрового фильтра выделить полосу около 50 Гц и получить отфильтрованную выборку из 64-х точек.

Когда-то в 1993г. такую задачу решал для измерения токов в рельсах метро на контроллере 8051 на частотах 50, 75, 125 ... Гц.
Решалась просто: сигнал умножался на косинус, синус 50 Гц и накапливался как скользящее среднее. Результат - как фильтрация узким фильтром. Если сигнал строго =50 Гц. то выходной сигнал - постоянная, если 49 Гц - вых. сигнал 1 Гц, если 51 Гц - сигнал 1 Гц, но крутящийся в противоположную сторону. И т.п. Соответственно, если на выходе синусоида с периодом 100 сек. то сигнал 50,01 Гц.

[Krys]

И т.п. Соответственно, если на выходе синусоида с периодом 100 сек. то сигнал 50,01 Гц.

дак это надо тогда 100 секунд ждать, пока хоть 1 период наловишь? А чтобы принять надёжное решение, нужен ещё и не один период? Видимо хотелось бы побыстрее...

[serjj]

Что-то я не знаю приборов измеряющих арктангенс произвольного аналогового сигнала, может просветите?

Кроме того, арктангенс - разрывная функция и брать производную от него...

Просвещайтесь на здоровье

[blackfin]

Просвещайтесь на здоровье

А чего это они в формулах (4) и (7) после второго знака равенства "arctan" выкинули?

[serjj]

А чего это они в формулах (4) и (7) после второго знака равенства "arctan" выкинули?

Арктангенс равен аргументу комплексного числа, записанного в (3) как две отдельные квадратуры, обычное приравнивание. Если кому-то интересно, то вот о процессе угловой модуляции. Тогда всё сразу становится понятнее )

Сообщение отредактировал serjj - Sep 14 2015, 07:24

[blackfin]

Если кому-то интересно, то вот о процессе угловой модуляции. Тогда всё сразу становится понятнее )

Да читали мы все это. И даже экзамены сдавали по радиолокации.. Просто очень давно..

[rudy_b]

И где же тут прямое измерение "мгновенной" частоты? После фильтров ФНЧ?
Вы просто подтвердили мои слова:

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

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

И, как обычно, вы перепутали реально существующую и оцениваемую по достаточно длинной выборке разность фаз двух сигналов с мифической фазой произвольного ОДИНОЧНОГО аналогового сигнала.

Измеряется не фаза, а разность фаз, и за несколько периодов.

И, если вы прочитали всю приведенную ссылку, в ней, опять же, подтверждают мои слова о разрывности дифференцируемого вами арктангенса

Кроме того, арктангенс - разрывная функция и брать производную от него...

Поэтому, в реальной обработке, арктангенс старательно исключают.

Расскажите мне, как определить фазу сигнала по одному отсчету АЦП? Ну ладно, учитывая нарисованные вами производные - пусть будет целых три последовательных отсчета с интервалом 1 мксек: -1В; 0В; +1В.
Ну, хотя бы, любимую "мгновенную" частоту определите?

Не получается? Тогда прекратите увиливать, лучше читайте книжки и думайте.

[Fat Robot]

Нет ничего проще: На основании именно этих полученных от АЦП отсчетов можно смело сказать, что на его входе присутствует синусоидальный сигнал с частотой 250 kHz, c амплитудой 1В. Центральному отсчету соответствует фаза 0. Мы считаем, что амплитуда сигнала нормирована под полный "раскрыв" АЦП +/-1В.

Расскажите мне, как определить фазу сигнала по одному отсчету АЦП? Ну ладно, учитывая нарисованные вами производные - пусть будет целых три последовательных отсчета с интервалом 1 мксек: -1В; 0В; +1В.
Ну, хотя бы, любимую "мгновенную" частоту определите?

[thermit]

rudy_b:
подтверждают мои слова о разрывности дифференцируемого вами арктангенса

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

Сообщение отредактировал thermit - Sep 14 2015, 13:52

[serjj]

И, если вы прочитали всю приведенную ссылку, в ней, опять же, подтверждают мои слова о разрывности дифференцируемого вами арктангенса

И прочитал и в своё время сделал в железке, не поверите
Его там исключают потому что далее делается дифференцирование и можно вместо d/dt(arg(y(t))) рассчитывать конечное выражение и собрать эквивалентую схему. Никто не запрещает сделать схему "в лоб" с использованием арктангенса + unwrap и всё будет работать (проверено на практике, работает, alas!).
Да и где ж разрыв?

Поэтому, в реальной обработке, арктангенс старательно исключают.

Приехали. Преобразование координат и демодуляцию PSK-N созвездий наверное тоже без арктангеса делают, ага.

Расскажите мне, как определить фазу сигнала по одному отсчету АЦП? Ну ладно, учитывая нарисованные вами производные - пусть будет целых три последовательных отсчета с интервалом 1 мксек: -1В; 0В; +1В.
Ну, хотя бы, любимую "мгновенную" частоту определите?

А мне не нужно работать на уровне АЦП. И ТС'у никто не предлагал работать на уровне АЦП и по одному отчёту. Есть сигнал после гетеродина и ФНЧ, и он комплексный. Для него понятие фазы определено строго математически и применимо к каждому символу. А раз определено понятие фазы, то значит можно оценить мгновенную частоту. Точность такой оценки будет зависеть от сигнал-шума. Процесс демодуляции ЧМ сигналов привёл как простой и распространённый пример.

Что ещё Вы хотите тут опровергать?

Сообщение отредактировал serjj - Sep 14 2015, 14:20

[mcheb]

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

Наверное, подразумевалась неоднозначность.

[blackfin]

Как я уже говорил ранее, понятие "мгновенная частота" в метрологии отсутствует. Оно есть только в головах пишущих формулы теоретиков.

Можно вполне строго доказать невозможность измерения "мгновенной частоты" на практике.

Для этого нужно воспользоваться формулой для вычисления частоты по трем значениям гармонической функции S(t) из соседней темы:

f = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

Это абсолютно точная формула, позволяющая получить точное значение частоты "f" при условии,

что нам известны точные значения гармонической функции S(t) в моменты времени t₁,t₂ и t₃:

S₁ == S(t₁),
S₂ == S(t₂),
S₃ == S(t₃),

где:

Δt == t₃ - t₂ = t₂ - t₁,

и при этом:

f*(t₃ - t₁) < 1.

Чтобы строго доказать невозможность измерения "мгновенной частоты" на практике, необходимо вычислить полный дифференциал функции трех переменных:

f(S₁,S₂,S₃) = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

А именно:

df == (δf/δS₁)*dS₁ + (δf/δS₂)*dS₂ + (δf/δS₃)*dS₃,

где частные производные функции "f" равны:

δf/δS₁ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

δf/δS₂ = +[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(S₁+S₃)/(2*S₂²)].

δf/δS₃ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

Соответственно, полный дифференциал равен:

df == -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(dS₁+dS₃-(S₁+S₃)*dS₂/S₂)/(2*S₂)].

Теперь мы можем найти погрешность измерения частоты "σ_F" связанную с погрешностью измерения сигнала "σ_S":

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*sqrt[2+(S₁+S₃)²/S₂²]/[2*abs(S₂)],

или:

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[(2*S₂)²-(S₁ + S₃)²]]*sqrt[2+(S₁+S₃)²/S₂²].

Знаменатель в этой формуле примечателен тем, что при стремлении t₃ к t₁ он стремится к нулю как const*(t₃-t₁)²

Как следствие, при стремлении t₃ к t₁ для любой погрешности "σ_S" измерения амплитуды сигнала S(t), погрешность измерения "мгновенной частоты" "σ_F" устремляется к бесконечности как σ_S/(t₃-t₁)²..

[rudy_b]

Нет ничего проще: На основании именно этих полученных от АЦП отсчетов можно смело сказать, что на его входе присутствует синусоидальный сигнал с частотой 250 kHz, c амплитудой 1В. Центральному отсчету соответствует фаза 0. Мы считаем, что амплитуда сигнала нормирована под полный "раскрыв" АЦП +/-1В.

Смелый вы человек. Но с чего вы взяли, что это синус? Боженька на ушко шепнул? А может это линейно нарастающее напряжение? Или пила, трапеция и т.д. с произвольным периодом?
Или вообще непериодический сигнал?

Но, даже если предположить чистый синус, то с чего вы решили, что амплитуда сигнала равна 1В?
Вот набор пар амплитуда/частота который дает именно эти отсчеты. Ваш вариант - это только один из них.
Амплитуда ----------- Частота --------- Напряжение при t=1 мксек
1.000000E+00 -- 2.500000E+05 1 -- Ваш вариант
1.500000E+00 -- 1.161398E+05 1
2.250000E+00 -- 7.329944E+04 1
3.375000E+00 -- 4.787579E+04 1
5.062500E+00 -- 3.164613E+04 1
7.593750E+00 -- 2.101973E+04 1
1.139063E+01 -- 1.399046E+04 1
1.708594E+01 -- 9.320293E+03 1
2.562891E+01 -- 6.211555E+03 1
3.844336E+01 -- 4.140452E+03 1
5.766504E+01 -- 2.760129E+03 1
8.649756E+01 -- 1.840034E+03 1
и далее ...

А если взять частоты более 250 кГц (не совсем корректно, но чего только в жизни не бывает) - то снова получим неограниченный набор пар частота/амплитуда, котрые дадут те же отсчеты.

Если вы собираетесь заниматься гаданием на кофейной гуще, то вы не туда попали, эта тема связана с метрологией.

А мне не нужно работать на уровне АЦП. И ТС'у никто не предлагал работать на уровне АЦП и по одному отчёту. Есть сигнал после гетеродина и ФНЧ, и он комплексный.

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

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

Действительно не совсем корректно выразился. Разрывность появляется при попытке определить именно фазовый угол (что, при грамотной обработке, практически не требуется) по соотношению Re и Im компонент применяя функцию арктангенса. Нужна, как минимум, atan2, но и с ней получаемый фазовый угол будет не линейно нарастающим (если есть постоянная разность частот), а скакать при переходе через +/-пи.

Наверное, подразумевалась неоднозначность.

Это как минимум.

[serjj]

f = [1/[pi*(t3 - t1)]]*arccos[(S1 + S3)/(2*S2)].
...
Как следствие, при стремлении t3 к t1 для любой погрешности "σS" измерения амплитуды сигнала S(t), погрешность измерения "мгновенной частоты" "σF" устремляется к бесконечности как σS/(t3-t1)2..

А тыкните носом в раздел анализа и теории оценок, где написано, что надо так считать. Просто интересно
В моём понимании f = [1/[pi*(t3 - t1)]]*arccos[(S1 + S3)/(2*S2)] при t3-t1->0 примет вид acos(1)/0, что даёт ноль в числителе и ноль в знаменателе, нормальное поведение для производной.
Кстати ваш вывод будет справедлив только для гармонического сигнала, но всё равно любопытно.

С другой стороны для идеального (нет шума) действительного гармонического сигнала можно написать такой простой скрипт:

Код

clear all;
Ns = 1e4;
fs = 100e5;
Ts = 1/fs;
t = 0:Ts:(Ns-1)*Ts;
f0 = 1333.756;
x = sin(2*pi*f0*t);
x1 = sqrt(1 - x.^2);
ph = angle(x1 + 1i*x);
index = randi([2 Ns]);
fcalc = abs(ph(index)-ph(index-1))/(Ts * 2 * pi);
str = sprintf('real frequency = %E -- estimated frequency = %E', f0, fcalc);
disp(str);

который вычисляет мгновенную частоту по двум точкам (по одной точке можно оценить фазу, но для действительного сигнала она не имеет смысла, т.к. в сущности может быть любой в зависимости от выбора начальной точки; для частоты достаточно двух точек). Легко проверяется в матлабе.
Увеличивая частоту дискретизации, мы сокращаем время между двумя точками, с другой стороны разница напряжений в этих точках тоже уменьшается, т.о. мы приближаемся к определению мгновенной частоты. Если шума нет, то мгновенная частота будет определяться с точностью до fs, если добавить шум, то погрешность, которую он внесёт, превысит погрешность от fs (при fs >> искомой частоты). Погрешность такого эстиматора уже будет оцениваться с помощью CRB.

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

Увиливаете тут только вы. Если вы хотите, чтобы ваши агрументы были приняты, опишите это математически или, если удобнее скриптом. Мы не на уроке литературы. А за одним посмотрите, чем отличается неоднозначность обратных тригонометрических функций от понятия разрыва и непрерывности. Это довольно познавательно.

Сообщение отредактировал serjj - Sep 15 2015, 08:56

[Fat Robot]

Ваши рассуждения довольно любопытны для человека, который декларирует, что хоть сколько-то занимается метрологическим оборудованием. Сначала вы ставите условие и просите получить по данным из условия оценку частоты. После того, как оценка частоты получена, вы дополняете условия, и оценка оказывается не верной и/или не может быть получена. Знаменитый принцип ведения дискуссии: "так да не так".

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

А второе: любой оценщик, как использующий 3 отсчета, так и 333 отсчета, можно обмануть и придумать бесконечное множество "а если", которые заставят его показывать не то, что вы подразумевали. Это не повод опускать руки. Нужно ограничить область применения доп. условиями, предоставить больше априорной информации о входном воздействии к началу измерений.

Выводы простые в общем-то:
Можно ли оценивать частоту сигнала по 3м отсчетам? Можно, если имеется другая необходимая априорная информация или предположения о сигнале.
Можно ли "обмануть" любой, наперед заданный оценщик? Можно, если выйти из базиса гипотез об оцениваемом сигнале.

А существование понятия "мгновенная частота".. Угловое положение вращающегося диска можно измерить. Время тоже.

Смелый вы человек. Но с чего вы взяли, что это синус? Боженька на ушко шепнул? А может это линейно нарастающее напряжение? Или пила, трапеция и т.д. с произвольным периодом?
Или вообще непериодический сигнал?

Но, даже если предположить чистый синус, то с чего вы решили, что амплитуда сигнала равна 1В?
Вот набор пар амплитуда/частота который дает именно эти отсчеты. Ваш вариант - это только один из них.

[...]

А если взять частоты более 250 кГц (не совсем корректно, но чего только в жизни не бывает) - то снова получим неограниченный набор пар частота/амплитуда, котрые дадут те же отсчеты.

Если вы собираетесь заниматься гаданием на кофейной гуще, то вы не туда попали, эта тема связана с метрологией.

[blackfin]

А тыкните носом в раздел анализа и теории оценок, где написано, что надо так считать. Просто интересно

Читайте про преобразования случайной величины, Якобиан и проч.

Или: Г.Корн Т.Корн, стр. 561, Гл. 18.5-4, "Функции от случайных величин. Замена переменных."

Ну и имейте ввиду, что случайная величина "f" есть однозначная функция трех независимых случайных величин: f = f(S1,S2,S3)..

[thermit]

С другой стороны для идеального (нет шума) действительного гармонического сигнала можно написать такой простой скрипт:
Код
clear all;
Ns = 1e4;
fs = 100e5;
Ts = 1/fs;
t = 0:Ts:(Ns-1)*Ts;
f0 = 1333.756;
x = sin(2*pi*f0*t);
x1 = sqrt(1 - x.^2);
ph = angle(x1 + 1i*x);
index = randi([2 Ns]);
fcalc = abs(ph(index)-ph(index-1))/(Ts * 2 * pi);
str = sprintf('real frequency = %E -- estimated frequency = %E', f0, fcalc);
disp(str);

Для неопределенной амплитуды скрипт, естественно мертвый. 3 точки нужно минимум.

[serjj]

Для неопределенной амплитуды скрипт, естественно мертвый. 3 точки нужно минимум.

Да, точно. Ну и на практике никто так не делает (x1 = sqrt(1 - x.^2)), если нужно получить аналитический сигнал, то ставят Гильберта. Там никакой зависимости от амплитуды + класс входных сигналов не ограничен только гармоническими. Ну или гетеродин + ФНЧ - если диапазон измерения априори задан.

Сообщение отредактировал serjj - Sep 15 2015, 12:47

[Pridnya]

Можно вполне строго доказать невозможность измерения "мгновенной частоты" на практике.

Для этого нужно воспользоваться формулой для вычисления частоты по трем значениям гармонической функции S(t) из соседней темы:

f = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

Это абсолютно точная формула, позволяющая получить точное значение частоты "f" при условии,

что нам известны точные значения гармонической функции S(t) в моменты времени t₁,t₂ и t₃:

S₁ == S(t₁),
S₂ == S(t₂),
S₃ == S(t₃),

где:

Δt == t₃ - t₂ = t₂ - t₁,

и при этом:

f*(t₃ - t₁) < 1.

Единственная формула, которую сразу можно посчитать в C++ (IDE NetBeans, компилятор MinGW ). Два результата: для 16-ти разрядного АЦП и для идеального АЦП (см. результат в конце C++ кода). Для хорошего АЦП метод точный, т.е. если фильтрануть основную гармонику, то будет хорошо.

Другие способы (MUSIC, MLE - Метод максимального правдоподобия...), по крайней мере, на первый взгляд очень сложные. Какие-то для матлаба, какие-то в общем виде многоэтажные формулы. Кто такие формулы пишет? Шифротекст для потомков. Неужели один я такой тупой, что не понимаю.

Код

#include <cstdlib>
#include <iostream>   // std::cout
#include <math.h>

using namespace std;

#define SAMPLES 64
//short Data[SAMPLES] = {0}; // АЦП 16-ти разрядный.
float Data[SAMPLES] = {0}; // идеальный АЦП, но его же нет (сигнал 16-ти разрядного АЦП как-то нужно улучшить???).

int main(int argc, char** argv) {

// Создаем сигнал частотой 50 Гц (первая гармоника): косинус, амплитуда 1000, фаза 0.
// таблица содержит 64 выборки на интервале 0.02 сек (период промышленной частоты 50 Гц).    
for(int i=0;i<SAMPLES;i++)
    {Data[i]=1000*cos(2*M_PI*i*1/SAMPLES+2*M_PI*0/360.0F);} // 1000*sin(1wt+0).

// Три равноотсоящих точки выборки (первые три точки).
float S1 = Data[0];
cout << "S1 = " << S1 << endl;  // 1000
float S2 = Data[1];
cout << "S2 = " << S2 << endl;  // 995
float S3 = Data[2];
cout << "S3 = " << S3 << endl;  // 980
float dT = 1.0F/3200.0F;
cout << "dT = " << dT << endl;  // 0.0003125 сек (период между выборками).
float t1 = 0;
cout << "t1 = " << t1 << endl;  // 0 сек.
float t3 = 2*dT;
cout << "t3 = " << t3 << endl;  // 0.000625 сек.

//f = [1/[pi*(t3 - t1)]]*acos[(S1 + S3)/(2*S2)] // Оригинальная формула.
float f = (1/(M_PI*(t3 - t1)))*acos((S1 + S3)/(2*S2)); // Заменил квадратные скобки на круглые.

cout << "f = " << f << endl;    // 51.078 Гц ??? для массива short Data[64]
cout << "f = " << f << endl;    // 50.0001 Гц для массива float Data[64]

return 0;

Сообщение отредактировал Pridnya - Sep 15 2015, 14:23

[blackfin]

Единственная формула, которую сразу можно посчитать в C++.

"Сложные проблемы всегда имеют простые, легкие для понимания неправильные решения."

[serjj]

Другие способы (MUSIC, MLE - Метод максимального правдоподобия...), по крайней мере, на первый взгляд очень сложные. Какие-то для матлаба, какие-то в общем виде многоэтажные формулы. Кто такие формулы пишет? Шифротекст для потомков. Неужели один я такой тупой, что не понимаю.

Вы никак не улучшите ваш сигнал оцифрованный 16-битным АЦП. При таком соотношении символьной и искомой частот (3200/50) аргумент функции acos -> 1, а значит сама функция - к 0, т.е. некоторому малому числу. Далее вы делите на 2pi*dt, другое малое число. С вычислительной точки зрения мало устойчивая задача. Любое незначительное искажение сигнала внесёт значительное искажение в оценку. Пропустив сигнал через АЦП вы добавили шум квантования, который внёс ошибку в ваш эстиматор. А теперь представте - в реальности при оцифровке стараются сделать так, чтобы оцифрованный аддитивный шум (тепловой шум, шум канала, кароче говоря неустронимый физический шум) был много больше шума квантования - условие необходимое для сохранения SNR на том уровне, который был до оцифровки. И посчитайте, какую ошибку внесёт он. Будете неприятно удивлены, т.к. даже незначительный шум даст полностью бесмыссленные значения частот. И даже узкополосный полосовой фильтр на 50 Гц вряд ли существенно улучшит результат.
Так что вам придётся применять какой либо из "сложных" методов.

Сообщение отредактировал serjj - Sep 15 2015, 15:40

[blackfin]

При таком соотношении символьной и искомой частот (3200/50) аргумент функции acos -> 1, а значит сама функция - к 0, т.е. некоторому малому числу.
Далее вы делите на 2pi*dt, другое малое число. С вычислительной точки зрения мало устойчивая задача.

Так Вы, оказывается, не поняли основной идеи предложенной формулы..

Я же там где-то раньше указывал, что для получения точного значения частоты должно выполняться условие: abs(S₁ + S₃) < abs(S₂).

Это требование фактически означает, что для вычисления "f" годятся не любые моменты времени t₁,t₂ и t₃, а только те,

в которых значения S₁ и S₃ находятся "вблизи" нуля и при этом значение S₂ "автоматически" оказывается вблизи максимума функции S(t).

Это примерно соответствует вычислению частоты "вручную", то есть когда мы просто измеряем моменты времени в которые гармонический сигнал S(t) пересекает ось абсцисс,

потом находим соответствующий этим моментам времени период гармоники "T" и уже затем вычисляем её частоту f = 1/T. Просто в предложенном способе не нужно точно вычислять моменты времени, когда S(t) = 0.

PS. А про шумы полностью согласен. Формула с arccos применима только для сферических функций в вакууме..

[Pridnya]

...
PS. А про шумы полностью согласен. Формула с arccos применима только для сферических функций в вакууме..

Формула сама интересная и, если сигнал будет чистый, то должна дать хороший результат.
Попробую сигнал с выхода АЦП обработать НЧ фильтром с КИХ и полосой 75Гц -3 дБ, 100 Гц -40 дБ, 64 коэффициента float.
А там можно разные способы попробовать. Что-то похожее советовал мне petrov.

Pridnya

Выделяйте комплексным полосовым КИХ фильтром интересующий диапазон, АЧХ в полосе пропускания должна быть как можно более равномерная, гармоники(в том числе на отрицательных частотах) должны быть задавлены как можно сильнее, затем вычисляете несколько бинов ДПФ, далее по статье считаете частоту:

http://electronix.ru/forum/index.php?s=&am...t&p=1141831

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

Только вот не пойму, откуда бины возьмутся после фильтрации сигнала фильтром с частотой среза 75 Гц, вроде там должен быть один бин, т.е. 50 Гц. Или я что-то не понимаю. Или он имеет ввиду посчитать 1024 точечное ДПФ от сигнала на интервале 0,02 секунды и по тем бинам посчитать?

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

 

[blackfin]

А там можно разные способы попробовать. Что-то похожее советовал мне petrov.

Эту задачу на форуме решали уже десятки раз..

См. например:

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

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

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

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

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

И там ещё в окрестностях есть что почитать..

[Pridnya]

Эту задачу на форуме решали уже десятки раз..

См. например:

И там ещё в окрестностях есть что почитать..

Спасибо! В любом случае потребуется много экспериментов, чтобы понять, как измерять, какой метод применить.

[EvgenyNik]

Pridnya, если с частотой семплирования и разрешающей способностью всё понятно, то на каком минимально допустимом интервале времени (в периодах основной гармоники) Вы бы хотели измерить частоту?

[Pridnya]

Pridnya, если с частотой семплирования и разрешающей способностью всё понятно, то на каком минимально допустимом интервале времени (в периодах основной гармоники) Вы бы хотели измерить частоту?

Хотелось бы измерять за 1, максимум за 5 периодов. Не знаю пока сколько достаточно и на сколько достоверно будет при большом интервале, а вдуруг там какое качание частоты, поэтому хотелось бы за один период.

[Corner]

Делов то. Фурье на 8к. Самая яркая палка это она))). Если предварительно сузить полосу прикинув частоту и лисперсию с помощью какого-нибудь капчюр модуля МК. Можно сделать перенос на 0 вокруг, то можно понизить частоту дискретизации и Фурье станет менее требовательным к числу точек и скорости вычислений.

[serjj]

Делов то. Фурье на 8к. Самая яркая палка это она))). Если предварительно сузить полосу прикинув частоту и лисперсию с помощью какого-нибудь капчюр модуля МК. Можно сделать перенос на 0 вокруг, то можно понизить частоту дискретизации и Фурье станет менее требовательным к числу точек и скорости вычислений.

Почитайте внимательно тему.
Фурье 8192 вы собрались на мк сделать?
На частоте дискретизации, обозначенной ТС'ом в 3200 Гц, 8192 точки это 2560 мс. 1 - 5 периодов 50 Гц сигнала это 20 - 100 мс.

[Pridnya]

Делов то. Фурье на 8к. Самая яркая палка это она))). Если предварительно сузить полосу прикинув частоту и лисперсию с помощью какого-нибудь капчюр модуля МК. Можно сделать перенос на 0 вокруг, то можно понизить частоту дискретизации и Фурье станет менее требовательным к числу точек и скорости вычислений.

Накидайте ка на Си, а я код прогоню на ARM и оценим:
а) реализуемость;
б) время выполнения кода в целевом устройстве;
в) глубину бездны знаний.

Сообщение отредактировал Pridnya - Sep 21 2015, 14:00

[Tiro]

Накидайте ка на Си, а я код прогоню на ARM и оценим:
а) реализуемость;
б) время выполнения кода в целевом устройстве;
в) глубину бездны знаний.

Сколько можно уже фигней страдать. Вот сюда посмотрите база сигнала
А потом положите B=1, дельта F=0.01 и вычислите дельта t. Затем примите как данность, что Вы проектируете показометр и рассчитывайте значение частоты хоть по 3 точкам на период хоть до 10 знаков после запятой.

[Pridnya]

Сколько можно уже фигней страдать. Вот сюда посмотрите база сигнала
А потом положите B=1, дельта F=0.01 и вычислите дельта t. Затем примите как данность, что Вы проектируете показометр и рассчитывайте значение частоты хоть по 3 точкам на период хоть до 10 знаков после запятой.

Δ t = B/Δf = 1/0.01 = 100 (секунд?). Не мойму, что-ли 100 секунд предлагается ждать? Подобный метод уже предлагали (если на выходе биение 0,01 Гц, то частота 49,99 Гц).

[Santik]

Хотелось бы измерять за 1, максимум за 5 периодов. Не знаю пока сколько достаточно и на сколько достоверно будет при большом интервале, а вдуруг там какое качание частоты, поэтому хотелось бы за один период.

Код

            М                  Погрешность, Гц
            1                               1
            3                               0.4
            5                               0.25
           10                              0.06
           15                              0.025
           20                              0.015
           30                              0.0065

Сигнал: синусоида f=50 Гц со случайной начальной фазой. fd=44100 Гц.
М- число периодов синусоиды.

[Pridnya]

Код

            М                  Погрешность, Гц
            1                               1
            3                               0.4
            5                               0.25
           10                              0.06
           15                              0.025
           20                              0.015
           30                              0.0065

Сигнал: синусоида f=50 Гц со случайной начальной фазой. fd=44100 Гц.
М- число периодов синусоиды.

Не понятно, как вы получили результат (матлаб посчитал с помощью dll-ки или вы на С код написали)? Понимаю, что частота дискретизации для оцифровки звука. Хотелось бы попробовать ваш метод расчета на С или С++ при отклонении частоты.

Интервал 5 периодов сети 50 Гц - это 0,1 секунды. При частоте дискретизации 44100 Гц, получаем 4410 точек выборки на том интервале. А дальше что делать?

Сообщение отредактировал Pridnya - Sep 22 2015, 07:16

[Santik]

Не понятно, как вы получили результат (матлаб посчитал с помощью dll-ки или вы на С код написали)? Понимаю, что частота дискретизации для оцифровки звука. Хотелось бы попробовать ваш метод расчета на С или С++ при отклонении частоты.

Интервал 5 периодов сети 50 Гц - это 0,1 секунды. При частоте дискретизации 44100 Гц, получаем 4410 точек выборки на том интервале. А дальше что делать?

Это скорее для иллюстрации, что на 5 периодах погрешность получится 0.25Гц, а вам надо 0.01 :-)
Алгоритм сложный, в МК точно не поместится.
Советую не заморачиваться, поставить 1-бит АЦП (0 - компаратор) и счётчик.
На этом принципе все протонные магнитометры работают. (Погрешность 0.005 Гц, частота 2-3 кГц, время 1-2 сек)

[Pridnya]

Это скорее для иллюстрации, что на 5 периодах погрешность получится 0.25Гц, а вам надо 0.01 :-)
Алгоритм сложный, в МК точно не поместится.

Да мне хотя бы 0,25 Гц программно получить, а там я бы посмотрел что дальше делать. Все говорят "возьмите...", а вот сколько точек на каком интервале взять и как дальше, вот что главное.
Хотя бы с пяти периодов начнем (сигнал на интервале 0,1 секунды). Частота выборки пока 3200 гц. Сколько точек БПФ брать на этом интервале?
И правильно ли я понимаю, что я получу разрешение по частоте 10 Гц? 64точки*5периодов = 320 точек. 3200Гц/320 = 10 Гц. Так?

Советую не заморачиваться, поставить 1-бит АЦП (0 - компаратор) и счётчик.
На этом принципе все протонные магнитометры работают. (Погрешность 0.005 Гц, частота 2-3 кГц, время 1-2 сек)

У меня использовался компаратор, модуль захват-сравнение, можно вытащить 0,01 Гц, с той же звуковой карты компа подаю и вижу 49.00 Гц, 50,00 Гц...Но мне интересно, как это сделать программно (без матлабов и прочих коммерческих продуктов для настольконго компа) для микроконтроллера.
Заодно и пойму, как получить высокое разрешение по частоте.

Сообщение отредактировал Pridnya - Sep 22 2015, 08:12

[Santik]

64точки*5периодов = 320 точек. 3200Гц/320 = 10 Гц. Так?

Да так. При данной частоте дискретизации и N=320 получится именно такое разрешение по частоте.
Весьма наивно полагать, что проведя параболу через 3 точки, близкие к соответствующему максиммуму спектра, можно получить хорошее разрешение по частоте.

[Pridnya]

Да так. При данной частоте дискретизации и N=320 получится именно такое разрешение по частоте.
Весьма наивно полагать, что проведя параболу через 3 точки, близкие к соответствующему максиммуму спектра, можно получить хорошее разрешение по частоте.

Уже лучше, можно смоделировать: задать чистую синусоиду 50Гц на интервале 0,1 сек, а потом от этого интервала взять 1024-х точечное ДПФ для нескольких частот, близких к искомой, а там посмотреть какая кривая получится.

А вот как считать? В смысле спектр будет симметричный относительно искомой частоты (чтобы парабола получилась), т.е. отрицательные и положительные частоты? Или от 0 до частоты выборки Fs.

Библиотека CMSIS_DSP для ARM выдает результат от 0 до Fs, т.е. спектр симметричный относительно Fs/2.

Сообщение отредактировал Pridnya - Sep 22 2015, 08:42

[Santik]

Уже лучше, можно смоделировать: задать чистую синусоиду 50Гц на интервале 0,1 сек, а потом от этого интервала взять 1024-х точечное ДПФ для нескольких частот, близких к искомой, а там посмотреть какая кривая получится.
А вот как считать? В смысле спектр будет симметричный относительно искомой частоты (чтобы парабола получилась), т.е. отрицательные и положительные частоты? Или от 0 до частоты выборки Fs.
Библиотека CMSIS_DSP для ARM выдает результат от 0 до Fs, т.е. спектр симметричный относительно Fs/2.

Можно сразу задать частоты 49,50,51 Гц и найти значения спектра на этих частотах. Потом провести параболу и найти её максимум. Это будет первое приближение по частоте. Затем взять +-0.5 Гц от найденой частоты, снова найти значения спектра и снова проводить параболу, находить ее максимум, делить интервал на 2 (т.е. до 0.25 Гц) и т.д. Такая вот итерационная процедура... Но это долго получится.

[Pridnya]

Можно сразу задать частоты 49,50,51 Гц и найти значения спектра на этих частотах. Потом провести параболу и найти её максимум. Это будет первое приближение по частоте. Затем взять +-0.5 Гц от найденой частоты, снова найти значения спектра и снова проводить параболу, находить ее максимум, делить интервал на 2 (т.е. до 0.25 Гц) и т.д. Такая вот итерационная процедура... Но это долго получится.

Буду пробовать. Спасибо! Открытие для себя сделал: оказывается можно больше одного периода сети анализировать, мне доступно 16, 32 ,64, 128, 256, 512, 1024, 2048 и 4096-точечное БПФ (спектр посмотрю). А я почему-то думал, что БПФ нужно на интервале периода брать, т.е. на интервале 0,02 секунды 64, 128, 256 и т.д. точек.

PS: Смущает только, как оценить скорость изменения частоты при больших интервалах. Но это вопрос следующий.

Сообщение отредактировал Pridnya - Sep 22 2015, 09:09

[petrov]

Но мне интересно, как это сделать программно (без матлабов и прочих коммерческих продуктов для настольконго компа) для микроконтроллера.
Заодно и пойму, как получить высокое разрешение по частоте.

Да хоть исходники тут выложить для вашего микроконтроллера, к пониманию это вас не приблизит. Ключ к пониманию - декомпозиция, упрощение, визуализация, симулинк для этого отличнейший инструмент. Давно бы уже комплексный дискретный тон сгенерировали бы, подали на 3 КИХ фильтра, соответствующие бинам ДПФ, посчитали бы частоту по статье Эрика Якобсена.

[Pridnya]

Да хоть исходники тут выложить для вашего микроконтроллера, к пониманию это вас не приблизит. Ключ к пониманию - декомпозиция, упрощение, визуализация, симулинк для этого отличнейший инструмент. Давно бы уже комплексный дискретный тон сгенерировали бы, подали на 3 КИХ фильтра, соответствующие бинам ДПФ, посчитали бы частоту по статье Эрика Якобсена.

До меня только сегодня дошло, как увеличить разрешение по частоте, используя преобразование Фурье и длинные интервалы, но это при условии, что частота основной гармоники не меняется. А вот что делать если она медленно меняется и нужно еще знать скорость изменения частоты? Вроде, интервал нужно уменьшать, но тогда разрешение по частоте уменьшится. Фурье не подходит, вроде. Может вейвлеты?

А суть идеи я понял (я не понимал, что вы считаете бинами, я их гармониками называл и думал, откуда там бины рядом с частотой основной гармоники, а там абстракция), и статья Эрика Якобсена у меня отложена, спасибо, вы мне её и давали.

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

http://www.ibm.com/developerworks/ru/libra...rl_1/index.html

А про "давно бы уже...", так есть ещё и другие вопросы.

Сообщение отредактировал Pridnya - Sep 22 2015, 11:51

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

 

[petrov]

До меня только сегодня дошло, как увеличить разрешение по частоте, используя преобразование Фурье и длинные интервалы, но это при условии, что частота основной гармоники не меняется. А вот что делать если она медленно меняется и нужно еще знать скорость изменения частоты? Вроде, интервал нужно уменьшать, но тогда разрешение по частоте уменьшится. Фурье не подходит, вроде. Может вейвлеты?

Должно дойти, что нужна точность, а не разрешение.

Картинку посмотреть и в командной строке можно, видел статью, там график в текстовом виде выводится. Главное знать алгоритм, а посмотреть - дело второе.

А Ньютону и листка бумаги хватит. Что не понятно в алгоритме, комплексный тон не можете сгенерировать, не понятно как 1 бин ДПФ вычисляется, слишком сложна самая простая из возможных формула Якобсена?

[Pridnya]

Что не понятно в алгоритме, комплексный тон не можете сгенерировать, не понятно как 1 бин ДПФ вычисляется, слишком сложна самая простая из возможных формула Якобсена?

Элементарно:
берем отсчеты АЦП и у каждого отсчета мнимую часть приравниваем нулю, далее обрабатываем комплексный входной сигнал библиотечными функциями CMSIS DSP.

Можно и без комплексного руками посчитать в тетрадке (для каждого бина ДПФ). Раньше я брал 64-х точечное ДПФ на интервале 0,02 сек и получал частоту 1 гармоники 50 Гц, второй 100 Гц и т.д. И не понял про какие бины возле 50-ти Гц вы мне говорили (откуда думаю там бины, если сигнал всего один, 50 Гц, если задать только один сигнал) .Теперь понятно.

Сообщение отредактировал Pridnya - Sep 22 2015, 12:18

[serjj]

До меня только сегодня дошло, как увеличить разрешение по частоте, используя преобразование Фурье и длинные интервалы, но это при условии, что частота основной гармоники не меняется. А вот что делать если она медленно меняется и нужно еще знать скорость изменения частоты? Вроде, интервал нужно уменьшать, но тогда разрешение по частоте уменьшится. Фурье не подходит, вроде.

На мк всю жизнь использовали рекурсивный алгоритм Герцеля, если нужно было в real-time считать отдельные частотные бины. Есть его модификация на тему скользящего Фурье преобразования (sdft), про это есть и на форуме и в гугле. Хотите отслеживать флуктуации частоты в окресностях 50 Гц - заменяете вычисление отдельных бинов dft на sdft, а далее - как вам советовали. Получаете скользящую модификацию алгоритма. Разрешение при этом будет зависеть от задержки преобразования (задержка sdft аналог количества точек dft), а точность - от сигнал-шума, нижнее значение которого вы нам кстати так и не озвучили (с этого нужно начинать проектирование любого эстиматора).

[Pridnya]

На мк всю жизнь использовали рекурсивный алгоритм Герцеля, если нужно было в real-time считать отдельные частотные бины. Есть его модификация на тему скользящего Фурье преобразования (sdft), про это есть и на форуме и в гугле.

Я его и использовал раньше, когда нужно несколько нечетных гармоник найти.

Хотите отслеживать флуктуации частоты в окресностях 50 Гц - заменяете вычисление отдельных бинов dft на sdft, а далее - как вам советовали. Получаете скользящую модификацию алгоритма. Разрешение при этом будет зависеть от задержки преобразования (задержка sdft аналог количества точек dft), а точность - от сигнал-шума, нижнее значение которого вы нам кстати так и не озвучили (с этого нужно начинать проектирование любого эстиматора).

dft - ДПФ?
sdft - буква s что означает? Может sign (знак?) Google не знает.
https://www.google.ru/search?ie=UTF-8&h...&gws_rd=ssl

А отношение сигнал/шум большое, поэтому не критично.

Сообщение отредактировал Pridnya - Sep 22 2015, 12:57

[serjj]

dft - ДПФ?
sdft - буква s что означает? Может sign (знак?) Google не знает.
https://www.google.ru/search?ie=UTF-8&h...&gws_rd=ssl
А отношение сигнал/шум большое, поэтому не критично.

dft=discrete fourier transform
s=sliding.
Вы заявили что хотите точность 0,01 Гц. Для такой точности нет некритичного SNR. Это большая точность. Тем более, если вы хотите отслеживать изменение частоты во времени -> статичный спектральный анализ уже отметается.

Я его и использовал раньше, когда нужно несколько нечетных гармоник найти.

Вот тоже самое. dft и sdft имеют одинаковый результат каждые N отчётов. sdft позволяет отслеживать изменение комплексных значений выбранных частотных бинов во времени с задержкой N. Обрабатывая их соответствующим образом (интерполяция например), вы получите слежение за частотой выбранной гармоники.
Это если вам так хочется через Фурье делать. Слежение можно сделать и чисто временными способами. Зависит от ситуации. Если вы уже знакомы с Герцелем, то возможно через него вам будет проще.

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

Сообщение отредактировал serjj - Sep 22 2015, 13:17

[Pridnya]

dft=discrete fourier transform
s=sliding.
Вы заявили что хотите точность 0,01 Гц. Для такой точности нет некритичного SNR. Это большая точность. Тем более, если вы хотите отслеживать изменение частоты во времени -> статичный спектральный анализ уже отметается.

Спасибо!

[rudy_b]

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

Я уже писал про параметры измерения частоты (25600 Гц, 320 мсек, окно Гаусса, FFT 8к), но вам такая точность, похоже, не нужна.
Вот результаты грубого моделирования с разной частотой самплинга на том же интервале 320 мсек Fft8k(25600Гц), Fft1k(3200Гц), Fft256(800Гц).
На входе синус 51 Гц + 10% третьей гармоники + 1% случайного шума.
1. FFT для трех случаев (8к, 1к и 256 точек)

2. Обрезали все точки выше 75 Гц (чтобы не хватать вторую и третью гармоники, можно отбросить и точки ниже 25 Гц, но поленился это делать) и провели простейший фиттинг по гауссу. Показано для 8к.


А вот численные результаты фиттинга. Обратите внимание на xc (положение центра, т.е. измеренная частота синуса) и его ошибку
Value Standard Error
A8k y0 0.23205 0.06367
A8k xc 50.99956 4.13929E-4
A8k w 9.05248 8.87164E-4
A8k A 28998.67448 2.76099
A8k sigma 4.52624 4.43582E-4
A8k FWHM 10.65848 0.00104
A8k Height 2555.94038 0.20737

Value Standard Error
A1k y0 1.03749 0.16077
A1k xc 50.99733 0.00105
A1k w 9.04926 0.00224
A1k A 28990.07709 6.9703
A1k sigma 4.52463 0.00112
A1k FWHM 10.65469 0.00264
A1k Height 2556.0913 0.52375

Value Standard Error
A256 y0 3.45667 0.4462
A256 xc 50.99655 0.0029
A256 w 9.04252 0.00622
A256 A 28920.49775 19.3375
A256 sigma 4.52126 0.00311
A256 FWHM 10.64675 0.00733
A256 Height 2551.85819 1.45431