Re: [Discuss-gnuradio] new peakfinding and averaging wxgui fft code (and dummy block has huge impact on usability fir filter)

[Martin Dvh]

James Cooley wrote:
> Speaking of averaging and peak-detection, I'm interested in doing so 
> outside of the gui.
>
> Could you build an averager like this, with fir filters applied to the 
> fft of the signal? Would it work? Is there a better way?
I thought about this too, it would probably perform better also.
> -take fft of signal (much as fft_sink does)
>
> -Run fft output through an averager to smooth it out a bit:
> FIR filter with taps [1/N,1/N,1/N,1/N,1/N,1/N...] up until N number of 
> samples you want to average over. (A running average FIR filter)
> averager = gr.fir_filter_fff(1, taps)
You have to take into account the way the fft-output-data is organized.
You want to average over several fft's.
The current algorithm in the gui does the following
absfft1=abs(fft1)
absfft2=abs(fft2)
...
absfftN=abs(fftN)

averagefft[1]=(absfft1[1]+absfft2[1]+absfft3[1]+.....+absfftN[1])/N
averagefft[2]=(absfft2[1]+absfft2[2]+absfft3[2]+.....+absfftN[2])/N
....
averagefft[FFTSIZE]=(absfft1[FFTSIZE]+absfft2[FFTSIZE]+absfft3[FFTSIZE]+......+absfftN[FFTSIZE])/N

So If you have a stream where all the ffts come in one after another you have 
to consider them as a parallel source.
The current FIR filters can not (yet) handle anything else as complex or float.

You could reorganise the single parallel datastream into a total number of 
FFTSIZE seperate streams and do a FIR average on all of them
after that recombine the output values of the seperate streams into a single 
fft.

Note also that the abs is not optional. As the fft produces complex values which have a magnitude and an arg (angle, phase) you cannot just add 
them.
You can add the magnitudes (hence the abs) or you have to convert them to all 
having the same phase.
The fftw site has info about doing this (precalculate some complex numbers which you use to convert the phases of the seperate output ffts and 
then combine them)

http://www.fftw.org/pruned.html
snippet:
"In this case, one simply computes N/K FFTs of size K, with stride N/K in the data, and then combines the output by summing with appropriate 
phase ("twiddle") factors. Essentially, you are performing "by hand" a single radix-K decimation-in-frequency Cooley-Tukey step, where you only 
compute one size-K block of the output."

If you take this path you get an actual averaged complex fft in stead of only 
the magnitude.

Ofcourse you can also compute a veryvery big FFT and only use the upper frequencies. This is in fact the same result as averagiing several FFTS, 
but it  is much more computational intensive.


Here some pseudocode which could be used/built for calculating an abs averaged 
fft using fir:

parallelize=gr.serial_to_parallel(FFTSIZE)
ffts=gr.gr_fft_vcc(FFTSIZE)
parallelmag =gr.parallel_complex_to_mag()   //We need a version which understands the parallel stream and just gets the abs value of every 
complex value in
seperator=gr.seperate_streams(FFTSIZE) //We need something to split up the ffts
firaverager1=gr_fir_filter_fff(averaging_taps)
firaverager2=gr_fir_filter_fff(averaging_taps)
firaverager3=gr_fir_filter_fff(averaging_taps)
....
firaveragerFFTSIZE=gr_fir_filter_fff(averaging_taps)

combiner=gr.combiner(FFTSIZE)  //we need something that recombines the streams

fg.connect(src,parallelize)
fg.connect(paralellize,ffts)
fg.connect(ffts,parallelmag)
fg.connect(parralelmag,seperator)
fg.connect([seperator,0],firaverager0)
fg.connect([seperator,1],firaverager1)
....
fg.connect([seperator,FFTSIZE],firaveragerFFTSIZE)

fg.connect(firaverager0,[combiner,0])
fg.connect(firaverager1,[combiner,1])
.....
fg.connect(firaveragerFFTSIZE,[combiner,FFTSIZE])

the combiner output now has the averaged fft.
The trick to use the complex values and using phase ("twiddle") factors could 
also put into this in stead of the parallel_complex_to_mag block.


This is all pseudocode as several of the blocks mentioned are not there (yet)
If we do all of this in one block you don't need this many new blocks.
But you still need 1024 fir filters for a 1024 FFT, but the rate is also a 
factor 1024 lower then the source.
The current gui fft code also does a decimation which is perfectly ok, only 
your output rate gets lower.


