[Pridnya 0]

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

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

 

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

Изменено 22 сентября, 2015 пользователем Pridnya

[petrov 7]

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

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

 

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

[Pridnya 0]

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

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

 

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

 

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

 

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

 

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

Изменено 22 сентября, 2015 пользователем Pridnya

[petrov 7]

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

 

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

 

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

 

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

 

[Pridnya 0]

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

Элементарно:

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

 

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

Изменено 22 сентября, 2015 пользователем Pridnya

[serjj1333 0]

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

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

[Pridnya 0]

На мк всю жизнь использовали рекурсивный алгоритм Герцеля, если нужно было в 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

 

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

Изменено 22 сентября, 2015 пользователем Pridnya

[serjj1333 0]

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

Изменено 22 сентября, 2015 пользователем serjj

[Pridnya 0]

dft=discrete fourier transform

s=sliding.

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

 

Спасибо! :rolleyes:

 

[rudy_b 1]

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

При понижении частоты самплинга (и, соответственно, снижении числа точек 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

[blackfin 23]

При понижении частоты самплинга (и, соответственно, снижении числа точек FFT) растет влияние шумов, соседних гармоник и т.п.

Вот результаты грубого моделирования с разной частотой самплинга на том же интервале 320 мсек Fft8k(25600Гц), Fft1k(3200Гц), Fft256(800Гц).

На входе синус 51 Гц + 10% третьей гармоники + 1% случайного шума.

1. FFT для трех случаев (8к, 1к и 256 точек)

2. Обрезали все точки выше 75 Гц и провели простейший фиттинг по гауссу.

Есть сильное подозрение, что "1% случайного шума" в первом случае (Fft8k-25600Гц) у Вас "размазан" в полосе 12800 Гц, во втором (Fft1k-3200Гц) - в полосе 1600 Гц, а в третьем (Fft256-800Гц) - в полосе 400 Гц. Соответственно, и спектральная плотность шума в полосе 75 Гц тоже отличается как (1/12800) : (1/1600) : (1/400)..

[Santik 0]

...На входе синус 51 Гц + 10% третьей гармоники + 1% случайного шума.

 

    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

 

Интересные результаты. А можно со случайной начальной фазой повторить? И разрядность тестового сигнала ограничить 10 бит?

 

[Tiro 0]

Интересные результаты. А можно со случайной начальной фазой повторить? И разрядность тестового сигнала ограничить 10 бит?

 

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

 

P.S.

1 Вместо случайной фазы можно взять пример с кратными соотношениями частот, проще считать. И ситуацию рассмотреть ближе к реальности, когда отрезки сигнала сшиты последовательно. Частота 50+1 Гц. Формируется сигнал длительностью 1 сек. Очевидно, что на этом интервале реализация закончится как для частоты 50 Гц, так и для частоты 51 Гц. АЦП даст 3200 отсчетов в любом случае. Затем делим сигнал на отрезки, либо по периодам частоты 50 Гц (к примеру, 5 периодов х 0,02 с = 0,1 с; либо по периодам БПФ, например 3 х БПФ-1024, остальные отсчеты оставим на след цикл). Получаем либо 10 значений аппроксимации, либо 3. (Хотя для оценки ошибки число измерений мало в любом случае). Считаем, находим матожидание, ско. Задаемся доверительной вероятностью и получаем точность результата. (Причем в условиях отсутствия систематической и случайной погрешностей измерения). Можно взять интервал 10 сек, результатов будет побольше.

 

2 А еще у ТС кроме 10 бит АЦП еще и кварц на МК скорее всего обычный. То есть порядки примерно 10^-4 начальное отклонение и 10^-6 на градус К температурный дрейф. Ну и джиттер частоты 10Гц/Мгц.

[rudy_b 1]

Есть сильное подозрение, что "1% случайного шума" в первом случае (Fft8k-25600Гц) у Вас "размазан" в полосе 12800 Гц, во втором (Fft1k-3200Гц) - в полосе 1600 Гц, а в третьем (Fft256-800Гц) - в полосе 400 Гц. Соответственно, и спектральная плотность шума в полосе 75 Гц тоже отличается как (1/12800) : (1/1600) : (1/400)..

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

Я хотел лишь прояснить основные моменты для ТС и на высокую точность не претендую - это просто иллюстрация. Если ему нужно - проделает все сам, надеюсь, что объяснил достаточно понятно.

 

А 1% шума - это на всякий случай, в том числе и учет погрешности и разрядности АЦП.

Интересные результаты. А можно со случайной начальной фазой повторить? И разрядность тестового сигнала ограничить 10 бит?

Многократно проверялось - при правильном окне гаусса влияние фазы пренебрежимо мало (порядка 10^-5). Но можете сами проверить, сигма гаусса - 900 для 8к и соответственно меньше для 1к и 256.

Разрядность АЦП - это тот же 1% шума. Он, конечно, случайный, а не систематический, но, если есть желание - легко считается любой. Хотя неслучайный шум - это, фактически, нелинейные искажения порождающие верхние гармоники, а они отсекаются.

 

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

В данном случае стандартные отклонения - это честно посчитанные ошибки определения каждого из параметров при фиттинге. Сигнал генерировался независимо для каждого случая и 1% вклад rnd() ессно варьировался.

 

2 А еще у ТС кроме 10 бит АЦП еще и кварц на МК скорее всего обычный. То есть порядки примерно 10^-4 начальное отклонение и 10^-6 на градус К температурный дрейф. Ну и джиттер частоты 10Гц/Мгц.

Уход кварца это уже систематическая ошибка. А джиттер... Такой кварц еще поискать нужно. Хотя, на плохом генераторе...

[Pridnya 0]

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

При понижении частоты самплинга (и, соответственно, снижении числа точек FFT) растет влияние шумов, соседних гармоник и т.п.

 

Я уже писал про параметры измерения частоты (25600 Гц, 320 мсек, окно Гаусса, FFT 8к), но вам такая точность, похоже, не нужна.

Вот результаты грубого моделирования с разной частотой самплинга на том же интервале 320 мсек Fft8k(25600Гц), Fft1k(3200Гц), Fft256(800Гц).

На входе синус 51 Гц + 10% третьей гармоники + 1% случайного шума.

...

Спасибо за метод и результаты моделирования! :laughing: