И опять же, этот пост для меня есть не что иное, как крик души. Меня поражает тот факт, что некоторые вещи из повседневной жизни используются нами на веру. Взял формулу, подставил значения и всё работает (пока, вроде бы:). Не нужно так! «Но ведь всё же верно!» — скажете вы. «Дяденька из форума сказал, что всё будет работать!». Возможно «дяденька» и прав, но как говорил классик: «Верить нельзя никому!». Мне же — можно (после перепроверки). Ведь задача инженера не только в создании. Задача инженера также и в глубоком понимании запущенного процесса. Ярким примером того, что «это» будет работать, стала рекурсивная формула генератора синусоидального сигнала, которую я неоднократно встречал на многих сайтах и блогах, посвящённым цифровой обработке сигналов. Естественно, я «испытывал» эту формулу на практике и был действительно удивлен такой простотой. Я был настолько поражен, что не удержался и задался целью понять, откуда берётся данная чудо-формула. Меня не сильно удивило и то, что толкового, ясного и простого способа получения формул в сети я не нашёл. Нашёл множество несуразностей, эмпирико-дедуктивных методов и прочего «мусора» но для себя — НИЧЕГО.

Для знатоков (почему-то читающих эту статью) скажу сразу, что метод получения формул, который базируется на размещении 2 полюсов сразу на Z плоскости и симметрично относительно реальной координатной оси конкретно для меня был абсолютно непонятным методом, с множеством предположений выданных автором и взятых мною опять же на веру. Новички в ЦОС сразу же задают вопрос о том, зачем брать 2 полюса, а не скажем четыре или всего один. Следующим вопросом может стать вопрос о количестве и размещении нулей… Ну да ладно. Это я так, для знатоков!

Ну а если серьезно, то дело до конца мне всё же удалось довести. Результатом делюсь с вами. Что вы получите в результате? Простую рекурсивную формулу, с помощью которой можно генерировать синусоиду, обходясь без вызова тригонометрической функции из библиотеки, которая к слову, будет работать в разы медленнее предложенного способа (а это, скажем, для контроллера даже очень важно). Скажу сразу, что зная два предыдущих значения синусоиды можно однозначно вычислить следующее.

Приступим!

Конструирование для себя я начал с так называемого технического задания ТЗ. Вот как оно выглядит:

«Нужно создать устройство, на вход которого будет подаваться единичный импульс (функция единичного скачка) в результате чего на выходе будет появляться (генерироваться) синусоида».

Сознаюсь, функция единичного скачка (или же функция Хэвисайда) мною была выбрана не случайно. Нужно было подобрать наиболее простой для преобразования Лапласа  входной сигнал. Ну, а что может быть проще? Ничего! Ведь этот тип сигнала расположен первым в любой таблице преобразования Лапласа. Для нас этот сигнал можно объяснить так: «На вход устройства не подавалось напряжение и вдруг, мы его включили!». Как видите, пока всё довольно просто. Функция синусоиды, в свою очередь, тоже удобное значение и найти её преобразование в таблице преобразования Лапласа не составит труда. Перейдём к составлению выражения, связывающего выходной сигнал с входным. Это выражение удобно представить в виде передаточной функции:

H(s)=\frac{\L(sin(t))}{\L(H(t)) }

H(s)=\L(sin(t)) и \L(H(t)) обозначены преобразования Лапласа для синусоиды и функции Хэвисайда. Подставим значения взятые из таблицы преобразования Лапласа:

H(s)=\frac{\L(sin(t))}{\L(H(t)) }=\frac{\frac{1}{s^2+1}}{\frac{1}{s}}=\frac{s}{s^2+1}

Собственно говоря, нами с вами было получено выражение, описывающее передаточную функцию генератора синусоидального сигнала в частотной области комплексной переменной s. Нам же, людям, не живущим в измерении мнимых чисел, хочется получить выражение, которое можно было бы использовать на практике. Как бы ни прискорбно это звучало, но для достижения этой цели нам придется использовать билинейное преобразование (которое мы уже использовали в предыдущих статьях по проектированию фильтров). Напомню формулу билинейного преобразования:

s={\frac{2}{T}}{\frac{z-1}{z+1}}

Подставим данное выражение в нашу формулу, тем самым преобразуем его в дискретное Z-преобразование:

H(z)=\frac{({\frac{2}{T}}{\frac{z-1}{z+1}})}{({\frac{2}{T}}{\frac{z-1}{z+1}})^2+1}=\frac{{\frac{2}{T}}{\frac{z-1}{z+1}}{(T(z+1))^2}}{4(z^2-2Z+1)+T^2(z^2+2z+1)}=\frac{2(z-1)T(z+1)}{{z^2}(4+T^2)+z(2T^2-8)+(4+T^2)}

Разделим числитель и знаменатель на z^2:

