Формирование околосинусоидальных сигналов

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

Если есть возможность вызвать встроенную функцию и получить результат, то все довольно тривиально и не вызывает затруднений. Однако, если это требуется например на какой-либо микроархитектуре, то задача сильно усложняется. 
В этой заметке я предлагаю поближе взглянуть на эту всем известную функцию и обсудить различные способы ее формирования. Кроме того, давайте условимся заранее, что автор не самый лучший рассказчик, да и ошибок может наделать. Поэтому не стесняемся и пишем свое мнение в комментариях.

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

Теперь предлагаю рассмотреть основные способы формирования наших сигналов. Первый и самый очевидный - табличный. Берем MATLAB, ну или то, что есть под рукой. Задаем количество точек и получаем четверть периода. 

angle = 0: (pi/2)/90 : pi/2; %brackets to show quater of periodogram
values = sin(angle);
figure('Name','Quater of sine');
plot(values);

Количество точек будет определяться не только необходимой частотой дискретизации, требованием к шумам квантования и т.д., но и объемом памяти, который мы хотим занять относительно бесполезными данными, служащими только для одной цели. Если вспомнить, что я говорил про микроархитектуры в начале заметки, то может выясниться, что память, занимаемая массивом, гораздо больше исполняемого кода программы...

Кроме всего прочего, будет здорово вывести выражение для определения выходной частоты такого сигнала:

Ft = 1000; %Hz
Foutput = Ft/(length(angle) - 1)/4;

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

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

x = linspace(0,4*pi,10);
y = sin(x);

p = polyfit(x,y,7);

x1 = linspace(0,4*pi);
y1 = polyval(p,x1);
figure
plot(x,y,'o')
hold on
plot(x1,y1)
hold off

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

Вычислять степенные функции задача довольно трудоемкая, а если еще предположить, что у нас используется фиксированная точка, то лучшим выходом будет написать небольшую модель, которая позволит отследить переполнения и рассчитать СКО при различных степенях полинома. Давайте для примера понизим степени полинома хотя бы до 5. 

Очень красивый синус вышел. Нельзя не отметить, что в MATLAB есть хороший и полезный инструмент для анализа поведения и портирования алгоритмов в фиксированную точку. Видел презентации, нравится, но только он предназначен для профессионального использования и в хобийный бюджет даже при большом желании не вписывается. 

Теперь предлагаю обсудить то, к чему мы так долго шли, третий способ формирования сигналов, очень похожих на синусоидальные. Фамилию математика я все-таки рискну руссифицировать, но для порядка приведу запись на латинице - "Bhaskara I's sine approximation formula", в общем будем называть предложенное решение аппроксимацией Бхаскары 1го. Да, был еще и второй известный математик из Индии с такой же фамилией. Итак, сама формула:

Здесь аргумент - это угол в градусах. Предлагаю построить и посмотреть, как она выглядит. 

rad = 0:pi/180:pi/2;
x = 0:1:90;

approx = 4.*x.*(180.-x)./(40500.-x.*(180.-x));
sine = sin(rad);

plot(sine, '*')
hold on
plot(approx, 'x')

На картинке вряд ли что-то будет понятно, точки очень близко, когд есть так, что можно без труда посмотреть на своем пк в удобном масштабе. Отмечу, что аппроксимация работает только до 180 градусов, но как я ранее упоминал, нам и четвертинки периода достаточно.

Перед тем как посмотреть еще одну вариацию данной формулы, построим график ошибки.

rad = 0:pi/180:pi/2;
x = 0:1:90;

approx = 4.*x.*(180.-x)./(40500.-x.*(180.-x));
sine = sin(rad);

diff = sine.-approx;

plot(diff, '*')

Ошибка лежит в пределах +- 0.0015, соответственно при переезде в фиксированную точку, например в Q12, она прячется за шумами квантования. Предлагаю взглянуть на еще более упрощенную аппроксимацию, правда для ее использования можно применять только целочисленные значения угла в градусах. 

Хочу еще сказать пару слов про переход к фиксированной точке. В целом, все более-менее тривиально: во-первых, выбираем где она будет находиться. Затем пишем функцию преобразования констант и входных аргументов с детекстирование переполнений и насыщением. Затем разбиваем формулу на элементарные шаги и устанавливаем триггеры переполнения. Далее пишем тестовое окружение, гоним массив данных и считаем ошибку, максимальную или ско, как больше нравится. Думаю, однажды я соберусь с силами и мы этот путь пройдем вместе. А пока всем успехов и до скорого.