Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

welch_pgram is easily misused without a default window #152

Open
c42f opened this issue Feb 16, 2017 · 4 comments · May be fixed by #153
Open

welch_pgram is easily misused without a default window #152

c42f opened this issue Feb 16, 2017 · 4 comments · May be fixed by #153

Comments

@c42f
Copy link

c42f commented Feb 16, 2017

When n != nfft, getindex(::ArraySplit, i) fills the fft buffer with n elements rather than nfft. This can result in a lot of trailing zeros in the buffer, so periodograms on irregular length data are likely to be contaminated with artifacts from the frequency response of a boxcar window.

On a related note, it would be great if we could change the default windowing function to something better. I'm aware that there can be subtle tradeoffs when choosing a window, but just about any choice will be a better default than what we have now: currently it's pretty much mandatory to set the window by hand, or the results will contain a lot of spurious ringing. I don't know a lot about window choice, but matlab seems to default to the Hamming window for pwelch().

@simonster
Copy link
Member

What would you expect the behavior to be if n != nfft?

@c42f
Copy link
Author

c42f commented Feb 16, 2017

A very good question, I now suspect this issue is more of a misunderstanding on my part about what should be expected from the function than any real problem with the code. A standard case of pebkec I suppose :-/

To expand on what I was trying to do - I've got samples of a non-stationary time series of varying lengths (actually the pitch of an aircraft). I've estimated the power spectrum using welch_pgram() to understand some vibration issues we're having. I don't do a lot of signal processing, so I got the windowing wrong to start with. In fact, I ended up having to implement my own my_welch() function to understand how to use the optimized and nicely abstracted version in this package.

So, what went wrong?

  • It's quite likely that n != nfft, due to the way that defaults are chosen for the function, which means there will be some of zero padding on average; the amount of zero padding can change a lot as n is changed.
  • In my case the signal has a long term trend, so chunks of the signal have a large DC offset, and sensible windowing is particularly important to avoid the spectrum being rubbish in the high frequencies. In particular, this seems to mean choosing a window which drops to zero (ruling out the Hamming window as a reasonable choice, and certainly ruling out the default of no window).

I'll double check this tomorrow.

@simonster
Copy link
Member

It's likely that we should apply some window by default if none is specified, since, as you say, this is what other implementations do.

@c42f c42f changed the title welch_pgram broken for n != nfft welch_pgram is easily misused without a default window Feb 17, 2017
@c42f
Copy link
Author

c42f commented Feb 17, 2017

The corresponding code in octave and scipy uses the Hamming and Hann(ing) windows respectively. Octave has some good commentary on the effect of window choice, zero padding, etc, and also mentions removal of the mean by default, which seems like something that might matter a lot in some cases.

https://sourceforge.net/p/octave/signal/ci/default/tree/inst/pwelch.m#l197

It looks like choosing good default settings is a tricky matter.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants