# Compute Hann window without cos function

In an environment with limited memory and computing power it is interesting to be able to generate a Hann window without using a cache or repetitive calling of expensive functions such as sine and cosine functions.

The following document demonstrate the use of an alternative approach in order to compute a pure tone by storing only a coefficient and the two last results:

Rewrite using python code:

import math
import matplotlib.pyplot as plt

tone = [0.0] * 50
tone_freq = 1720.0
sample_rate = 44100.0
ffs = tone_freq / sample_rate;
k1 = math.cos(2 * math.pi * ffs)
k2 = math.sin(2 * math.pi * ffs)
tone[1] = k2
for i in range(2, len(tone)):
tone[i] = 2 * k1 * tone[i-1] - tone[i-2]

plt.plot(tone)
plt.show()


Using the same approach is it possible to generate the Hann window ?

Note:

You can generate a Hann window with sin function like that:

tone = [(0.5 - 0.5 * math.sin(0.5 * math.pi +( 2 * math.pi * i) / (len(tone) - 1))) for i in range(len(tone))]


Are you seeking a solution (you seem to have found one) or understanding as well?

Do you know complex numbers?

Your technique is based on the properties of the complex unit circle. Suppose you have a point on the complex unit circle:

$$u = a + bi,\; a^2+b^2=1$$

If you square that point:

$$u^2 = ( a^2 - b^2 ) + 2abi$$

The result will also lay on the unit circle twice as far along the perimeter. If you cube, three times, etc.

$$z_n = u^n$$

If you generate a sequence, you will generate a set of evenly spaced points around the circle. If your interval is a whole number fraction $$(1/N)$$ of the circumference, then the $$N$$th point will land back on the origin and $$u$$ is then called the $$N$$th Root of Unity.

This is also the principle behind "When you multiply two complex numbers, your multiply the magnitudes adn add the angles." So each time you multiply by $$u$$ you are adding its angle.

RB-J's answer is a slick recursive improvement on the principle which is the basis for Goertzel's bin calculation formula.

You will need to do one or two trig calculations before the loop. Since the angle will be small for those, Taylor series approximations will only take a few terms.

See these articles of mine for further explanation.

P.S. If you are using the DSP version of the VonHann for a DFT, you should be using N in the denominator, not (N-1).

The proof/derivation of RB-J's and your solution:

$$y[n] = \cos( \omega n )$$

$$y[n-1] = \cos( \omega (n-1) ) = \cos( \omega n )\cos( \omega ) + \sin( \omega n )\sin( \omega )$$

\begin{aligned} y[n-2] &= \cos( \omega n )\cos( 2\omega ) + \sin( \omega n )\sin( 2\omega ) \\ &= \cos( \omega n ) \left( 2\cos^2( \omega ) - 1 \right) + \sin( \omega n )2\sin( \omega )\cos( \omega ) \\ &= \cos( \omega n ) \left( 2\cos^2( \omega ) - 1 \right) + \left( y[n-1] - \cos( \omega n )\cos( \omega ) \right) 2\cos( \omega ) \\ &= y[n] \left( 2\cos^2( \omega ) - 1 \right) + \left( y[n-1] - y[n] \cos( \omega ) \right) 2\cos( \omega ) \\ &= -y[n] + y[n-1] 2\cos( \omega ) \\ \end{aligned}

So the recursion becomes:

$$y[n] = 2\cos( \omega ) y[n-1] - y[n-2]$$

Starting with:

\begin{aligned} y[0] &= 1\\ y[1] &= \cos( \omega ) \\ \end{aligned}

Here is an updated program and results. Notice that the unit circle approach is more accurate. Worth it? You decide.

Program:


import numpy as np

#==========================================================
def main():

#---- Set parameter

N = 20

#---- Initialize

slice = 2.0 * np.pi / N

u_real = np.cos( slice )
u_imag = np.sin( slice )

z_real = 1.0
z_imag = 0.0

#---- Goertzel generation

k = u_real
y = np.zeros( N )
y[0] = 1.0
y[1] = k

k2 = 2.0 * k

for n in range( 2, N ):
y[n] = k2 * y[n-1] - y[n-2]

