Comments 20
Такая штука на дифурах неустойчива – амплитуда или падает в ноль, или улетает в бесконечность. Надо или калиброваться (например, в момент перехода "синуса" через ноль задавать косинусу фиксированное значение), или сделать контур с определённой добротностью и возбуждать его периодическим (да хоть прямоугольным) сигналом.
Кстати, удобнее думать не о синусе и косинусе, а о комплексной экспоненте. Тогда и шаг по времени – просто умножение на константу, и калибровка очевидна.
Ну и если вам нужна свёртка с синусом-косинусом – возможно, удобнее будет просто подавать исходный сигнал на тот самый фильтр, который вы используете для генерации синуса и косинуса.
Именно это я и делаю). Получилось устойчиво. Впрочем не проверял на очень длинных выборках. Но свертка работает прекрасно. Что касается устойивости, да, возможно ошибки при вычислении конечных разностей приведут к накоплению лишних значений, что эквивалентно энергии. То же самое будет с интегральной ошибкой в любой системе, но если внести небольшую поправку в добротность, то колебания будут очень медленно затухать. Но период синхронизации можно выбрать очень большим и он может быть любым. В обычном случае значения синуса придется интерполировать или просчитывать для каждого отсчета.
В цифровой форме это будет приближенно. Интересно узнать как аппроксимировали дробно-рациональной функцией передаточную функцию КК. Там в z-пространстве иррациональная функция.
Про Cordic снова напишу. Генерирует сразу и синус и косинус.
CORDIC хорош для произвольных углов, а у автора последовательный перебор, тут можно проще, главное погрешность слишком большую не накопить.
Вот именно, не накопить. Может, лучше иметь таблицу со значениями и от одного-двух ближайших уже рассчитывать через ряды? (не простая апроксимация).
У меня есть расчет sin(x), cos(x) на базе ряда Тейлора с точностью около 1LSB 32-битного числа. Используется математика в целых числах на Arm Cortex M. Шесть умножений кажется 32bit x 32bit => 64bit
Понимаю, что таблица + линейная интерполяция может оказаться немного быстрее, однако не уверен что это даст 1LSB.
В программах БПФ, ЦФ используется табличный способ и тригонометрические формулы Sin, Cos суммы(разности) углов. Всегда использую эти способы в программах цифровой обработки данных особенно на микроконтроллерах.
Какое преимущество имеет Ваш вариант вычисления ?
Есть такая штука как Numerically Controlled Oscillator, NCO. И там есть момент, связанный со спектральной чистотой генерируемого колебания.
NCO, на сколько я знаю, считает sin/cos по таблице с интерполяцией, или размыванием спектра. Первый способ точнее, и оптимален по соотношению затрат мощности к точности в ущерб площади кристалла. Тем более, что sin и cos - производные друг-друга. Второй способ проще, но менее точный.
TI GC4016 datasheet например
С размыванием спектра не понял. А так, да, по таблице.
Наверное, я неправильно написал - всё-таки оба способа применяют одновременно. LUT тут и так немалой разрядности, или разрядность под капотом всё-таки увеличивают с помощью интерполяции.
Но вот, чтобы по-максимуму уйти от нелинейности ЦАП, дополнительно водят широкополосный фазовый шум. Шума в спектре становится больше, но, как TI считает, не критично. А вот высота палок гармоник и комбинационных частот уменьшается. Как-то так.
Я думаю там просто LUT. А дизеринг вводят из-за недостатков аккумулятора фазы. Я это говорю, ориентируясь на Matlab. Там никакого ЦАП нет. Блок NCO самодостаточный.
Я думаю там просто LUT
Возможно. Я не готов утверждать, но почти мегабайтный LUT может занимать места на кристалле больше, чем схема линейной интерполяции. А может, и тупой большой LUT будет жрать меньше энергии. В эквивалентной схеме TI могут нарисовать всё, что угодно, но как оно реализовано под капотом - не знаю.
По сути тут техасские товарищи просто улучшают SFDR, не увеличивая разрядности данных. Для узкополосной фильтрации при обработке сигнала SFDR часто важнее, чем SNR. И даже не важно, генерируется ли физически это сигнал через ЦАП, или просто используется в расчётах.
Догадываюсь о чем идет речь. У меня есть FPGA от GoWin и Altera. Как-то, в порядке эксперимента, хотел синтезировать "честные" квадратуры минимальными усилиями и недорого. Использовал четыре резистора по два на выход. Побоялся разброса, по этому использовал не R - 2R, а просто пару одинаковых резисторов. Получил три логических уровня 0, 1, 2. И таких два канала.
Купил в радиодеталях пару недорогих кварцевых фильтра на промежуточную частоту 455КГц. На FPGA, использовал сдвиговые регистры как самый быстрый элемент на LUT + триггер (планировалось использовать 10.7MHz потом). Раскачивал кварцевые фильтры на резонансе фактически в двух точках + "0". Хотелось использовать их для сглаживания фазы, но не как полосовой фильтр для прямоугольника а на резонансе.
Далее, по фигурам Лиссажу, смотрел на окружность, хотя скорее, все-же, немного эллипс, который образован двумя одинаковыми каналами.
Далее, хотел использовать два двойных балансных смесителя для получения квадратур. Этакая AD9361 с на фиксированной частоте и "на рассыпухе". Стандартный "прямоугольник" не хотелось из-за энергии нечетных гармоник. Это работает, но конечно, хуже чем цифровая копия. Макет рядом с резонатором трогаешь и фаза плывет превращая почти окружность в эллипс.
Но я системный программист-микроконтроллерщик / BSP-шник и реверс-инженер. ЦОС от меня далековат, впрочем кое что интересно. Для хобби нужны деньги а в РФ сейчас по этому направлению не очень. Или мне так невезет просто.
Да, похоже. С 1,5-битным ЦАП вам как раз помогло бы внесение доли цифрового шума. И аккуратный выбор частоты дискретизации, чтобы паразитные палки в спектре не попадали в полосу фильтра. Применяют же 1-битный АЦП при приёме сигналов с SNR<1, например, GNSS. Так что и обратное вполне реально.
А вот с кварцевым фильтром всё очень просто: у него гарантируется полоса, и линейность фазы. Но вот абсолютное значение фазы, или задержки - не гарантируется. И по факту оно случайное, и зависит от погоды на Марсе.
Один фильтр сам по себе - работает. Два фильтра вместе - нет. Поэтому для генерации точных квадратур приходится полосу пост-ЦАП-фильтра сильно расширять.
Ну алгоритм довольно нестабильный на реальных разрядностях
Вот так выглядят амплитуды при разных способах округления в младших разрядах
Программа на Python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Drift caused by biased rounding (Q1.15 fixed-point oscillator)
"""
import math
import numpy as np
import matplotlib.pyplot as plt
# ------------------------------------------------------------
# ПАРАМЕТРЫ ЭКСПЕРИМЕНТА
# ------------------------------------------------------------
Fs = 1_000.0 # частота дискретизации, Гц
f0 = 5.0 # генерируемая частота, Гц
seconds = 60 # длительность теста, с
N = int(Fs * seconds)
omega = 2 * math.pi * f0 / Fs
kappa_f = 2 * math.sin(omega / 2) # float-коэффициент
Q = 15 # формат Q1.15
S = 1 << Q # масштаб 2^Q
kappa = int(round(kappa_f * S)) # фикс-пойнт коэффициент
# ------------------------------------------------------------
# ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ
# ------------------------------------------------------------
def one_step(c: int, s: int, rounding: str = 'trunc'):
"""
Один рекурсивный шаг генератора:
c_{n+1} = c_n - kappa*s_n
s_{n+1} = s_n + kappa*c_{n+1}
rounding:
'trunc' – обычный сдвиг вправо (предвзятое округление)
'sym' – симметричное округление к ближайшему
"""
if rounding == 'trunc':
prod1 = (kappa * s) >> Q
else: # symmetric
tmp = kappa * s
prod1 = (tmp + (1 << (Q-1))) >> Q if tmp >= 0 else (tmp - (1 << (Q-1))) >> Q
c_next = c - prod1
if rounding == 'trunc':
prod2 = (kappa * c_next) >> Q
else:
tmp = kappa * c_next
prod2 = (tmp + (1 << (Q-1))) >> Q if tmp >= 0 else (tmp - (1 << (Q-1))) >> Q
s_next = s + prod2
return c_next, s_next
def run(rounding: str):
"""Запускаем генератор и возвращаем массив амплитуд"""
c, s = S, 0 # cos(0)=1, sin(0)=0 в Q1.15
amp = np.empty(N)
for n in range(N):
amp[n] = math.hypot(c / S, s / S) # мгновенная амплитуда
c, s = one_step(c, s, rounding)
return amp
# ------------------------------------------------------------
# МОДЕЛИРУЕМ ОБЕ СХЕМЫ И РИСУЕМ
# ------------------------------------------------------------
amp_trunc = run('trunc')
amp_round = run('sym')
t = np.arange(N) / Fs
plt.figure(figsize=(10, 6))
plt.plot(t, amp_trunc, label='truncate-shift', lw=0.8, color='orange')
plt.plot(t, amp_round, label='symmetric-round', lw=0.8, color='orangered')
plt.title('Drift caused by biased rounding (Q1.15)')
plt.xlabel('Time [s]')
plt.ylabel('Instant amplitude')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
Если нужна квадратурная детекция частот, то напрямую использовть Goertzel гораздо практичнее.
На самом деле все интереснее.
Вообще, я делитант в ЦОС. Но кое-какие остаточные знания имеются. Еще на втором курсе университета, меня попросили в качестве кусовой работы сделать цифровой однополосный модулятор. Предлагали использовать преобразование Гильберта чтобы получить квадратуры сигналов.
Что такое левая/правая боковые полосы? По сути это же вращающиеся в разные стороны вектора с меняющимися начальными фазами. Я предположил следующее. Если разбить сигнал на его спектральные составляющие и каждую из составляющих домножить на соответствующий вектор вращения (cos(Wn*t) - i*sin(Wn*t)), то получим I/Q отсчеты для квадратурного модулятора для данной составляющие. Берем N таких составляющих. А затем складываем получившиеся вектора. По идее это просто суммирование их координат. В итоге имеем координату вектора для модулятора.
О чем я подумал? Колебательный контур имеет две реактивных составляющих и в идеальных условиях на выходе дает синусообразную составляющую. Производная от этой составляющей есть та же функция, но сдвинутая на некоторый угол. В идеале,- PI/2.
Если из модели колебательного контура получить обе составляющих, то получаем сразу и фильтр и вращающийся вектор с нужной нам частотой. Таким образом остается сложить составляющие в от разных контуров и получить квадратуры для однополосного модулятора.
А из математики, там лишь пара умножений. Т.е. если семплировать речь с частотой 128КГц, получается порядка 32семплов на отсчет. Это уже достаточно для того чтобы использовать "решение уравнения" для получения квадратур.
Теоретически, можно использовать в частотном уплотнении каналов. Правда не знаю, использует ли это сейчас кто-то кроме военных, ведь проще ортогональные сигналы синтезировать, а передавать в цифре.
Кроме того на вход дифференциального уравнения можно подать сигнал, который будет раскачивать его как реальный колебательный контур
Ну вообще для каких-то применений - очень неплохой способ, плюсую. БИХ-фильтром смоделировать колебательный контур и заявить, что частота вынужденных колебаний равна частоте вынуждающей силы. И устойчивость в этом случае как раз гарантирована.
Главное, не забыть, что разрядность вычислений узкополосного БИХ-фильтра должна быть не маленькая.
Метод синтеза синусоидальных колебаний, с применением цифрового колебательного контура