[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Discuss-gnuradio] new peakfinding and averaging wxgui fft code (and
From: |
Martin Dvh |
Subject: |
Re: [Discuss-gnuradio] new peakfinding and averaging wxgui fft code (and dummy block has huge impact on usability fir filter) |
Date: |
Wed, 02 Feb 2005 03:49:13 +0100 |
User-agent: |
Mozilla Thunderbird 0.9 (X11/20041124) |
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