JuliaDSP / DSP.jl

Filter design, periodograms, window functions, and other digital signal processing functionality
https://docs.juliadsp.org/stable/contents/
Other
381 stars 110 forks source link

welch_pgram is easily misused without a default window #152

Open c42f opened 7 years ago

c42f commented 7 years ago

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 commented 7 years ago

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

c42f commented 7 years ago

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?

I'll double check this tomorrow.

simonster commented 7 years ago

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 commented 7 years ago

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.