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

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

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

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

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

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

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

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

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

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

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



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

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

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

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

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

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


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

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

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

Спасибо!

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

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

Интересные результаты. А можно со случайной начальной фазой повторить? И разрядность тестового сигнала ограничить 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Гц/Мгц.

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

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

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


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


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

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