Комментарии 15
Такая штука на дифурах неустойчива – амплитуда или падает в ноль, или улетает в бесконечность. Надо или калиброваться (например, в момент перехода "синуса" через ноль задавать косинусу фиксированное значение), или сделать контур с определённой добротностью и возбуждать его периодическим (да хоть прямоугольным) сигналом.
Кстати, удобнее думать не о синусе и косинусе, а о комплексной экспоненте. Тогда и шаг по времени – просто умножение на константу, и калибровка очевидна.
Ну и если вам нужна свёртка с синусом-косинусом – возможно, удобнее будет просто подавать исходный сигнал на тот самый фильтр, который вы используете для генерации синуса и косинуса.
Именно это я и делаю). Получилось устойчиво. Впрочем не проверял на очень длинных выборках. Но свертка работает прекрасно. Что касается устойивости, да, возможно ошибки при вычислении конечных разностей приведут к накоплению лишних значений, что эквивалентно энергии. То же самое будет с интегральной ошибкой в любой системе, но если внести небольшую поправку в добротность, то колебания будут очень медленно затухать. Но период синхронизации можно выбрать очень большим и он может быть любым. В обычном случае значения синуса придется интерполировать или просчитывать для каждого отсчета.
В цифровой форме это будет приближенно. Интересно узнать как аппроксимировали дробно-рациональной функцией передаточную функцию КК. Там в 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 например
Ну алгоритм довольно нестабильный на реальных разрядностях
Вот так выглядят амплитуды при разных способах округления в младших разрядах
Программа на 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семплов на отсчет. Это уже достаточно для того чтобы использовать "решение уравнения" для получения квадратур.
Теоретически, можно использовать в частотном уплотнении каналов. Правда не знаю, использует ли это сейчас кто-то кроме военных, ведь проще ортогональные сигналы синтезировать, а передавать в цифре.
Кроме того на вход дифференциального уравнения можно подать сигнал, который будет раскачивать его как реальный колебательный контур
Ну вообще для каких-то применений - очень неплохой способ, плюсую. БИХ-фильтром смоделировать колебательный контур и заявить, что частота вынужденных колебаний равна частоте вынуждающей силы. И устойчивость в этом случае как раз гарантирована.
Главное, не забыть, что разрядность вычислений узкополосного БИХ-фильтра должна быть не маленькая.
Метод синтеза синусоидальных колебаний, используя цифровой колебательный контур