#---- Unit circle powers generation

for n in range( N ):
angle = n * slice
c = np.cos( angle )
s = np.sin( angle )

print "%3d  %6.3f %6.3f  %6.3f %6.3f %6.3f  %11.4e %11.4e"  \
% ( n, c, s, z_real, z_imag, y[n], c - z_real, c - y[n] )

next_real = u_real * z_real - u_imag * z_imag    # cos(A) cos(B) - sin(A) sin(B)
next_imag = u_real * z_imag + u_imag * z_real    # cos(A) sin(B) + sin(A) cos(B)

z_real = next_real
z_imag = next_imag

#==========================================================
main()

#
# Complex multiplication
#
#    A = cA + i sA
#    B = cB + i sB
#
#   AB = ( cA + i sA )( cB + i sB )
#      = cA cB + i cA sB + i sA cB + i^2 sA sB
#      = (cA cB - sA sB ) + i ( cA sB + sA cB )
#      = cAB + i sAB
#
#   cAB = cA cB - sA sB
#   sAB = cA sB + sA cB
#


Results:

  0   1.000  0.000   1.000  0.000  1.000   0.0000e+00  0.0000e+00
1   0.951  0.309   0.951  0.309  0.951   0.0000e+00  0.0000e+00
2   0.809  0.588   0.809  0.588  0.809   1.1102e-16  2.2204e-16
3   0.588  0.809   0.588  0.809  0.588   0.0000e+00  3.3307e-16
4   0.309  0.951   0.309  0.951  0.309  -1.1102e-16  4.4409e-16
5   0.000  1.000   0.000  1.000 -0.000  -1.0530e-16  5.0532e-16
6  -0.309  0.951  -0.309  0.951 -0.309  -1.6653e-16  4.9960e-16
7  -0.588  0.809  -0.588  0.809 -0.588  -3.3307e-16  4.4409e-16
8  -0.809  0.588  -0.809  0.588 -0.809  -4.4409e-16  1.1102e-16
9  -0.951  0.309  -0.951  0.309 -0.951  -5.5511e-16 -2.2204e-16
10  -1.000  0.000  -1.000  0.000 -1.000  -5.5511e-16 -6.6613e-16
11  -0.951 -0.309  -0.951 -0.309 -0.951  -5.5511e-16 -1.2212e-15
12  -0.809 -0.588  -0.809 -0.588 -0.809  -4.4409e-16 -1.6653e-15
13  -0.588 -0.809  -0.588 -0.809 -0.588  -1.1102e-16 -1.9984e-15
14  -0.309 -0.951  -0.309 -0.951 -0.309   1.6653e-16 -2.3315e-15
15  -0.000 -1.000  -0.000 -1.000  0.000   3.7141e-16 -2.5152e-15
16   0.309 -0.951   0.309 -0.951  0.309   6.1062e-16 -2.4425e-15
17   0.588 -0.809   0.588 -0.809  0.588   6.6613e-16 -2.1094e-15
18   0.809 -0.588   0.809 -0.588  0.809   8.8818e-16 -1.4433e-15
19   0.951 -0.309   0.951 -0.309  0.951   9.9920e-16 -6.6613e-16


Appendix I. A memory preserving version of the Geortzel

        N = 20

slice = 2.0 * np.pi / N

k = np.cos( slice )
k2 = 2.0 * k

yn2 =  1.0
yn1 =  k

print 0, yn2
print 1, yn1

for n in range( 2, N ):
yn = k2 * yn1 - yn2
yn2 = yn1
yn1 = yn
print n, yn



Appendix II. Cosine approximation that could be done with paper and pencil

Using Taylor series to generate the initial cosine value for various values of N.

If you need more accuracy at low N, just increase the number of terms following the pattern.

import numpy as np  # needed for cosine comparison

#==========================================================
def main():

pi = 3.14159265358979323846264

f2 =       0.5  #   1 / 2
f4 = f2 / 12.0  # / 3 / 4
f6 = f4 / 30.0  # / 5 / 6
f8 = f6 / 56.0  # / 7 / 8

for N in range( 10, 100, 10 ):

x = 2.0 * pi / N   # Angle in Radians

