[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Hamming etc. windows are wrong
From: |
Jerry |
Subject: |
Hamming etc. windows are wrong |
Date: |
Fri, 26 Sep 2014 01:41:55 -0700 |
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 symmetric.” 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. The window can also be imagined to have the
same kind of periodic extension, so that what would be the N-th point of the
window matches the N-point of the data with the same weight as was applied to
the 0-th point of the data but without duplicating that point at N-1.
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.”
Note that he limits himself therefore to even-numbered "DFT even" windows, but
the rule remains--honor the implied periodicity even for windows of odd length.
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"
xRight = x .* hRight; # shown 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 useful for some
applications. 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** does. 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, NumPy, 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 and not
blindly using the default.
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.
***** http://docs.scipy.org/doc/numpy/reference/generated/numpy.hamming.html
- Hamming etc. windows are wrong,
Jerry <=