Как-то плохо вникаете, в статье же всё предельно ясно и кратко изложено. Формула(3) не даёт смещения. Идеальный синус - комплексный синус. Сначала надо отладить на идеальном и получить результат не хуже чем в статье. Реальный синус уже состоит из суммы двух комплексных на положительной и отрицательной частотах, можно задавить комплексный синус на отрицательной частоте фильтром, не вносящим искажений в области 45 Гц до 55 Гц, тогда формула(3) будет работать.

А насколько идеален этот "сигнал с генератора сигналов"?

И какова неравномерность этого "антиалисингова фильтра" в полосе от 45 Гц до 55 Гц?

Вникаю по мере понимания материала. Все же нужно понять, протестировать.
Отлаживаю на идеальном сигнале, программно на С и С++.
Суть метода Эрика Якобсена ясна, но не ясно как с помощью МНК подобрать оптимальную кривую, да еще чтобы формула (или способ) работала в диапазоне частот.
Он там формулы привел, считайте. А я бы хотел визуализировать способ в диапазоне частот, с таблицами и сравнением, чтобы там побольше материала было. В общем эта публикация не раскрывает всей сути.
Слишком много математики нужно, чтобы это все потестирвоать. Нужен не один человеко-день. Поэтому подождем с результатами по Эрику Якобсену.

Далее вы пишете

Что значит реальный синус? Сигнал с выхода датчика напряжения или трансформатора тока, нагруженного на резистор, можно считать "реальным" в вашей терминологии?
Для меня реальные сигналы, это те, которые можно задать приборами, фильтровать, измерять.

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

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



Для тестирвоания формул я программно задаю синус определенной частоты, с шагом 0,01 Гц. Т.е. можно сказать, что генератор сигналов идеальный.

Имею возможность проверить на реальном сигнале любой метод, пригодный для практических расчетов: 0-100 вольт, 50,00 Гц. Сигнал получается с выхода ЦАП и усиливается прецизионными ОУ.
Антиалисинговые фильтры (ФНЧ Баттерворта 4 порядка ) имеют частоту среза 1600 Гц, поэтому в полосе частот герц до 400 не ослабляют сигнал. Посчитаны в FilterPro.

Вы вообще представляете себе, что вычисляет ПФ от конечного отрезка сигнала? Какой предполагается исходный сигнал?
Отрезок реального? Нет. ПФ вычисляет коэффициенты от периодического повторения Вашего отрезка сигнала.
Это точно Ваш сигнал? Вряд ли

Наложение окон во временной области позволяет уменьшить краевые эффекты при ПФ. Как-то так.

А по теории рядов Фурье - Воробьев НН Теория рядов. УДК.517.52(075.8)

Предполагается, что исходный сигнал существует от минус бесконечности до плюс бесконечности.
А вот анализируется отрезок реального, выборка из некоего количества точек, например 1024 для 1024-х точечного ДПФ или БПФ.

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



Это известное заученное наизусть предложение. Оно как бы позволяет уменьшить ошибку после того как мы перешли "от минус бесконечности до плюс бесконечности" к выборке "1024 точки", т.е. выдрали её из пространственно временного континиума.


Я хоть и про теорию рядов вообще не спрашивал и у вас в частности тоже, но просмотрел монографию.

Мои впечатления:
Об авторе ничего нет, зато еть бюрократическое "Допущено минобразования СССР для студентов ВТУзов", каких ВТУзов, не указывается, т.е. инженерам строителям, мостовикам или инженерам в области ЦОС. Понятно, что не последним.
Подскажите ссылку в интренете, а то об авторе ничего не известно, кем он был, где работал, какие у достижения были, есть.

Для практиков в конце книги есть глава "Применение рядов Фурье к теории изгиба Балок".
Я бы такие книги уничтожал, они для бедных, можно даже сказать сумасшедших людей. Их мало кто понимает (и проверить их крайне сложно), сами они не могут, используя свои знания заработать денег, но будут учить "как надо и как одобрено".
Берем, подставляем, отсюда следует, затем это нужно подставить в [2]...и на деревянных счетах посчитать (вряд ли у него в середине 70-х был доступ к ЭВМ, тогда бы он программу какую привел для расчета...)

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

Формула Эйлера:
https://ru.wikipedia.org/wiki/%D0%A4%D0%BE%...%B5%D1%80%D0%B0

А вот и комплексный сигнал на экране прибора наблюдают:
http://electronix.ru/forum/index.php?showt...p;#entry1379445

Я понимаю, что вы перешли на язык абстракций.
Куда думаю смотреть, а вы "про вообще".
И что вы называете "нулевой частотой"? Постоянную составляющую или основную гармонику?

