In this tutorial http://www.ee.ic.ac.uk/pcheung/teaching/ee3_Study_Project/Sinewave%20Generation(708).pdf
an algorithm to generate a sinewave through an IIR filter is described. On page 2 it is stated: "The formal proof can be found in the DSP related text book". Could anyone please suggest me a good book where the formal proof can be found?
There are different ways to look at this. Here are two.
The IIR interpretation is that you put a pole on the unit circle at the desired (normalized) frequency $\omega = 2\pi \cdot f_0/f_s$, where $f_0$ is the oscillator frequency and $f_s$ is the sample rate.
Since the output needs to be real you need a pair of conjugate poles and the transfer function becomes
$$H(z) = \frac{1}{(1-e^{j\omega}z^{-1})\cdot (1-e^{-j\omega}z^{-1})} = \frac{1}{1-2cos(\omega)z^{-1}+z^{-2}}$$ That's an IIR filter with $a_1 = -2cos(\omega)$ and $a_2=1$ so the time difference equation becomes
$$y[n] = 2cos(\omega)y[n-1] - y[n-2]$$
Note that there is no input and that the filter has infinite gain at the oscillating frequency, so this only works if it's properly initialized. If you want a cosine you would do
$$y[0] = 1, y[1] = cos(\omega), y[n] = 2cos(\omega)y[n-1] - y[n-2], n > 1$$
You can just also proof the recursion using the trig identities. We need to show that
$$cos(\omega \cdot (n+1)) = 2 \cdot cos(\omega) \cdot cos(\omega\cdot n )-cos(\omega \cdot (n-1))$$
We can substitute $x = \omega \cdot n$ and $y = \omega$ and start from the right side of the equation:
$$2cos(x)cos(y)-cos(x-y) = \\ 2cos(x)cos(y)-[cos(x)cos(y)-sin(x)sin(y)] = \\cos(x)cos(y)+sin(x)sin(y) = cos(x+y) \\q.e.d. $$
Note that this algorithm is susceptible to numerical noise. Since it's purely recursive numerical issues like rounding or truncating can accumulate over time and result in amplitude drift.
An alternative algorithm would be based on phasor rotation. You use a complex state variable $X = a + jb$. It is initialized to $X[0] = 1$ and the recursion becomes
$$X[n] = X[n-1] \cdot e^{-j\omega}, n> 0$$
The real part is a cosine oscillator and the imaginary part a sine, so you get two sine waves that are 90 degrees apart.
Since it's also recursive it's also vulnerable to drift, but it can easily be "corrected". We know that the magnitude squared of the state variable should always be unity. We can apply an amplitude correction as $$X'[n] = X[n] \cdot \frac{1}{\sqrt{a^2+b^2}}$$
The magnitude will always be very close to unity, so we can approximate this through a Taylor expansion to
$$X'[n] = X[n] \cdot \frac{1}{\sqrt{a^2+b^2}} \approx X[n] \cdot 0.5 \cdot (3- (a^2+b^2) )$$
which is much cheaper to implement.
Since any drift would be fairly slow, this correction only needs to be applied every few hundred samples or so and it guarantees the oscillator to be stable forever.
