Again and again, this post for me is none other than the cry of my soul. I am struck by the fact that some things from everyday life used by us for granted. We believe even in lie. Take the formula, substitute the value and it works (so far, seems to be :). People, stop! “But all is working right!” – You might say. “Smart guy from the forum said that everything will work perfectly!”. Perhaps smart guy is right, but as the classic said “Do not trust anyone!”. Should you trust me? – It is really possible, but after rechecking. After all the task of the engineer not only to create something but also deeply understand the running process. As an example of “it should work!”, we will consider the recursive formula for the sine wave generation, which I have repeatedly seen on many websites and blogs dedicated to digital signal processing. Naturally, I checked this formula in practice and was really surprised by such simplicity. It really works! I was so amazed that could not resist and set out myself to understand the magic of the miracle formula. I was not very surprised when the clear and simple way of the formula deduction was not found in the network. Found a lot of absurdities, empirical and deductive methods and other “rubbish” but for myself – NOTHING.

For experts (somehow reading this) I have to say that the method of obtaining a formula that is based on 2 poles immediately from Z plane and taken symmetrically with respect to the real coordinate axis specifically for me was absolutely incomprehensible method, with many author’s assumptions taken by me again on faith. Newcomers to the DSP immediately ask about why take two poles and not, say four or just one. The next question may be a question about the number and location of zeros on Z plane… Oh well. It was just for experts!

But seriously, the formula was really deduced in simple way by me. The process of deduction I share with you. What do you get as a result? Simple recursive formula, which can be used to generate a sine wave without need to call trigonometric functions from the library, which by the way, will run slower at times comparing to the proposed method (which is, for example, for the controller is very important). This is really interesting fact that knowing two previous values of the sinusoid the next value can be easily calculated.

So, let us start!

Designing for myself, I started with the so-called technical specifications. Here’s how it looks:

“We need to create a device to the input of which unit step function will be supplied. As a result sine wave must appear on the output.”

It is some kind of Black box. Looks simple, isn’t it? I confess, the unit step (or Heaviside) function, I was not chosen at random. It was necessary to find the simplest input signal for the Laplace transform. Well, what could be easier? Nothing, at all! After all, this type of signal is first in any table of the Laplace transform. For us, this signal can be explained as follows: “The input device suddenly turns voltage on!”. As you can see, it’s pretty simple. Sine function, in turn, is also easy to find in Laplace transform conversion table. Let’s move on to the preparation of expression which shows us the relation between output and the input of black box. This expression can be conveniently represented in the form of the transfer function:

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

H(s)=\L(sin(t)) and \L(H(t)) are Laplace transform for the sine wave and Heaviside functions. Let us substitute expressions taken from the Laplace transform table:

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

In fact, we have obtained an expression that describes the transfer function of the sine wave signal generator in the frequency domain for the complex variable s. We are not the people who are living in the dimension of imaginary numbers! We have to get an expression that could be used in practice. No matter how sad it sounds, but to achieve this we need to use bilinear transformation (which we have already used in previous articles on designing filters). Let me remind you the formula of bilinear transform:

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

Substituting this expression into previous formula we are converting it into discrete Z-transform:

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)}

Divide the numerator and denominator by 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})]

Taking the inverse Laplace transform, we go into the time domain, i.e. we obtain the finite difference expression for the sine wave generator:

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}

Desired equation was achieved! In order to prove its validity to the reader here are a few examples of its implementation. For this purpose, we use Python programming language. In the beginning of this article I mentioned that knowing the two previous values of sine wave the next value can be calculated easily. The implementation of this would look like:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# -*- coding: utf-8 -*-
"""
Created on Sun Feb 15 17:02:08 2015
 
Calculating sin(51deg) from sin(49deg) & sin(50deg)
 
@author: piligrim
"""
 
import math
import numpy as np
 
# We are supposing step is equal to 1degree
T = math.radians(1.0/360.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))

The Python output is:

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

As you can see, we’ve got an accuracy of 3 decimal places, which is good.

For those who believe the correct result coincidence, I cite a more interesting example. Using the above formulas we derive the full period of a sine wave and display the data to the plot. On the same graph display the sine wave with the same parameters, but obtained using the sine function from the library. It is not difficult to deduce the graph error. Below you can find the source code:

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()

As a result, we obtain the following plot:

Here, white dots denote a sine wave obtained using a recursive formula. The curve marked in red, obtained by sine function from the library. The yellow curve shows the difference between the two curves. As you can see, the difference is present, but not significant.

Perhaps the only confusing to the reader is the initialization string:

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)

But here, it’s pretty easy! As recursive formula has been successfully deduced, to work with it we need to have the two previous values. As a result I just took them manually (but using the same formula).

Thank you for your attention! Good luck! As always, I welcome your feedback at the bottom of the page :)