Я же могу создать идеальный сигнал: синус, амплитудой 100 вольт с фазовым сдвигом 30 градусов, первая гармоника.
Затем с помощью 64-х точечного ДПФ вычислить реальную и мнимую части, модуль и фазовый сдвиг для первой гармоники (могу аналогично для 2...63 гармоник).
Т.е. результат ДПФ совпадает с исходным сигналом.
Где у меня там отрицательная частота? И зачем вы мне формулу Эйлера напомнили?

Вот мой код для расчета ДПФ (IDE NetBeans компилятор MinGW):


Аналогичный результат расчета будет для выборки с 63-й гармоникой, но это уже эффект наложения спектра. От этого мы избавляемся с помощью антиалисинговых фильтров (в расчетах не учавствуют). Т.е. я буду учитывать только 0, 1...31 гармоники (до половины спектра).

Т.е. я что задал, то и получил. Можно сказать "отфильтровал". Что там еще можно отфильтровать? ведь был задан только один сигнал, представьте, что это сигнал с выхода датчика напряжения.

PS: Или вы как-то по другому обсчитываете сигнал и то что я вообще не учитываю (вернее учитываю, но коэффициентом 2 при расчете модуля) из-за симметрии спектра относительно половины частоты дискретизации вы переносите в область отрицательных частот? Мне не понятно.

PPS: Очень хотелось бы не на словах, а на примере Си-шного кода, аналогичного моему увидеть ваши расчеты (создаем один сигнал, считаем и получаем). Желательно, чтобы сигнал был тот же самый: первая гармоника 100 вольт, 30 градусов.

Нулевая частота = постоянка. Просто посчитайте FFT от вашего идеального синуса, постройте график, убедитесь что в спектре 2 гармошки симметричных относительно постоянки.




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

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


Можно и на ассемблере для Cortex M4F, там модуль FPU хороший. Есть и покруче микроконтроллеры, там вообще матричные вычисления поддерживаются, но я не думаю, что мне сейчас это нужно. К тому же совремменные компиляторы с их уровнями оптимизации сделают лучше чем я.
Я понял, что вы хотите мне результат матлаба представить?

PS: Когда есть понимание, то и реализация есть. Я за 10 минут свой готовый код для вас подправил. Хотел ваш увидеть. Не можете на Си, покажите на любом другом языке (кроме матлабовского, хочу максимум понимания).

Нажмите для просмотра прикрепленного файлаЯкобсен-3. Максимальная ошибка измерения 45-55 Гц. ЧД=3200 дпф-512. Шум 0.1%
Красный - дополнительное подавление отрицательных частот.

То же, дпф-1024Нажмите для просмотра прикрепленного файла

Нажмите для просмотра прикрепленного файла дпф-1024 Шум 10%

Код сишный? Я бы хотел проверить его в работе (программно и/или в железе, лучше в железе). Варианта два: в виде EXE-шника или в виде *.lib библиотеки для STM32F407, это если хотите сохранить код закрытым.

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



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

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

PS: Из библиотечных сишных функций я использовал только sinf(), cosf(), sqrtf(), atanf(), т.е. синус, косинус, корень квадратный и арктангенс для чисел типа float.


Предубеждения у меня нет, просто чтобы мне посмотреть ваш пример, мне нужен такой же симулинк как у вас.
Лучше бы было без симулинка в коде, чтобы этот код можно было без симулинка использовать и чтобы этот код был оптимизирован для ARM MCU.

Матлаб можно на рутрекере взять посмотреть.
http://electronix.ru/forum/index.php?s=&am...t&p=1380116

Спасибо за труд, за желание помочь. Воздастся.
Для этого файла из архива похоже нужна последняя версия. Я её так и не смог установить. Попробовал открыть седьмым, не открывает.
С этими коммерескими продуктами всегда какая-нибудь бяка вылезает, то версия не та, то обновил версию, что-нибудь из старого не так работает.

Я уже и про комплексный синус узнал из источника, вот скриншот одной из страниц, никогда нужно не было. Там указано, что если комплексная часть равна нулю, то комплесный синус превращается в обычный.
Интересен следующий момент, что данные с выхода АЦП не являются комплексными (например 16-ти разрядный АЦП) и при использовании комплексного БПФ комплексная выборка представлена только Re частью, а Im часть приравнивается нулю. Т.е. N выборок АЦП представляются N комплексными числами, у которых действительная часть отлична от нуля, а мнимая равна нулю. Так работает CMSIS DSP библиотека для Cortex M4F MCU.
Что-нибудь можете прояснить, что это и зачем? Т.е. я не понимаю, существует ли комплексный синус и используется ли он в практических расчетах вне матлаба.

Это когда Баскакова обратно склеите. В этом топике информации для реализации более, чем достаточно.

