GNU Radio, One Step at a Time: [Discuss-gnuradio] Some misconceptions about the "peak

/* -*- c++ -*- */
/*
* Copyright 2007,2010,2013 Free Software Foundation, Inc.
*
* This file is part of GNU Radio
*
* GNU Radio is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3, or (at your option)
* any later version.
*
* GNU Radio is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Radio; see the file COPYING. If not, write to
* the Free Software Foundation, Inc., 51 Franklin Street,
* Boston, MA 02110-1301, USA.
*/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "peak_detector2_fb_impl.h"
#include <gnuradio/io_signature.h>
#include <string.h>

namespace gr {
namespace blocks {

peak_detector2_fb::sptr
peak_detector2_fb::make(float threshold_factor_rise,
int look_ahead, float alpha)
{
return gnuradio::get_initial_sptr
(new peak_detector2_fb_impl(threshold_factor_rise,
look_ahead, alpha));
}

peak_detector2_fb_impl::peak_detector2_fb_impl(float threshold_factor_rise,
int look_ahead, float alpha)
: sync_block("peak_detector2_fb",
io_signature::make(1, 1, sizeof(float)),
io_signature::make2(1, 2, sizeof(char), sizeof(float))),
d_threshold_factor_rise(threshold_factor_rise),
d_look_ahead(look_ahead), d_alpha(alpha), d_avg(0.0f), d_found(false)
{
}

peak_detector2_fb_impl::~peak_detector2_fb_impl()
{
}

int
peak_detector2_fb_impl::work(int noutput_items,
gr_vector_const_void_star &input_items,
gr_vector_void_star &output_items)
{
float *iptr = (float *)input_items[0];
char *optr = (char *)output_items[0];

assert(noutput_items >= 2);

memset(optr, 0, noutput_items*sizeof(char));

printf("\nEnter work(): noutput_items= %d \n", noutput_items);
printf("look ahead remaining=%d\n",d_look_ahead_remaining);

for(int i = 0; i < noutput_items; i++) {
if(!d_found) {
// Have not yet detected presence of peak
if(iptr[i] > d_avg * (1.0f + d_threshold_factor_rise)) {
d_found = true;
d_look_ahead_remaining = d_look_ahead;
d_peak_val = -(float)INFINITY;
printf("crossed threshold i=%d\n",i);
}
else {
d_avg = d_alpha*iptr[i] + (1.0f - d_alpha)*d_avg;
}
}
else {
// Detected presence of peak
if(iptr[i] > d_peak_val) {
d_peak_val = iptr[i];
d_peak_ind = i;
printf("update peak i=%d, val=%f\n",i,d_peak_val);
}
else if(d_look_ahead_remaining <= 0) {
optr[d_peak_ind] = 1;
d_found = false;
d_avg = iptr[i];
printf("PEAK=%f\n",d_peak_val);
}

// Have not yet located peak, loop and keep searching.
d_look_ahead_remaining--;
}

// Every iteration of the loop, write debugging signal out if
// connected:
if(output_items.size() == 2) {
float *sigout = (float *)output_items[1];
sigout[i] = d_avg;
}
} // loop

if(!d_found) {
printf("returning (state=not found) items =%d\n",noutput_items);
return noutput_items;
}

// else if detected presence, keep searching during the next call to work.
int tmp = d_peak_ind;
d_peak_ind = 1;

printf("returning (state=found) items =%d\n",tmp - 1);
if(tmp - 1<0) printf("BUGGGGGGGGGG in peak_detector2!!!!!!!!!!!!\n");

return tmp - 1;
}

} /* namespace blocks */
} /* namespace gr */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "peak_detector2_fb_impl.h"
#include <gnuradio/io_signature.h>
#include <string.h>

