0
$\begingroup$

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:

http://ww1.microchip.com/downloads/en/appnotes/00543c.pdf

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))]
| improve this question | |
$\endgroup$
2
$\begingroup$

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

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$$

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – nicolas-f Jun 22 at 16:22
  • $\begingroup$ 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? $\endgroup$ – Knut Inge Jun 23 at 22:02
  • $\begingroup$ Yes the error is accumulating, this is why double precision is required. $\endgroup$ – nicolas-f Jun 24 at 7:30
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson Jun 24 at 20:29
1
$\begingroup$

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

enter image description here

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()
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.