4

What is the best method for finding impulse peaks (dirac delta) in a 2d matrix.

More specifically, I would like to find the harmonic frequencies of a given image and so I need to find impulse peaks in the image absolute value DFT.

I thought of using findpeaks but there's no 2d version. I also saw earlier posts regarding finding ordinary peaks using imdilate and/or imextendedmax but those find all the peaks in a 2d matrix whereas I am only interested in impulse peaks. I am sure DSP people have a common recipe for this...

Please Help,

Thanks

smichak
  • 4,716
  • 3
  • 35
  • 47
  • I'm looking for frequencies in which the Fourier Transform absolute value takes a value of a Dirac delta function at that frequency (infinity). Obviously, for a DFT you don't expect to see infinity but a major peak. Problem is how to find those peaks in 2d and how to distinguish them from normal (lower, non impulse) peaks. – smichak Nov 22 '10 at 05:28

3 Answers3

6

What you want to do is find peaks with high contrast. Thus, you need a way to identify local maxima, plus a way to measure the difference between the peak and the surrounding values. Thresholding on this difference will identify the impulse peaks for you.

Assuming your input signal is called signal

%# dilate to find, for every pixel, the maximum of its neighbors
dilationMask = ones(3);
dilationMask(5) = 0;
dilSignal = imdilate(signal, dilationMask);

%# find all peaks
%# peaks = signal > dilSignal;

%# find large peaks peaks by thresholding, i.e. you accept a peak only 
%# if it's more than 'threshold' higher than its neighbors
peaks = (signal - dilSignal) > threshold;

peaks is a logical array with 1's wherever there is a good peak. You can use it to read peak heights from signal with signal(peaks), and to find coordinates using find(peaks).

Jonas
  • 74,690
  • 10
  • 137
  • 177
  • OK - that's a start... but how do I choose threshold? I assume it should depend on the input signal. Is there a good method /rule of thumb for choosing it? – smichak Nov 22 '10 at 20:50
  • They're Dirac peaks, so the threshold is infinity :P. Joking aside, there is no single best way of choosing a threshold. You can try and run `graythresh` on `signal-dilSignal` if there are two clear peaks, one for signal and one for noise. You can also set the threshold as 3 robust standard deviations (median absolute deviation) of the signal. Maybe you can even decide on a fixed value if all your data is very consistent, and the height of the peak has an easily identifiable meaning. – Jonas Nov 22 '10 at 22:26
  • OK - that sounds promising... I will go with that. – smichak Nov 23 '10 at 05:53
2

This paper I wrote contains Matlab source code for fast local peak detection in 2D. It works similar to imregionalmax() in Mathworks Image Processing Toolbox but allows you to specify a local neighborhood radius: bigger radius -> sparser peaks.

Since you expect sparse impulses, the nonmaxsupp_scanline() function may be suitable for you.

Jonas Heidelberg
  • 4,984
  • 1
  • 27
  • 41
Tuan
  • 21
  • 1
0

The findpeaks algorithm is pretty trivial; if an element is bigger than both its neighbours, then it is a peak. Writing a 2D version of this should be pretty simple.

George Hilliard
  • 15,402
  • 9
  • 58
  • 96
Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680