Добрый вечер, уважаемые коллеги. Статья будет короткой и возможно натолкнет кого-то на реализацию. У меня получилось.
Недавно разрабатывал код, который рассчитывает значения sin(x), cos(x) и квадратного корня из x, на arm. По сути это была апроксимация рядом Тейлора. Но код написан на Assembler и выполнялся относительно быстро. Изначально предполагалось использовать его для своего станка. Немного позже задался вопросом одновременной генерации sin(x) и cos(x) как сигналов для свертки с исходным сигналом. Некоторые, предлагали cordic, но я пошел другим путем.
Изначально, моделировал систему дифференциальных уравнений, использующих два умножения и конечные разности для получения sin(x) / cos(x) на каждый временной отсчет. Фактически, это отклик дифференциального уравнения на дельта-функцию (единичный импульс) который запускает колебательный процесс. Результат вычисления и есть "синусойды". Они получились не ортогональными при небольшом количестве отсчетов на период, из-за неточности вычисления сдвига фазы. Однако, при использовании длительных последовательностей данный метод может заменить классические методы генерации sin(x)/cos(x) поскольку конечные разности становятся небольшими.
Кроме того на вход дифференциального уравнения можно подать сигнал, который будет раскачивать его как реальный колебательный контур. Например, у Вас возникнет желание свернуть сигнал с длинноволновыми последовательностями. В этом случае, можно, просчитывать каждый отсчет sin(x)/cos(x), а можно просто подать входной сигнал с АЦП на вход такого дифференциального уравнения и получить "синхронный интегратор", настроенный на определенную частоту.
Поскольку добротность подобного контура может быть любой, все зависит от распределения шума. Сколько кореляторов можно реализовать для СДВ? Сколько коэффициентов требуется для КИХ/БИХ фильтра? В моем случае, требуется две производных (конечных разности).
Кроме того, используя фильтр основанный на моделе колебательного контура, в моем случае, происходит сглаживание PSK-модулированного сигнала, который формируется как синусойда со сдвигом фазы, домноженная фильтр приподнятого косинуса, с перекрытием в 1/2.
Комментарии (14)
sci_nov
27.06.2025 19:13В цифровой форме это будет приближенно. Интересно узнать как аппроксимировали дробно-рациональной функцией передаточную функцию КК. Там в z-пространстве иррациональная функция.
mozg37
27.06.2025 19:13Про Cordic снова напишу. Генерирует сразу и синус и косинус.
aamonster
27.06.2025 19:13CORDIC хорош для произвольных углов, а у автора последовательный перебор, тут можно проще, главное погрешность слишком большую не накопить.
dyadyaSerezha
27.06.2025 19:13Вот именно, не накопить. Может, лучше иметь таблицу со значениями и от одного-двух ближайших уже рассчитывать через ряды? (не простая апроксимация).
EvgenySbl Автор
27.06.2025 19:13У меня есть расчет sin(x), cos(x) на базе ряда Тейлора с точностью около 1LSB 32-битного числа. Используется математика в целых числах на Arm Cortex M. Шесть умножений кажется 32bit x 32bit => 64bit
Понимаю, что таблица + линейная интерполяция может оказаться немного быстрее, однако не уверен что это даст 1LSB.
dyadyaSerezha
27.06.2025 19:13Я же специально написал, что не простая апроксимация/интерполяция. Возможно, можно использовать тот же ряд Тейлора, но уже с в 2 раза меньшим кол-во операций. Ну и что-то мне подсказывает, что в инете должна быть куча вариантов таких расчётов для разных условий.
nikolz
27.06.2025 19:13В программах БПФ, ЦФ используется табличный способ и тригонометрические формулы Sin, Cos суммы(разности) углов. Всегда использую эти способы в программах цифровой обработки данных особенно на микроконтроллерах.
Какое преимущество имеет Ваш вариант вычисления ?
sci_nov
27.06.2025 19:13Есть такая штука как Numerically Controlled Oscillator, NCO. И там есть момент, связанный со спектральной чистотой генерируемого колебания.
alcotel
27.06.2025 19:13NCO, на сколько я знаю, считает sin/cos по таблице с интерполяцией, или размыванием спектра. Первый способ точнее, и оптимален по соотношению затрат мощности к точности в ущерб площади кристалла. Тем более, что sin и cos - производные друг-друга. Второй способ проще, но менее точный.
TI GC4016 datasheet например
Indemsys
27.06.2025 19:13Ну алгоритм довольно нестабильный на реальных разрядностях
Вот так выглядят амплитуды при разных способах округления в младших разрядах
Программа на 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 гораздо практичнее.
EvgenySbl Автор
27.06.2025 19:13На самом деле все интереснее.
Вообще, я делитант в ЦОС. Но кое-какие остаточные знания имеются. Еще на втором курсе университета, меня попросили в качестве кусовой работы сделать цифровой однополосный модулятор. Предлагали использовать преобразование Гильберта чтобы получить квадратуры сигналов.
Что такое левая/правая боковые полосы? По сути это же вращающиеся в разные стороны вектора с меняющимися начальными фазами. Я предположил следующее. Если разбить сигнал на его спектральные составляющие и каждую из составляющих домножить на соответствующий вектор вращения (cos(Wn*t) - i*sin(Wn*t)), то получим I/Q отсчеты для квадратурного модулятора для данной составляющие. Берем N таких составляющих. А затем складываем получившиеся вектора. По идее это просто суммирование их координат. В итоге имеем координату вектора для модулятора.
О чем я подумал? Колебательный контур имеет две реактивных составляющих и в идеальных условиях на выходе дает синусообразную составляющую. Производная от этой составляющей есть та же функция, но сдвинутая на некоторый угол. В идеале,- PI/2.
Если из модели колебательного контура получить обе составляющих, то получаем сразу и фильтр и вращающийся вектор с нужной нам частотой. Таким образом остается сложить составляющие в от разных контуров и получить квадратуры для однополосного модулятора.
А из математики, там лишь пара умножений. Т.е. если семплировать речь с частотой 128КГц, получается порядка 32семплов на отсчет. Это уже достаточно для того чтобы использовать "решение уравнения" для получения квадратур.
Теоретически, можно использовать в частотном уплотнении каналов. Правда не знаю, использует ли это сейчас кто-то кроме военных, ведь проще ортогональные сигналы синтезировать, а передавать в цифре.
alcotel
27.06.2025 19:13Кроме того на вход дифференциального уравнения можно подать сигнал, который будет раскачивать его как реальный колебательный контур
Ну вообще для каких-то применений - очень неплохой способ, плюсую. БИХ-фильтром смоделировать колебательный контур и заявить, что частота вынужденных колебаний равна частоте вынуждающей силы. И устойчивость в этом случае как раз гарантирована.
Главное, не забыть, что разрядность вычислений узкополосного БИХ-фильтра должна быть не маленькая.
aamonster
Такая штука на дифурах неустойчива – амплитуда или падает в ноль, или улетает в бесконечность. Надо или калиброваться (например, в момент перехода "синуса" через ноль задавать косинусу фиксированное значение), или сделать контур с определённой добротностью и возбуждать его периодическим (да хоть прямоугольным) сигналом.
Кстати, удобнее думать не о синусе и косинусе, а о комплексной экспоненте. Тогда и шаг по времени – просто умножение на константу, и калибровка очевидна.
Ну и если вам нужна свёртка с синусом-косинусом – возможно, удобнее будет просто подавать исходный сигнал на тот самый фильтр, который вы используете для генерации синуса и косинуса.
EvgenySbl Автор
Именно это я и делаю). Получилось устойчиво. Впрочем не проверял на очень длинных выборках. Но свертка работает прекрасно. Что касается устойивости, да, возможно ошибки при вычислении конечных разностей приведут к накоплению лишних значений, что эквивалентно энергии. То же самое будет с интегральной ошибкой в любой системе, но если внести небольшую поправку в добротность, то колебания будут очень медленно затухать. Но период синхронизации можно выбрать очень большим и он может быть любым. В обычном случае значения синуса придется интерполировать или просчитывать для каждого отсчета.