namespace gr {
namespace blocks {

peak_detector2_fb_impl::~peak_detector2_fb_impl()
{
}

printf("\nEnter work(): noutput_items= %d \n", noutput_items);

if(output_items.size() == 2)
sigout = (float *)output_items[1];

memset(optr, 0, noutput_items*sizeof(char));

if(d_found==false) { // have not crossed threshold yet
for(int i=0;i<noutput_items;i++) {
d_avg = d_alpha*iptr[i] + (1.0f - d_alpha)*d_avg;
if(output_items.size() == 2)
sigout[i]=d_avg;
if(iptr[i] > d_avg * d_threshold_factor_rise) {
d_found = true;
d_peak_val = -(float)INFINITY;
printf("crossed threshold at i=%d\n",i);
printf("return %d items\n",i);
return i;
}
}
printf("returning (below threshold) items =%d\n",noutput_items);
return noutput_items;
} // end d_found==false
else if(noutput_items>=d_look_ahead) { // can complete in this call
printf("Can complete in this call\n");
for(int i=0;i<d_look_ahead;i++) {
d_avg = d_alpha*iptr[i] + (1.0f - d_alpha)*d_avg;
if(output_items.size() == 2)
sigout[i]=d_avg;
if(iptr[i] > d_peak_val) {
d_peak_val = iptr[i];
d_peak_ind =i;
printf("peak found at i=%d, val=%f\n",i,d_peak_val);
}
}
optr[d_peak_ind] = 1;
printf("PEAK=%f\n",d_peak_val);
d_found = false; // start searching again
printf("returning (above threshold and finished searching) items =%d\n",d_look_ahead);
return d_look_ahead;
} // end can complete in this call
else { // cannot complete in this call
// this is the point of contention:
//
// srictly speaking the block can always get in here and
// it will never return for a finite input sequence.
// This is not an issue, since this block is typically used
// with a continuous stream of inputs so eventually the condition
// noutput_items>=d_look_ahead
// will be met.
//
// However we can choose to bailout from this situation even
// with a finite sequence. This can be done by counting how many times
// we enter here continuously and bailout if a max number of entries
// has been recorder.
// A smoother way to do that is to reduce the d_look_ahead
// parameter by some number every time we are here.
// Observe that this changes somewhat the intended
// behaviour of the block
//
//d_look_ahead--; // this is an example...
//
//Alternatively, a user-defined parameter can be added in the block
//that switches between these two behaviours.

printf("returning (above threshold, but not enough inputs). New lookahead=%d\n",d_look_ahead);
return 0; // ask for more
}
}

} /* namespace blocks */
} /* namespace gr */
#!/usr/bin/env python2
##################################################
# GNU Radio Python Flow Graph
# Title: Test Peak2
# Generated: Thu Apr 9 17:24:11 2015
##################################################

from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser

class test_peak2(gr.top_block):

def __init__(self):
gr.top_block.__init__(self, "Test Peak2")

##################################################
# Variables
##################################################
self.n = n = 100
self.data = data = (0,1,2,3,4,5,6,7,8,9,10,9,8,7,6,5,4,3,2,1,0)+ n*(0,)
self.length = length = len(data)

##################################################
# Blocks
##################################################
self.blocks_vector_source_x_0 = blocks.vector_source_f(data, False, 1, [])
self.blocks_vector_sink_x_0 = blocks.vector_sink_b(1)
self.blocks_peak_detector2_fb_0 = blocks.peak_detector2_fb(5.0, 15, 0.1)

##################################################
# Connections
##################################################
self.connect((self.blocks_vector_source_x_0, 0), (self.blocks_peak_detector2_fb_0, 0))
self.connect((self.blocks_peak_detector2_fb_0, 0), (self.blocks_vector_sink_x_0, 0))

if __name__ == '__main__':
parser = OptionParser(option_class=eng_option, usage="%prog: [options]")
(options, args) = parser.parse_args()
tb = test_peak2()
tb.start()
tb.wait()
#!/usr/bin/env python2
##################################################
# GNU Radio Python Flow Graph
# Title: Test Peak2
# Generated: Thu Apr 9 17:24:11 2015
##################################################

from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser

class test_peak2(gr.top_block):

def __init__(self):
gr.top_block.__init__(self, "Test Peak2")

