x-coordinate of peaks in 1D plot

Hello,

I have the following 1D plot:

http://i46.tinypic.com/35clydu.jpg

I want to find all the x-coordinate values of the peaks in a above plot. I have used the following lines of code and end up having only the y-coordinates(maginutude) values sorted in descending order:

out; %row vector containing the magnitude values(dependent variable)
omega  %row vector containing the frequency axis of FFT (independent variable)

peaks=imregionalmax(out);
[label n]=bwlabel(peaks);
peak = imdilate(peaks,strel('disk',1));
G=regionprops(label,'PixelIdxList');
maxForThisBlob=zeros(1,n);

for blob = 1 : n
  % Get pixel intensities of this particular blob.
  thisBlob = out(G(blob).PixelIdxList);
  % Find the max intensity out of those pixel intensities.
  maxForThisBlob(blob) = max(thisBlob(:));
end
% Sort them in descending order.
sortedMaxes = sort(maxForThisBlob, 'descend');

%----plotting----%
N=size(omega,2); 
figure;
plot(omega(1:N/2),out,'r') %plot one-sided spectrum

The above variable sortedMaxes contains the Magnitude (y-coordinate) values of the peaks(that is 0.02542,0.0191,......so on). How can I extract the corresponding x-coordinates at each of the peaks(that is 3.6,7.2,......so on) ?

please help

Answer 1

[Walter Roberson]

[sortedMaxes, maxidx] = sort(maxForThisBlob, 'descend');
maxx = omega(maxidx);

[zozo]

@walter: the values in 'maxx' do not match with those on the graph..is there offset factor?

[Walter Roberson]

Hmmmm....

Looking again I see that maxidx would be a blob number. I do not have any feel for how blobs correspond to peak positions in your graph.

Change your G=regionprops(label,'PixelIdxList'); to

G=regionprops(label,'PixelIdxList', 'Centroid');

Then

for blob = 1 : n

thiscentroid = G(maxidx(blob)).Centroid;

blobx(blob) = thiscentroid(1);

end

Then blobx will be a vector of x coordinates of the centroids (for lack of anything better to use). The x coordinates would be pixel coordinates.

Then you have the problem of converting pixel coordinates to graph x coordinates. Doing that automatically (without user input) is a challenge. Programs can help, such as http://www.mathworks.co.uk/matlabcentral/fileexchange/7173-grabit

or http://www.mathworks.com/matlabcentral/fileexchange/4316-reverseplot

Though, you could probably automatically detect the pixel coordinates corresponding to the beginning and ending of the x axis without so much difficulty. I guess you could use interp1() to interpolate those to omega values.

XLen = Xaxisend - Xaxisstart + 1; %might need to tweaked as right end of plot does not touch axis

blobomega = interp1( Xaxisstart : Xaxisend, linspace( omega(1), omega(end), XLen ), blobx );

Or more simply,

XLen = Xaxisend - Xaxisstart + 1; %might need to tweaked as right end of plot does not touch axis

blobomega = omega(1) + (blobx - Xstart) / XLen * (omega(end) - omega(1));

[zozo]

@walter: that looks very challenging..!! Is there anyway I can combine the results from sortedMaxes (from my piece of code) and findpeaks() as illustrated below by other experts to achieve my goal?

Answer 2

[Honglei Chen]

You can use findpeaks

http://www.mathworks.com/help/toolbox/signal/ref/findpeaks.html

[zozo]

@Chen: findpeaks is detecting all the peaks in the graph. I need only the ones which are dominant, that is the fundamental frequency and its harmonics as shown in my plot above.

[Matt Kindig]

Use findpeaks to get the frequencies corresponding to the peak amplitudes, and then identify which of these frequencies are harmonics (integer multiples of maxx). You should use a threshold to determine whether each of the peak frequencies is within a suitably small tolerance of the desired harmonic frequency.

[zozo]

i got our idea, but Iam having problems 'filtering' the fundamental freq. and its harmonics from the rest of the peaks detected by findpeaks().

Can u help me with that?

[zozo]

do I have to use the results from sortedMaxes and findpeaks() together to achieve my target? please guide.

[Honglei Chen]

The findpeaks() has an option to sort the peaks. You probably want to set a minimum peak height and then get the 3 maximum peaks.

[Walter Roberson]

findpeaks() applies to datasets, not to images of plots. imdilate() and so on apply to images, not to datasets.

Which is it that you have: the data for the plot, or the image of the plot?

[zozo]

data-set...It is row vector _out(1x10000)_ containing the Magnitude of frequencies.

[zozo]

@Chen: Yes, I have done that using 'MINPEAKHEIGHT' handle and it works. But the point is, those peaks should be automatically detected without me entering the 'minpeakheight' value. The reason for this is I have different signals which produce different plots with different scales of magnitude. So, cant I do it somehow that no matter what the signal, my program should check if it has a peak frequency 'f' and atleast 3 harmonics(its integer multiples), then 'f' should be declared as dominant or fundamental frequency and display its value as output.

[Honglei Chen]

It seems that you can just use 'SORTSTR' to 'descend' and then pick the first 3 peaks.

However, this only works if you are reasonably confident for your data set because as long as you have noise, there will be a lot of small peaks. Hence, if there are no big peaks, the small peaks will be chosen. To solve this issue, you will need some sort of adaptive minimum peak height. You will need to estimate the noise floor of your data and then use that to come up with a reasonable minimum height

[zozo]

okay..thank you @honglei and @walter