Что-то Вы залезли в такие дебри математики... Вот результаты моделирования. 3200 Гц дискретизации 1024 точки Фурье, окно Гаусса. При отсутствии шума во входном сигнале точность лучше 6 знаков после запятой. Далее привожу файлы результатов при 5% белого шума и 10%.
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла

Странное у вас условие. Если есть знания, то их нужно реализовывать. На память приходит товарищ Перельман и доказательство гипотезы Пуанкаре.
Странные люди, эти математики (с хорошим знанием математики, если оно есть, если это не пользователи какого-либо матпакета), нет бы все на благо Родины, а они "я знаю, но не скажу...мне ничего не нужно, я особенный, а вы балбесы еще лет 30 нничего сделать не сможете...". . Григорий, если это вы, то так нельзя, не допустимо и не позволительно, потомки вам этого не простят.

]
Результат хороший. Каким методом считали? Эрик Якобсен, или по Гауссу но как-то по своему...?
У вас в обоих файлах входной сигнал задается через каждые 0,10 Гц, а хотелось бы увидеть через 0,01 Гц. Требуется же разрешение показометра 0,01 Гц и ошибку для такой сетки частот. Если *.txt файл будет большой, то открою с помощью notepad++.

Вот уже несколько пользователей писали про сигнал с подмешанным шумом в процентах от сигнала. Вы создаете такое сигнал, наверное, в матлабе, а вот как мне оценить уровень шума в моей системе? Я могу лишь на частоте дискретизации взять выборки АЦП. Вот что получится. Какой здесь уровень шума, кто-нибудь может ответить.

Все сигналы отображаются в отсчетах АЦП. При отсутствии сигнала, например, для 16-ти разрядного АЦП (-32768...32767) в канале Ic амплитудное значение с выхода АЦП равно от -4 до 1 (это шум системы). Сколько процентов равен шум в моем случае? Какое отношение сигнал/шум возможно при амплитуде сигнала 1000 разрядов АЦП и этом шуме?

Мне кажется, что тема довольно неплохо была раскрыта вот тут:
http://electronix.ru/forum/index.php?showtopic=84446
По крайней мере, можно использовать как стартовую точку.

Ну и моделирование. Матлаб- "наше все".

Тогда, учитывая, что вы "там" увидели раскрытие темы и пока не забыли, подсобите простым крестьянам, где там "ключ" и как им воспользоваться? И еще интересует, где начинается "стартовая точка", т.е. откуда начинать. На одном языке разговариваем, а дело сделать не можем. Кто-то утверждает, что уже все сделал, но в матлабе. Ну не генерит же матлаб Си-шный код для Cortex M4F. А без матлаба мы никак, во какая зависимость.

PS: Так мы до профтехучилищ дойдем: специальность "оператор ЭВМ", специализация "пользователь матлаб".

Почему это не генерирует? Еще как. И не один год уже.

Ну, мне "там" помогли. В процессе "раскрытия" я освоил Матлаб, изучил его язык, сделал матмодель, проверил ее на реальных данных (залитых в компьютер с помощью того же Матлаба), соорудил си-программу (FFT использовал готовый, но без закрытых либ- все на чистом стандартном си), протестировал в С++Билдере, перевел на мое тогдашнее железо (PIC32), проверил в симуляторе на подсунутых известных данных, проверил в отладчике от генератора, проверил от реального датчика. Уф! Шаг за шагом.

Коротко по методу: FFT (у меня-4096 точек), окно(у меня-Гаусс), интерполяция (у меня-квадратичная интерполяция по трём точкам).
Но чудес не бывает, нужно хоть немного понимать что делаете и что от чего и как зависит, иначе результат случаен. А для этого нужно немного теорию почитать и много моделировать (в матлабе).

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

Вы тут уже утомили всех своими претензиями..

Вот этот код в MATLAB'e дает ошибку меньше чем 10^-3 [Hz]:


Ошибка измерения:
Нажмите для просмотра прикрепленного файла

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

Спасибо!


Какие могут быть претензии, это от беспомощности.

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

Это как раз то, что мне нужно. Совсем другое дело! Остается перенести на Си и добиться аналогичного результата! Спасибо!
Вообще, лучшие импортные приборы измеряют частоту основной гармоники как раз до 0,001 Гц, показометр так у них показывает.



Круто! Этот код в седьмом работает и такую же картинку выдал (ошибку). Значит это готовая модель.

У меня программа на С, сигнал на входе 1, я подмешиваю просто RND, деленный на фиксированный коэффициент. В примерах выше на 20 и 10, соответственно. Это означает, что максимальное отклонение будет +-0.05 или +-0.025. Среднеквадратичное меньше, какое именно - считать надо.
Для шага 0.01 посчитаю вечером, кода домой доберусь. При отсутствии шума все будет хорошо, там нули до 6 знака после запятой. Какой шум Вам добавить, чтобы было похоже на Ваш реальный? И я еще попробую ввести Вашу дискретность на 1024 отсчета, посмотрим как отразится на результате. Сейчас-то синус на входе точный.

