How to find the modes of a biphasic histogram?

Hi there, I'm plotting histograms that occasionally turn out to be bimodal. When they are bimodal, I want to find the dual peak points of the histogram so I can characterise them accordingly. Are there any neat ways I can do this? Currently I'm stuck at finding the FWHM for the max peak and then calculating the mid point for each peak for this.

Ideally I would prefer something that characterised the two distributions within the histogram and have the nagging feeling someone will have already written some nice code for this...

Many thanks

[Image Analyst]

When they are bimodal, how well separated are the peaks? Can you upload a screenshot of the bimodal histogram? In the meantime, look up Gaussian Mixture Models http://en.wikipedia.org/wiki/Mixture_model#Gaussian_mixture_model and K-S test http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test

Answer 1

[Star Strider]

To the best of my knowledge, there aren’t any functions to specifically do what you want. In their absence, my approach is to fit a polynomial to the histogram data and then do the ‘Calculus I’ calculations on it to identify the x-axis (or more appropriately ‘bin’-axis) locations of the peaks. It worked on my test data, and lacking yours to test it on, I leave it to you to see if it works on yours:

% Create Data: 
D = 1E+4*exp(-([randi(15,1,100)  randi(15,1,100)+20]/10).^2);
B = 1:0.1:15;
C = histc(D, B);

% Find Peaks: 
prm = polyfit(B, C, 7);             % Polynomial Fit To Histogram
d1prm = polyder(prm);               % First Derivative (D1)
d2prm = polyder(d1prm);             % Second Derivative (D2)
rd1prm = roots(d1prm);              % All Extrema Locations
xtr = polyval(d2prm,rd1prm);        % Evaluate D2 At All Extrema Locations
pks = rd1prm(xtr < 0);              % Peaks Are Where D2 < 0

Note that here ‘C’ are the histogram data and ‘B’ are the bin locations. My ‘pks’ variable should return the locations of the peaks in the histogram. Experiment with it to get the results you want with your data (the order of the polynomial likely being the most relevant).

If you have the Signal Processing Toolbox, the findpeaks function is quite likely more robust, and I suggest you try it first if it’s available to you.