x2 = x  * x
x4 = x2 * x2
x6 = x2 * x4
x8 = x4 * x4

c = 1 - x2*f2 + x4*f4 - x6*f6 + x8*f8
#         s = x - x3*f3 + x5*f5 - x7*f7 + x9*f9

k = np.cos( x )

print "%3d  %14.10f  %14.10f %8.1e" % ( N, c, k, c - k )

#==========================================================
main()


Look, it's actually a dot product!

Results:

 10    0.8090169970    0.8090169944  2.6e-09
20    0.9510565163    0.9510565163  2.6e-12
30    0.9781476007    0.9781476007  4.5e-14
40    0.9876883406    0.9876883406  2.6e-15
50    0.9921147013    0.9921147013  2.2e-16
60    0.9945218954    0.9945218954  1.1e-16
70    0.9959742940    0.9959742940  0.0e+00
80    0.9969173337    0.9969173337  0.0e+00
90    0.9975640503    0.9975640503  0.0e+00

• What column is the Hann window ? – nicolas-f Jun 24 '20 at 8:40
• Ok I understand (missing the 0.5 - 0.5 * y(n) ). Yes your solution is much more accurate and require less memory. Thank you. – nicolas-f Jun 24 '20 at 9:35
• @nicolas-f You're welcome. I appended an arrayless version. – Cedron Dawg Jun 24 '20 at 14:31
• @nicolas-f To complete the library-less ensemble, I appended a cosine calculation as well. Have fun! – Cedron Dawg Jun 24 '20 at 15:26

i think this generates a sinusoid:

$$y[n] = 2 \cos(\omega_0) y[n-1] - y[n-2]$$

if you initialize

$$y[-2] = \cos(2 \omega_0)$$ $$y[-1] = \cos(\omega_0)$$

then the result is

$$y[n] = \cos(\omega_0 n) \qquad n \ge 0$$

the period of this cosine function is $$N=\frac{2\pi}{\omega_0}$$ and you want $$N-1$$ to be the non-zero width of the window.

finally you construct your window with

$$w[n] = \tfrac{1}{2} (1 - y[n]) \qquad 0 \le n \le N$$

• I'm sorry, I am maths novice, I don't know what to do with your answer. k1 and k2 must use both cos function ? I can't generate any usable window with your formulas. – nicolas-f Jun 22 '20 at 16:22
• In Matlab-esque: w0 = 2*pi/100; a = [2*cos(w0) -1]; b = 1; y(1:2)=cos(w0*[2 1]); T=1e2; for t=3:T+2, y(t)=2*cos(w0)*y(t-1)-y(t-2); end, y=y(3:end); max(abs(y-cos(w0*(0:T-1))')) I am curious what happens when implementing this in floating point if one lets the sinoid oscillate for a long time? – Knut Inge Jun 23 '20 at 22:02
• Yes the error is accumulating, this is why double precision is required. – nicolas-f Jun 24 '20 at 7:30
• how big is $N$? and how many different values of $N$ are you going to need? $\cos(\frac{2\pi}{N})$ and $\cos(\frac{4\pi}{N})$ and $2\cos(\frac{2\pi}{N})$ can all be pre-calculated and put into a table if there are only a few different values of $N$ you are dealing with. – robert bristow-johnson Jun 24 '20 at 20:29

With the help of a colleague we come up to this answer:

import math
import matplotlib.pyplot as plt
import numpy

window_size = 31
hann = numpy.zeros(window_size)
expected = numpy.hanning(window_size)

window_wT = math.pi / (window_size - 1)
window_k1 = math.cos(window_wT)
window_k2 = math.sin(window_wT)
window_k3 = 0
hann[1] = window_k2 * window_k2
for i in range(2, window_size):
w_nk2 = window_k2
window_k2 = 2 * window_k1 * window_k2 - window_k3
window_k3 = w_nk2
hann[i] = window_k2 * window_k2

fig, axs = plt.subplots(2, 1)
axs[0].plot(expected)
axs[0].plot(hann)
axs[0].set_ylabel('Value')
axs[1].plot(numpy.abs(expected - hann))
axs[1].set_ylabel('Error')
plt.show()