Polyphase FFT

Started by Alberto November 1, 2005
I have a few programs that compute and display a running waterfall
of an audio signal input to the PC sound card. The waterfall is
computed via ordinary FFTs up to 1048576 points for the highest
resolution settings.
Recently I have seen a similar program that uses polyphase FFTs
and it has a spectral leakage much lower than mine (I use a Blackman
window). The subject of polyphase FFT is new to me, and from a cursory
investigation with Google, it looks like (but I am certainly wrong)
the signal is passed through a bank of FIR filters whose width is
comparable to a single FFT bin. This looks almost impossible to me,
as the time needed to compute 1048576 such FIR filters would make
the program not suitable for real time operation.

Is there anybody willing to clarify a little this subject, maybe 
with hints or examples on how to implement in C code a polyphase
FFT whose goal is to compute magnitudes to be used for the waterfall ?

Thanks indeed

Alberto
Alberto wrote:
> I have a few programs that compute and display a running waterfall > of an audio signal input to the PC sound card. The waterfall is > computed via ordinary FFTs up to 1048576 points for the highest > resolution settings. > Recently I have seen a similar program that uses polyphase FFTs > and it has a spectral leakage much lower than mine (I use a Blackman > window). The subject of polyphase FFT is new to me, and from a cursory > investigation with Google, it looks like (but I am certainly wrong) > the signal is passed through a bank of FIR filters whose width is > comparable to a single FFT bin. This looks almost impossible to me, > as the time needed to compute 1048576 such FIR filters would make > the program not suitable for real time operation. > > Is there anybody willing to clarify a little this subject, maybe > with hints or examples on how to implement in C code a polyphase > FFT whose goal is to compute magnitudes to be used for the waterfall ? > > Thanks indeed > > Alberto
I'm not familiar with the term polyphase FFT, but I'll suggest a technique that can reduce leakage beyond what is possible using an ordinary window. It is called Weighted Overlap Add (WOLA) and is described in an out-of-print book called Multirate DSP by Crochiere and Rabiner. For an FFT size of M, a 4:1 WOLA is computed by windowing off a chunk of 4M samples, then stacking and adding each length M segment of the windowed data together to make a length M time-aliased vector -- your FFT input. 50% overlap is used by sliding 2M new samples into the next window calculation. John
I may understand the reason about doing the N:1 weighted-overlapping.
Since the book is out-of-print, i do not have more information about
detail.
The question in my mind is how to make sure the weighting do not
corrupt the continuity? also not to introduce unwanted spectrum from
next frame?

john wrote:
> Alberto wrote: > > I have a few programs that compute and display a running waterfall > > of an audio signal input to the PC sound card. The waterfall is > > computed via ordinary FFTs up to 1048576 points for the highest > > resolution settings. > > Recently I have seen a similar program that uses polyphase FFTs > > and it has a spectral leakage much lower than mine (I use a Blackman > > window). The subject of polyphase FFT is new to me, and from a cursory > > investigation with Google, it looks like (but I am certainly wrong) > > the signal is passed through a bank of FIR filters whose width is > > comparable to a single FFT bin. This looks almost impossible to me, > > as the time needed to compute 1048576 such FIR filters would make > > the program not suitable for real time operation. > > > > Is there anybody willing to clarify a little this subject, maybe > > with hints or examples on how to implement in C code a polyphase > > FFT whose goal is to compute magnitudes to be used for the waterfall ? > > > > Thanks indeed > > > > Alberto > > I'm not familiar with the term polyphase FFT, but I'll suggest a > technique that can reduce leakage beyond what is possible using an > ordinary window. It is called Weighted Overlap Add (WOLA) and is > described in an out-of-print book called Multirate DSP by Crochiere and > Rabiner. For an FFT size of M, a 4:1 WOLA is computed by windowing off > a chunk of 4M samples, then stacking and adding each length M segment > of the windowed data together to make a length M time-aliased vector -- > your FFT input. 50% overlap is used by sliding 2M new samples into the > next window calculation. > > John
cuteworm@wildmail.com wrote:
> I may understand the reason about doing the N:1 weighted-overlapping. > Since the book is out-of-print, i do not have more information about > detail. > The question in my mind is how to make sure the weighting do not > corrupt the continuity? also not to introduce unwanted spectrum from > next frame?
The 4M samples are multiplied by a window function that tapers at the endpoints, the first 25%, second 25%, third 25%, and fourth 25% of the windowed samples are added together as vectors. It is fairly easy to model this process in Matlab to demonstrate the time aliasing that WOLA does permits a more optimal bin shape than a standard windowed FFT. John
> > john wrote: > > Alberto wrote: > > > I have a few programs that compute and display a running waterfall > > > of an audio signal input to the PC sound card. The waterfall is > > > computed via ordinary FFTs up to 1048576 points for the highest > > > resolution settings. > > > Recently I have seen a similar program that uses polyphase FFTs > > > and it has a spectral leakage much lower than mine (I use a Blackman > > > window). The subject of polyphase FFT is new to me, and from a cursory > > > investigation with Google, it looks like (but I am certainly wrong) > > > the signal is passed through a bank of FIR filters whose width is > > > comparable to a single FFT bin. This looks almost impossible to me, > > > as the time needed to compute 1048576 such FIR filters would make > > > the program not suitable for real time operation. > > > > > > Is there anybody willing to clarify a little this subject, maybe > > > with hints or examples on how to implement in C code a polyphase > > > FFT whose goal is to compute magnitudes to be used for the waterfall ? > > > > > > Thanks indeed > > > > > > Alberto > > > > I'm not familiar with the term polyphase FFT, but I'll suggest a > > technique that can reduce leakage beyond what is possible using an > > ordinary window. It is called Weighted Overlap Add (WOLA) and is > > described in an out-of-print book called Multirate DSP by Crochiere and > > Rabiner. For an FFT size of M, a 4:1 WOLA is computed by windowing off > > a chunk of 4M samples, then stacking and adding each length M segment > > of the windowed data together to make a length M time-aliased vector -- > > your FFT input. 50% overlap is used by sliding 2M new samples into the > > next window calculation. > > > > John
Hi Alberto

