bugGNU Octave - Bugs: bug #43305, Hamming etc. windows are wrong

 
 

bug #43305: Hamming etc. windows are wrong

Submitter:  Oscar Ruitt <oscarruitt>
Submitted:  Sat 27 Sep 2014 12:23:41 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Jerry Open/Closed:  * Closed
Release:  * 3.8.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 22 Feb 2015 10:56:21 PM UTC, comment #10: 

I have fixed this for the following window functions in Octave: blackman, hamming, and hanning. These functions now take the optional argument "periodic" or "symmetric", which is the Matlab way to specify which form of the window should be returned.

http://hg.savannah.gnu.org/hgweb/octave/rev/26fb4bfa4193

Closing this bug as fixed (in Octave). I will open another bug report about adding the periodic types to the window functions provided by the signal package.

And if someone wants to add other options in place of "periodic" and "symmetric", this should be made yet another bug report.

Mike Miller <mtmiller>
Group Member
Tue 30 Dec 2014 04:55:57 PM UTC, comment #9: 

I seem to remember this coming up on the list before, and I made a note somewhere to revisit the window functions and add the "periodic" option to the relevant functions. Obviously hasn't happened yet :)

I think that the default should be compatible with Matlab, which seems to be "symmetric", with the "periodic" option added for at least those that have it in Matlab.

Mike Miller <mtmiller>
Group Member
Tue 30 Dec 2014 04:15:32 PM UTC, comment #8: 

The default setting is a very important decision and worthy of a discussion on the Octave Maintainer's list.

@Doug: Could you post something on the problem, the resolution, and the pros/cons?

Rik <rik5>
Group administrator
Tue 30 Dec 2014 11:27:31 AM UTC, comment #7: 

I vote to do the mathematically correct way ( opposite Matlab )
and have a very clear statement in the help section.

Doug Stewart <dastew>
Tue 30 Dec 2014 09:07:15 AM UTC, comment #6: 

I agree—same length and provide both types. I think that is what Matlab does, and defaults to symmetric. Like I said, this is unfortunate since one would like the default to be the correct window. I suppose this all comes down to how rigorously Octave follows Matlab and I’m not steeped enough in the Octave culture to be of any help there.

Oscar Ruitt <oscarruitt>
Tue 30 Dec 2014 12:02:06 AM UTC, comment #5: 

THe short answer is as you recommend: provide
both symmetric and periodic and return 8 for length
8 for both.

AND, an explanation in the Manual would help, with a
reference to Harris.

If anyone has an argument for another choice, it
should be considered.

Users may want both choices, if only to reproduce
someone else's results...

Michael Godfrey <godfrey>
Group Member
Mon 29 Dec 2014 08:20:18 PM UTC, comment #4: 

Oscar stated:

"
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.
"

So each window function now has 1 input  (length)

1) Should we now have 2 inputs  (length , type)
   where  type =  “periodic” or “symmetric”
   Which should be the default?

2) What length do we want?
   example:
     length =8
     Should we produce 8 for both periodic and symmetrical?
     or 7 for periodic and 8 for symmetrical
    8 for symmetrical is what we do now and would have to stay
    that way for compatibility with our old way and with Matlab.
     If we did 7 then the actual numbers would be the same
     between the 2 methods. With 8 the numbers are different.
  My preference is 8 for both. 


Mike what are your thoughts?

Doug Stewart <dastew>
Mon 29 Dec 2014 07:05:58 PM UTC, comment #3: 

I'm adding Mike Miller to the CC list because he also has a large interest in getting signal processing to be correct in Octave.

Rik <rik5>
Group administrator
Mon 29 Dec 2014 06:02:12 PM UTC, comment #2: 

I have an original of the Harris paper!!!

If this needs fixing in octave, then I will
try and do it.
I will look at it tomorrow.

Doug Stewart <dastew>
Sat 27 Sep 2014 10:57:26 AM UTC, comment #1: 

This is very useful, and the Fred Harris paper is a
good reference and
he quotes other useful contributions.

One minor point: the Harris paper is not under copyright,
but the IEEE trys to charge non-members for it.
A free source is:
https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf

Michael

Michael Godfrey <godfrey>
Group Member
Sat 27 Sep 2014 12:23:41 AM UTC, original submission:  

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://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.

Oscar Ruitt <oscarruitt>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by oscarruitt (Submitted the item)
  • -email is unavailable- added by oscarruitt
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-02-22 mtmiller StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2014-12-29 rik5 Carbon-Copy- Added mtmiller
    2014-12-29 rik5 StatusNone Confirmed
    2014-09-27 oscarruitt Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code