##################################################
# Variables
##################################################
self.n = n = 100 # put n=10 and run a couple of times to recerate the return -1 bug
self.m = m = 100
self.l = l = 8100
self.data = data = l*(0,)+ (10,)+ m*(0,)+(100,)+ n*(0,)
self.length = length = len(data)

if __name__ == '__main__':
parser = OptionParser(option_class=eng_option, usage="%prog: [options]")
(options, args) = parser.parse_args()
tb = test_peak2()
tb.start()
tb.wait()

Hi all,

recently there has been some discussion regarding the peak_detector2 block, both in the github/gnuradio (pull request 404) as well as in the issue tracker (issue 783).

It is now well accepted that this block is buggy: there are cases the work function returns -1, which is a bug (see issue 783 on how to recreate this bug).

I believe however that there is a DEEPER misconception about how this block works/should work that has resulted in some frustration on what an appropriate fix should be.
In particular there is an insistence that an appropriate bug fix should pass the qa_test of this block and it should be [in the spirit] of the existing algorithm.
In the following I will explain why passing the qa_test is a consequence of the buggy behaviour of this block and NOT its feature.

In addition I will suggest what a proper behaviour of this block should be, so that others who may want to write their own version of a peak detector find it useful.

--------

So the peak_detector block is very reasonable in its conception and its name is very informative and appropriate. It works as follows:

Reads the input and keeps track of a running average (through a single-pole iir filter)

When the current input crosses a threshold (= average * a user-defined factor) upwards the block enters a search state, where it looks for the maximum value of the input over a window of user-defined length.

This is clearly the intended behaviour of the block according to the documentation (I don't know who the original author is...).

----------

It should be obvious that this block is intended for a scenario where the inputs are coming continuously and the block is determining the peaks.
A prototypical application is finding the peak after a MF with a training sequence in order to determine the start of a packet
(BTW, there is some effort on creating a correlate and sync block, which is very relevant to this discussion).

It should also be obvious that both the pole of the IIR filter, the threshold, AND the window should be properly selected per application

(eg, for a training sequence of length N, it is appropriate to set the window ~N) otherwise the block will simply not work according to its specs.

-----------

So here are some questions and my suggested answers:

Q1: How should this block behave when it has already crossed the threshold upwards and is still inside the user-defined search window but there is not enough input data to process?

A1: Clearly the work function should return (consuming all inputs/outputs up to the most current peak--but not more) and hopefully the next time it is called, the scheduler would provide more samples so that the search over the remaining window is exhausted and the peak is reported.

---

Q2: And then, what if the scheduler does not give this block more input data? Maybe because (as in the qa_test) this block is connected to a source that outputs a finite (=21 !) number of samples?

A2: Since the intended use of this block is a continuous stream of data, the above situation is an inappropriate use of this block, and thus it cannot be considered a test case. This means that an appropriate test case should make sure that enough tail data is in the finite input so that the search window is exhausted.

---

Q3: What if the user insists in using the block in a way that is not intended, ie, by providing it with a finite (and small compared to the search window) number of inputs. Should there be a way that the block bails out?

A3: The current implementation has this behaviour but it is a BUG not a feature: ERRONEOUSLY it decrements the window even if it does not return any data, so eventually the window becomes smaller than the available input and the block exits!
This results in the block passing an ill-conceived qa_test (where the window is set to 1000 for a peak occurring within 10 samples of the threshold crossing). However the intended use of the block was never this (according to the documentation).

------------

From the above it should be clear that any fix of this block should NOT

be in the spirit of the original one, nor not passing the existing qa_test is a measure of its difference from the original block (at least from the stated intended behaviour of the original block).

I attach here a clean modification of this block that has the intended behaviour. It DOES NOT pass (correctly) the inappropriate qa_test,

but it does pass more appropriately designed test(s) (also attached).

------------

I am requesting some feedback for this approach so that I can finalize it and submit it for merging into the gnuradio trunk.

best

Achilleas

GNU Radio, One Step at a Time

Friday, April 10, 2015

[Discuss-gnuradio] Some misconceptions about the "peak_detector2" block

No comments:

Post a Comment