Если получится, то попробуйте добавить шум +/- 2 разряда 16-ти разрядного АЦП (диапазон преобразования -32768...+32767), как на скриншоте сигнал Ic из этого поста. И не стесняйтесь показывать свои знания.

Проблем с добавлением моей дискретности у вас не должно возникнуть, просто берете 1024 точки через интервал 1/3200 секунды, или вообще абстрагируетесь от времени.

Вот и результаты. Здесь амплитуда сигнала +-32768 и точность входного сигнала ограничена целым. Плюс шум со среднеквадратичным значением 1.7:
Нажмите для просмотра прикрепленного файла
Здесь амплитуда +-1000 плюс тот же шум:
Нажмите для просмотра прикрепленного файла
Даже при небольшой амплитуде точность расчетов превосходит требуемую на полтора порядка.

В моем возрасте я уже абсолютно точно знаю, что умею хорошо, что плохо и в чем мне вломину разбираться.

И что, дальше никому не интересно? Написали бы, чем дело кончилось.

Здравствуйте!
Вопрос очень похожий, задам тут же.

Изменились частоты и семплы, но цель такая же- как можно точнее считать частоту (точнее, чем бины FFT)
что не так как раньше:
1) Исходные данные поступают в 5 (или 10) раз чаще, чем нужно для FFT
2) нужно применить окно Ханна: 0.5-0.5*cos(2pi*n/N)

вопрос: при таком окне можно ли применить интерполяцию для уточнения частоты по вычисленному FFT? и какую?
До этого при использовании окна Гаусса применял квадратичную интерполяцию, и все отлично было.
Можно ли с Ханном ту же квадратичную интерполяцию применить?

И влияет ли метод децимации на качество результата финальной интерполяции в области максимального бина?
Что лучше применить для децимации с фактором 5 или 10? простое среднее, IIR, FIR ?


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

Форма спектра будет совпадать с видом окна. Для хорошей точности нужна аппроксимация соответствующей функцией. Так что квадратичная в этом случае не подойдет. Здесь нужен косинус. Можно попробовать полином (первых три члена от разложения в ряд), но нужно меделировать и смотреть какая точность получится. И я не помню наизусть каких степеней там ряд.
С децимацией фильтр, безусловно лучше. Но еще нужно посмотреть, даст ли он выигрыш в счете по сравнению с Фурье на весь массив сразу. Последний вариант должен быть лучше по точности и усреденения и фильтров.

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

Про FFT без децимации сразу на весь массив- сам бы так хотел. Но это у меня просто еще одна дополнительная опция в приборе, желательно бы обойтись без сотрясения основ уже сделанного. А CMSIS-овская уже используемая библиотека только FFT4096 умеет, вот и приходиться децимировать.
Мне бы вообще-то очень бы подошло что-то типа FFT длиной 16К или даже 32К точек float32, но это уже свой код применять вместо встроенных библиотек, про это тоже думал. Снижение быстродействия даже в 5 раз еще приемлемо, по сравнению с CMSIS. Тоже интересно подумать.

Попробуйте, я не смотрел, как это будет оптимизироваться на АРМ'е, но на DSP ложится прилично.
Нажмите для просмотра прикрепленного файла

Вы будете таки смеяться, но именно этот код я использую уже лет 5, про него и говорил.
Работает с небольшими вариациями в разных системах, начиная от мелкого 8-битного Майкрочипа, где данные храню во внешней микросхеме сериальной RAM. Но именно в АРМе я его тоже не пробовал.
Причем, действительно, работает всегда и везде. Нашел в Интернете.
Ну раз и Вы его тоже рекомендуете, значит он действительно неплох, не буду искать что-то еще.
Спасибо!

Интерполяцию окном Hanna смотри http://calculator2006.narod.ru/articles/haan.htm

отлично, спасибо!!!!
Скопировал себе на диск и пдф-ку сделал, а то в интернете бывает так- сегодня есть, а завтра -"404"

То есть расчет идет по трем точкам.
И странно там в тексте как-то окно задано, формула 1-cos(2pi*n/N) правильная? Она же двойку дает в максимуме, даже с рисунком в том же тексте не совпадает.
Думаю, просто опечатка, должно быть все-таки 0.5*(1-cos(2pi*n/N)), собственно так окно в первоисточнике и задается.

Будет спектр чуть выше.
И да скорей всего опечатка я у же и не помню почему так вышло.

Боюсь, что для Вашей жуткой точности не хватит интнрполяции по трем точкам. Лучше бы 5-7 точек или даже больше.