> Recently I have seen a similar program that uses polyphase FFTs > and it has a spectral leakage much lower than mine
Is this by any chance http://www.dxatlas.com/Rocky/? If so my Googling of Polyphase FFT revealed no results explaining the method either. John's explanation would intuitively produce the right results, and I am indebted to him for his explanation that it is likely to be simply a further enhancement to OLA techniques. I have seen this extension to the classic 50% OLA in a text elsewhere but I don't remember reading an explanation of the benefits, although it seems clear now. Kind Regards, Howard
Alberto

Just spoke with the Rocky author:

"The basic idea of polyphase FFT is that you take an FFT that is longer than
required, and then decimate it - except that calculations are arranged in
such a way that you save some CPU cycles by not computing the values that
will be discarded during decimation. The algorithm is described in detail,
though under a different name, in this article:
http://archive.chipcenter.com/dsp/DSP000315F1.html
"My implementation is different, in order to improve the program speed I
re-use the FFT already calculated for another purpose and just convert its
output to the polyphase form."

Kind Regards, Howard


Howard,

  many thanks indeed for your message. Yes, the very nice Rocky program is
what I had in mind when writing my first posting.
I will read that reference. But now I am just curious about the WOLA method,
which seems to be of easy implementation. Just one detail is not completely
defined (for me) : the weighting (windowing) function used in the first step
after concatenating 4 length M input samples, is just one of the customary 
windowing functions, or must it be a special one ? And many thanks to John
also for having introduced me to the WOLA algorithm.

BTW in my search with Google of the Polyphase FFT I found this article which
can be of some interest :
http://www.wsdmag.com/Articles/Index.cfm?ArticleID=7108&pg=2

Alberto
-----------------------------------------------

Howard Long wrote:
> Alberto > > Just spoke with the Rocky author: > > "The basic idea of polyphase FFT is that you take an FFT that is longer than > required, and then decimate it - except that calculations are arranged in > such a way that you save some CPU cycles by not computing the values that > will be discarded during decimation. The algorithm is described in detail, > though under a different name, in this article: > http://archive.chipcenter.com/dsp/DSP000315F1.html > "My implementation is different, in order to improve the program speed I > re-use the FFT already calculated for another purpose and just convert its > output to the polyphase form." > > Kind Regards, Howard > >
Howard,

  I have just read the reference cited in your message, and I must really
thank you, as it has made all the subject much, much more clearer.
I think to have now the needed ammunitions to code my Polyphase FFT, or
should I call it window-presum FFT ? :-)

Alberto
----------------------------------------

Howard Long wrote:
> Alberto > > Just spoke with the Rocky author: > > "The basic idea of polyphase FFT is that you take an FFT that is longer than > required, and then decimate it - except that calculations are arranged in > such a way that you save some CPU cycles by not computing the values that > will be discarded during decimation. The algorithm is described in detail, > though under a different name, in this article: > http://archive.chipcenter.com/dsp/DSP000315F1.html > "My implementation is different, in order to improve the program speed I > re-use the FFT already calculated for another purpose and just convert its > output to the polyphase form." > > Kind Regards, Howard > >
Alberto wrote:
> Howard, > > many thanks indeed for your message. Yes, the very nice Rocky program is > what I had in mind when writing my first posting. > I will read that reference. But now I am just curious about the WOLA method, > which seems to be of easy implementation. Just one detail is not completely > defined (for me) : the weighting (windowing) function used in the first step > after concatenating 4 length M input samples, is just one of the customary > windowing functions, or must it be a special one ? And many thanks to John > also for having introduced me to the WOLA algorithm. > > BTW in my search with Google of the Polyphase FFT I found this article which > can be of some interest : > http://www.wsdmag.com/Articles/Index.cfm?ArticleID=7108&pg=2 > > Alberto
The window function used on the 4M sample block depends on the desired bin response. What I have done is design my own window using the 'firls()' routine in Matlab. Since 'firls()' is limited in how big of a result it can produce, just interpolate to get up to 4M. John
"Alberto" <i2phdNOSPAMTHANKS@weaksignals.com> wrote in message
news:Jy4af.19395$65.599136@twister1.libero.it...
> Howard, > > I have just read the reference cited in your message, and I must really > thank you, as it has made all the subject much, much more clearer. > I think to have now the needed ammunitions to code my Polyphase FFT, or > should I call it window-presum FFT ? :-) >
I just got home and reached directly for my Rick Lyons (2nd edn) book and headed straight to the tricks section where I remembered a section on spectrum analysis, and there it is on page 544. Rick refers to it as either WOLA or window-presum depending on the reference. And I guess now it's also known as Polyphase FFT. Just to confuse us mere mortals ;-) Cheers, Howard