[Octave-bug-tracker] [bug #43305] Hamming etc. windows are wrong

[Oscar Ruitt]

URL:
  <http://savannah.gnu.org/bugs/?43305>

                 Summary: Hamming etc. windows are wrong
                 Project: GNU Octave
            Submitted by: oscarruitt
            Submitted on: Sat 27 Sep 2014 12:23:41 AM GMT
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Jerry
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.8.0
        Operating System: Any

    _______________________________________________________

Details:

I’ve noticed that the Hamming window (and I suppose other windows) provided
by Octave, SciPy****, and NumPy is wrong, returning a window that is symmetric
rather than one that is “DFT even.” This is easily spotted by looking at
the first and last points of the window; if they are the same values, then the
window is incorrect.

These windows are normally applied to data that are equally spaced, thus
implying that their DFTs are periodic. Let’s assume that we are windowing N
data points, indexed from 0 to N-1 as is normal in signal processing. Then the
DFT is also N points long, also indexed from 0 to N-1. The periodic extension
principal implies that the point indexed by N in the DFT domain has the same
value as the point indexed by 0. Thus any window applied to the data should
have length N and have the same period.

This mistake is widespread, infecting many software packages. (I am sending
this same note to the NumPy, SciPy, and Octave lists.) It is certainly wrong
in most textbooks, including the classic, Oppenheim & Schafer, 1975. However,
the definitive and widely referenced paper on windows by Fredric J. Harris,
“On the Use of Windows for Harmonic Analysis with the Discrete Fourier
Transform,” Proceedings of the IEEE, vol. 66, NO. 1 , January, 1978 *
discusses this problem, both its roots and its correction. For example,
quoting from the paper:

“Since the DFT essentially considers sequences to be periodic, we can
consider the missing end point to be the beginning of the next period of the
periodic extension of this sequence. In fact, under the periodic extension,
the next sample (at 16 s in Fig. 1.) is indistinguishable from the sample at
zero seconds.

“This apparent lack of symmetry due to the missing (but implied) end point
is a source of confusion in sampled window design. This can be traced to the
early work related to convergence factors for the partial sums of the Fourier
series. The partial sums (or the finite Fourier transform) always include an
odd number of points and exhibit even symmetry about the origin. Hence much of
the literature and many software libraries incorporate windows designed with
true even symmetry rather than the implied symmetry with the missing end
point!”

...and later...

“We will now catalog some well-known (and some not well-known windows. For
each window we will comment on the justification for its use and identify its
significant parameters. All the windows will be presented as even (about the
origin) sequences with an odd number of points. To convert the window to DFT
even, the right endpoint will be discarded and the sequence will be shifted so
that the left end point coincides with the origin.”

Lest one consider this a trivial difference, consider the following
Octave-Matlab code (corresponding NumPy and SciPy code would yield the same
results).


x = ones(8, 1); # An even sequence; its DFT should be real.
hWrong = hamming(8)
h9 = hamming(9);
hRight = h9(1:8)
xWrong = x .* hWrong; # Multiplication by ones shown
xRight = x .* hRight; # for clarity.
XWrong = fft(xWrong) # Contains nonzero imaginary part.
XRight = fft(xRight) # Real, as expected from evenness of x.
XWrongMag = abs(XWrong) # The magnitudes are
XRightMag = abs(XRight) # also different.


This causes the following output:


hWrong =

  0.080000
  0.253195
  0.642360
  0.954446
  0.954446
  0.642360
  0.253195
  0.080000

hRight =

  0.080000
  0.214731
  0.540000
  0.865269
  1.000000
  0.865269
  0.540000
  0.214731

XWrong =

  3.86000 + 0.00000i
 -1.76795 - 0.73231i
  0.13889 + 0.13889i
  0.01906 + 0.04602i
  0.00000 + 0.00000i
  0.01906 - 0.04602i
  0.13889 - 0.13889i
 -1.76795 + 0.73231i

XRight =

  4.32000 + 0.00000i
 -1.84000 + 0.00000i
  0.00000 + 0.00000i
  0.00000 + 0.00000i
  0.00000 + 0.00000i
  0.00000 - 0.00000i
  0.00000 - 0.00000i
 -1.84000 - 0.00000i

XWrongMag =

  3.86000
  1.91362
  0.19642
  0.04981
  0.00000
  0.04981
  0.19642
  1.91362

XRightMag =

  4.32000
  1.84000
  0.00000
  0.00000
  0.00000
  0.00000
  0.00000
  1.84000


This news will likely upset many and I will be interested to see what
arguments will flow in favor of keeping the status quo--I’m sure one will
be, “we’ve always done it this way” and others will be more substantial,
possibly expressing decent arguments why the current situation is correct for
some uses. In any case, I refer again to the Harris paper.

So rather than host an argument on this list, this is what I propose: Do what
Matlab** and SciPy*** **** do. Acknowledge both uses by making a flag to
handle the “periodic” or the “symmetric” cases. The Matlab default is
“symmetric” which is of course unfortunate but at least such inclusion in
Octave and SciPy would retain compatibility with the existing usage. Then
it’s up to the user whether to shoot him/herself in the foot, assuming that
such a decision is guided by actually referring to the documentation for the
package being used, since the default behavior is wrong.

Jerry Bauck


*
http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1455106&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D1455106

** http://www.mathworks.com/help/signal/ref/hamming.html

***
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.hamming.html

**** I am just learning Python-NumPy-SciPy but it appears as though the SciPy
situation is that the documentation page above *** mentions the flag but the
flag has not been implemented into SciPy itself. I would be glad to stand
corrected.




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?43305>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/