Another approach is just to forget using fir for the averaging and just stick 
with the python code in wxgui and make a non-gui block.
For me it is fast enough. (I average over 100 FFTs, do a peakfind and still get 
a few frames per second)
Or make/modify a new  fir averaging block so it can handle more complicated 
(parallel) streams.


> -run filter to find some peaks. Can this be done by constructing a FIR 
> with taps [1,-1] and then noting where the 0 crossings take place? These 
> taps should be equivalent to:
> y[n] = -1*x[n] + 1*x[n+1]
> which is like the slope.... derivative in continuous functions.... over 
> a pair of samples.
You really also want to look at the neighbouring values. Even an averaged fft has many many zero crossings in its derivative. You only want the 
 zerocrossings which have no nearby neigbours which are higher then itsself. (You want global and semi-global maxima, not all local maxima I 
would guess)
Even better would be to do a real gaussian peak fit, but don't expect 
interactive rates with that.
The nice thing about a gaussian peak fit would be that you also get an 
indication of the bandwith of the signal.

The peakfitting algorithm in the wxgui code is in its own file so it should be 
possible to use it in a non-gui block.
It it written in python but can be rewritten in C++ or anything you like. (The 
original was in java)


> -probably you'd want to also use a threshold on the original fft signal 
> to weed out insignificant peaks?
Yes, this is in the current peakfinding also used.

Greetings,
Martin


> Thanks,
>
> -jamie
>
>
>
>
> Martin Dvh wrote:
>
> > Eric I am Recompiling now (again;-).
> > I don't know if I will be doing the testing tonight (It is kind of 
> > late here in the Netherlands)
> >
> > I did however also finish my modified fft_sink in wxgui and put in on 
> > my website.
> > I added the (default) option to average over several fft frames 
> > (default 100)
> > This gives a much quiter spectrum
> >
> > The most interesting feature I also added is a peakfinding algorithm.
> > The spectrumwindow can now find peaks and shows the exact frequency 
> > and signal level of these peaks.
> >
> > def make_fft_sink_f (fg, parent, label, fft_size, 
> > input_rate,show_peak_info=False,show_peak_markers=False,navg_frames=100,peak_searchwidth=-1): 
> >
> >     #show_peak_info (default=False): if set to true the fft_window 
> > will show the frequency and level of every big peak found in the spectrum
> >     #show_peak_markers (default=False):  if set to true the fft_window 
> > will show a marker at every big peak found in the spectrum
> >     #navg_frames (default=100): number of fft frames to average before 
> > showing it. If set to 100 or higher you get a much quiter spectrum 
> > (less noisy and not changing very much every frame)
> >     #peak_search_width(-1 gets you the default of 10 for a fft_size of 
> > 1024): number of neigbours the peakfinding algorithm will look at for 
> > determing if something is a peak.
> >     #                                         A higher value means you 
> > get less peaks detected
> >     #Default the peakfinding algorithm determines a minimum 
> > signalstrength threshold below which no peaks are detected
> >     # you can override this with set_autopeakthreshold(False) and 
> > set_peakthreshold(mythreshold) to use your own threshold
> >     # the higher the threshold, the less peaks found
> >
> >
> > I modified Makefile.am,fftsink.py and added peakfind.py
> > You can all find them on:
> > http://www.olifantasia.com/projects/gnuradio/mdvh/wxgui/
> >
> > I ported peakfind.py fromjava code from:
> > #ported from  YALE - Yet Another Learning Environment
> > #  Copyright (C) 2001-2004
> > #      Simon Fischer, Ralf Klinkenberg, Ingo Mierswa,
> > #          Katharina Morik, Oliver Ritthoff
> > #      Artificial Intelligence Unit
> > #      Computer Science Department
> > #      University of Dortmund
> > #      44221 Dortmund,  Germany
> > #  email: address@hidden
> > #  web:   http://yale.cs.uni-dortmund.de/java
> > Which is GPL
> >
> > The files need to be put into:
> > gr-wxgui/src/python/gnuradio/wxgui
> >
> > I also tried to make a diff against cvs.
> > It is at 
> > http://www.olifantasia.com/projects/gnuradio/mdvh/wxgui/fft_avg_and_peakfind.diff 
> >
> > But cvs diff does not seem to add any new local files. So you have to 
> > still add peakfind.py manually.
> >
> > How do I tell cvs diff to include local new files ???
> >
> > greetings,
> > Martin
> >
> >
> >
> > _______________________________________________
> > Discuss-gnuradio mailing list
> > address@hidden
> > http://lists.gnu.org/mailman/listinfo/discuss-gnuradio