[SciPy-Dev] Function for finding relative maxima of a 1d array
Tue Sep 20 15:13:27 CDT 2011
On Tue, Sep 20, 2011 at 7:25 PM, Jacob Silterra <email@example.com> wrote:
> Hello all,
> I posted an idea for a simple function for finding relative maxima of 1D
> data in the numpy mailing list; the general consensus seemed to be that it
> was useless but more advanced algorithm would be useful in scipy. Since
> this is an addition of new functionality, I thought (re: was told in the
> case of numpy) it'd be best to discuss it over the mailing list before
> issuing a pull request. Full code available at
> https://github.com/jsilter/scipy/tree/find_peaks. The method is named
> find_peaks, located in optimize.py in the optimize module.
That looks promising.
> Copied from the docstring:
> The algorithm is as follows:
> 1. Perform a continuous wavelet transform on `vector`, for the supplied
> `widths`. This is a convolution of `vector` with `wavelet(width)` for
> each width in `widths`. See scipy.signals.cwt (pending).
There is also a CWT implementation at
http://projects.scipy.org/scipy/ticket/922 (the first code review link still
works). It didn't get merged but still contains examples and references that
may be useful. Also the author may still be around and willing to review
your code for example.
> 2. Identify "ridge lines" in the cwt matrix. These are relative maxima
> at each row, connected across adjacent rows. See identify_ridge_lines
> (in optimize.py)
> 3. Filter the ridge_lines using filter_ridge_lines (in optimize.py) by
> length and signal to noise ratio.
> This algorithm requires several parameters be specified, I chose what I
> thought were reasonable values to be defaults. This method was designed for
> use in mass-spectrogram analysis, which typically have large,sharp peaks on
> top of a flat and noisy baseline. With proper parameter selection, I believe
> it would function well for a wide class of peaks. The default wavelet for
> step 1 is a ricker <http://en.wikipedia.org/wiki/Mexican_hat_wavelet>wavelet, which is added in scipy.signal.wavelets.
> This is implemented in pure python, which is appropriate for steps 1 and 3.
> Step 2 involves a lot of looping, that may be better in C, but it's Python
> at the moment.
Optimizing can always be done at a later stage, no need to worry about that
> I used my best judgement on placing functions in appropriate modules, which
> led to some circular dependencies. As a temporary fix, the imports are
> within each function rather than at the top of each module. Once it's
> decided where each function should actually live that can be resolved.
I think ricker landed in the right place, cwt could probably also live in
wavelets.py. optimize.py is not the best place for peak finding algorithms,
it's already a very large file. Why not create a _peak_finding.py file under
signal? That should get rid of those circular imports.
> PS There's a good comparison different algorithms available here<http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631518/>
> PPS I have Eclipse configured to use PEP8 standards for spacing after
> commas, assignment operators, and so forth. This ended up changing a lot of
> code which may make this specific function hard to review (sorry). I put the
> new functions close to the bottom of each file.
> Code cleanups are very welcome, your changes mostly look good to me. It
would be helpful to have them as a separate commit though. The other
additions can also be split up more to make them more easy to review. For
example one commit for ricker and cwt with tests, one for argrelmin/max and
one for find_peaks and the ridge_lines functions.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the SciPy-Dev