H(z)=\frac{2T(1-z^{-2})}{(4+T^2)+z^{-1}(2T^2-8)+z^{-2}(4+T^2)}=\frac{Y(z)}{X(z)}

Y(z)[(4+T^2)+z^{-1}(2T^2-8)+z^{-2}(4+T^2)]=X(z)[2T(1-z^{-2})]

Взяв обратное преобразование Лапласа, перейдем во временную область, т. е. получим конечно-разностное выражение для генератора синусоиды:

y_t(4+T^2)+y_{t-1}(2T^2-8)+y_{t-2}(4+T^2)=2T(x_t-x_{t-2})

y_t=\frac{2T(x_t-x_{t-2})-y_{t-1}(2T^2-8)-y_{t-2}(4+T^2)}{(4+T^2)}=\frac{2(4-T^2)}{4+T^2}y_{t-1}-y_{t-2}

y_t=\frac{2(4-T^2)}{4+T^2}y_{t-1}-y_{t-2}

Получено желаемое уравнение! Дабы доказать читателю его жизнеспособность приведу несколько примеров его имплементации. Для этих целей воспользуемся языком программирования Python. В начале статьи я упоминал о том, что зная 2 предыдущие значения синусоиды можно рассчитать следующее. В реализации это будет выглядеть так:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 10 19:05:28 2015
 
 
@author:  Viktor Signaievskyi 
@blog:    www.HowToDoIt.com.ua 
@e-mail:  piligrim2007@meta.ua 
@details: Calculating sin(51deg) from sin(49deg) & sin(50deg)
 
Copyrighted by ©Viktor Signaievskyi 
"""
 
import math
import numpy as np
 
# We are supposing step is equal to 1degree
T = math.radians(1.0)
PrevBeforePrevValue = np.sin(math.radians(49.0))            # y(t-2)
PrevValue = np.sin(math.radians(50.0))                      # y(t-1)
ValueNowTrigonometricalFormula = np.sin(math.radians(51.0))
ValueNowRecursiveFormula = ((2.0*(4.0-T*T))/(4.0+T*T))*PrevValue - PrevBeforePrevValue
 
print("From RECURSIVE formula - " + str(ValueNowRecursiveFormula))
print("From TRIGONOMETRICAL formula - " + str(ValueNowTrigonometricalFormula))

 

Вывод отработавшего скрипта:

1
2
From RECURSIVE formula - 0.777145973303
From TRIGONOMETRICAL formula - 0.777145961457

Как видим, получена точность до 7го знака после запятой.

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import math
import numpy as np
import matplotlib.pyplot as plot
 
# It is important to remember that one period is equal to 2*pi = 6.28 radians!!!
 
N = 360 # Number of samples to generate (full period)
T = math.radians(360.0/N) # Sampling period
 
x = range(0, N-1)
y = range(0, N-1)
d = range(0, N-1)
 
y[0] = (2.0*T - 0     *(2*T*T-8) - 0*(4+T*T)) / (4.0+T*T)
y[1] = (2.0*T - y[1-1]*(2*T*T-8) - 0*(4+T*T)) / (4.0+T*T)
 
x[0] = np.sin(2.0*np.pi*0.0/N)
x[1] = np.sin(2.0*np.pi*1.0/N)
 
d[0] = y[0] - x[0]
d[1] = y[1] - x[1]
 
# Generating test signal
for t in range(2, N-1):
    x[t] = np.sin(2.0*np.pi*t/N)
    y[t] = ((2*(4-T*T))/(4.0+T*T))*y[t-1] - y[t-2]
    d[t] = y[t] - x[t]
 
# Drawing original and filtered signals
plot.clf()
plot.xlabel("Time")
plot.ylabel("Sample value")
plot.plot(x, 'r', lw=9)
plot.plot(y, 'ow', lw=3)
plot.plot(d, 'y')
 
plot.draw()

В результате, получаем следующий график:

Тут, белыми точками обозначена синусоида, полученная с помощью рекурсивной формулы. Кривая, обозначенная красным цветом, получена с помощью библиотечной функции синуса. Жёлтая кривая показывает разницу между двумя кривыми. Как видим, разница присутствует, но несущественная.

Возможно, единственными непонятными читателю строками является инициализация:

1
2
y[0] = (2.0*T - 0     *(2*T*T-8) - 0*(4+T*T)) / (4.0+T*T)
y[1] = (2.0*T - y[1-1]*(2*T*T-8) - 0*(4+T*T)) / (4.0+T*T)

Но и тут всё довольно просто! Поскольку была получена рекурсивная формула, для работы с ней нам необходимо иметь два предыдущих значения. Итого, я просто их рассчитал вручную.

Спасибо за внимание! Удачи в познаниях! Как всегда жду ваших отзывов внизу страницы:)