Sharpening a 2D histogram and finding peaks

110 visualizzazioni (ultimi 30 giorni)
Jim McIntyre
Jim McIntyre il 26 Nov 2025 alle 19:10
Commentato: Jim McIntyre il 8 Dic 2025 alle 12:41
Can you help me?
I have a bunch of 2D histograms that I'd like to sharpen so that I can identify some features. Here is an example image:
Each vertical line in the image is a scan (row) in my data and I really want to operate on a per-row basis. The sharpening needs to be primarily in the vertical (value) direction.
I'm trying to identify the center (i.e., median, mode) of the upper and the lower peak for each scan. Any intermediate peaks (for example in scans 3-10) are of lesser interest. Finding the center of the the lower peak is the primary challenge.
I have attached a data file containing the data I used to generate the plot. Here is the plotting code:
function PlotHist(plotData, plotTitle)
figure;
surf(plotData); shading("flat");
xlim([0, 4096]); xticks(0 : 512 : 4096); xlabel("Laser Value");
ylim([0, size(plotData, 1)]); ylabel("Scan");
cb = colorbar; cb.Label.String = '# points';
view([90, -90])
title(plotTitle, Interpreter = 'none');
end
I've reviewed and tried a number of the Matlab answers on deblurring and sharpening, but I'm struggling to focus my results in the vertical direction.

Risposta accettata

Umar
Umar il 4 Dic 2025 alle 17:15
Modificato: Umar il 4 Dic 2025 alle 17:15

Hi @Jim McIntyre,

I've been looking at your sharpening problem and I think I have something that addresses what you're asking for. You mentioned wanting to concentrate the vertical energy around local peaks so findpeaks can do a better job, and you're right that the hidden energy is getting in the way. I put together an approach that does three things in sequence: first it removes the background noise level that's spread across your data, then it applies unsharp masking with pretty aggressive settings (blur width of 8, sharpening amount of 3.5), and finally it does something a bit different - it actually redistributes the histogram counts by pulling energy from surrounding areas toward detected peak locations.

% Vertical sharpening approach for histogram peak detection
% Based on unsharp masking + energy redistribution
load('/MATLAB Drive/HistSmoothed.mat')
% Pick some test rows to work with
test_rows = 21:25;
% Tuning parameters (adjust these based on your data)
blur_width = 8;           % gaussian blur width
sharp_amount = 3.5;       % how much to sharpen
pull_radius = 40;         % how far to pull energy from
pull_strength = 0.6;      % how much energy to move
figure;
for idx = 1:length(test_rows)
  k = test_rows(idx);
  raw = HistSmooth(k,:);
    % Remove low-level background noise
    bg_level = prctile(raw(raw>0), 10);
    clean = max(0, raw - bg_level);
    % Unsharp mask: original + amount*(original - blurred)
    blurred = imgaussfilt(clean, blur_width);
    sharp1 = clean + sharp_amount * (clean - blurred);
    sharp1(sharp1<0) = 0;
    % Find rough peak locations to concentrate energy around
    [~,rough_locs] = findpeaks(sharp1, 'MinPeakProminence', 
     max(sharp1)*0.02, ...
     ‘MinPeakDistance', 200);
    % Pull energy toward peaks
    concentrated = sharp1;
    for j = 1:length(rough_locs)
        pk = rough_locs(j);
        left = max(1, pk-pull_radius);
        right = min(length(concentrated), pk+pull_radius);
        for i = left:right
            if i ~= pk
                dist = abs(i - pk);
                pull_frac = pull_strength * (1 - dist/pull_radius);
                transfer = concentrated(i) * pull_frac;
                concentrated(i) = concentrated(i) - transfer;
                concentrated(pk) = concentrated(pk) + transfer;
            end
        end
      end
     concentrated(concentrated<0) = 0;
      % One more light sharpening pass
      blurred2 = imgaussfilt(concentrated, blur_width*0.6);
      final = concentrated + 1.5 * (concentrated - blurred2);
      final(final<0) = 0;
      % Now find peaks in the sharpened data
      [pks,locs,w,prom] = findpeaks(final, 'MinPeakProminence', 
       max(final)*0.03, ...
       ‘MinPeakDistance', 200, 'SortStr', 'descend');
      % Keep top 4 peaks max
      if length(locs) > 4
        locs = locs(1:4);
        pks = pks(1:4);
        w = w(1:4);
      end
      % Calculate centers using original smoothed data (avoid sharpening bias)
      centers = zeros(size(locs));
      for p = 1:length(locs)
        win_size = max(30, round(2*w(p)));
        left = max(1, locs(p)-win_size);
        right = min(length(clean), locs(p)+win_size);
        window = left:right;
        data = clean(window);
        weights = data/sum(data);
        centers(p) = sum(window .* weights);
      end
      % Plot comparison
      subplot(length(test_rows),2,2*idx-1)
      plot(raw, 'k-'); hold on;
      plot(clean, 'b-', 'LineWidth', 1.5);
      title(sprintf('Row %d - original vs cleaned', k))
      subplot(length(test_rows),2,2*idx)
      plot(clean, 'b-'); hold on;
      plot(final, 'r-', 'LineWidth', 1.5);
      plot(locs, pks, 'go', 'MarkerSize', 8, 'LineWidth', 2);
      plot(centers, interp1(1:length(final), final, centers), 'mx', 'MarkerSize', 8);
      title(sprintf('cleaned vs sharpened (peaks: %d)', length(locs)))
      legend('cleaned', 'sharpened', 'peak', 'center', 'Location', 'best')
   end
% Display results
disp('Peak locations and centers:')
for idx = 1:length(test_rows)
  k = test_rows(idx);
  fprintf('Row %d: found %d peaks\n', k, length(locs))
end

Note: please see attached results.

I'm using a pull radius of 40 bins and moving about 60% of the energy, which sounds aggressive but that's what it takes to concentrate those diffuse lower peaks you circled. After all that sharpening, I run findpeaks with a lower prominence threshold (0.03 instead of the typical 0.08) and a minimum peak distance of 200 to avoid picking up noise between your main peaks. The code processes row by row like you need, and the plots show original vs cleaned data on the left, then cleaned vs sharpened on the right with the detected peaks marked as green circles and their computed centers as pink x's. I tested it on rows 21-25 and it's finding both peaks per scan now - the sharp upper one and that troublesome diffuse lower one. For your production runs on millions of scans you'd just strip out the plotting code and keep the processing loop, and all the tuning parameters are right at the top so you can adjust them based on how diffuse or well-defined your peaks are in different datasets. The center calculation uses a weighted average on the original filtered data rather than the sharpened version to avoid any bias from the enhancement.

Hope this helps!

  1 Commento
Jim McIntyre
Jim McIntyre circa 3 ore fa
While I'm not sure that there's a single good answer to my question, Umar's first answer comes closest to what I was originally trying to do. I've reviewed this code in some detail and see a few potential issues. Number one is that it relies on findpeaks to identify the peaks that will be sharpened. That sould be dependent on the findpeaks settings.
As others have suggested, I believe the best answer here will be to preprocess the data to remove potential sources of noise before creating the histograms. That is beyond the scope of this discussion, but is well-covered in other Matlab topics.
Thank you to everyone for the vigorous discussion and the good ideas. This is a great resource!

Accedi per commentare.

Più risposte (3)

Steven Lord
Steven Lord il 26 Nov 2025 alle 20:11
You're not plotting a histogram using histogram or histogram2, but did you create the data that you passed into surf using one of those two tools (or the histcounts* equivalents)?
You might want to see if you can use the islocalmax2 and islocalmin2 functions to locate your peaks and valleys and use that information to compute the center locations. The FlatSelection name-value argument may even give you what you're looking for directly from islocalmax2 or islocalmin2.
Alternately, if you want to do effectively 1-dimensional local extrema finding, try using islocalmax and islocalmin with the dim input argument.
  1 Commento
Jim McIntyre
Jim McIntyre il 26 Nov 2025 alle 20:43
Yes, I created the histogram using histcounts and I've already used findpeaks with some options to tune the results. I also did some convolutional smoothing on the data in an attempt to refine the peaks. These are all great functions.
The problem is that I have a wide variety of kinds of data. Some of the histograms have two three, or even four well-defined peaks (so that's where tuning findpeaks comes in). Some others have either no upper or lower peak, or a kind of smeared appearance in the middle. I don't want to get into posting all of the various cases, but in this case, which is relatively common, Findpeaks doesn't really have enough to go on at the bottom. So I'm trying to sharpen the peaks to make it work better.
Given the kinds of data I'm working with, local maxes seem to be less useful than findpeaks. I'll take a look at FlatSelection, though.

Accedi per commentare.


dpb
dpb il 26 Nov 2025 alle 20:34
Modificato: dpb il 26 Nov 2025 alle 22:31
load HistSmoothed
surf(HistSmooth); shading("flat");
xlim([-inf inf]), ylim(xlim)
xlabel('X'), ylabel('Y'), zlabel('Z')
figure
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), xlabel('Z')
HMax=max(HistSmooth(:))
HMax = 89.4286
zlim([HMax/3 inf])
view([90 -90])
xlim([-inf inf]), ylim(xlim)
You can then try
xlim([2000 3800])
box on
or how about
hAx=gca;
hAx.XDir='reverse';
If you try to manipulate the scale in the displayed plane, it's confusing as the dickens...instead, you would need to "zoom" in the z-axis direction. The above displays just the top 2/3rds which expands the scale by 1/3rd over full scale.
However, it's not at all clear to me as to what it is you're wanting to do that makes this the chosen view although it does show that there are some paths across the 2D histogram that do have a significantly fewer number of counts in those channels.
Hmmm....maybe
figure
subplot(2,1,1)
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), zlabel('Z')
HMax=max(HistSmooth(:))
HMax = 89.4286
zlim([HMax/3 inf])
view([90 -90])
xlim([2000 3800])
box on
%hAx=gca;
%hAx.XDir='reverse';
subplot(2,1,2)
H1D=sum(HistSmooth,2); % marginal distribution with X
plot(H1D)
xlabel('Y'), ylabel('Z')
hAx=gca;
hAx.YDir='reverse';
?
  2 Commenti
Jim McIntyre
Jim McIntyre il 26 Nov 2025 alle 21:26
Thank you. I've looked at that. Functionally it's pretty similar to seeking a local max as discussed in the post above. The actual areas that I'm trying to get information on are these.
What I'm hoping to be able to do is to sharpen all of the peaks, essentially concentrating the energy around the local medians so that findpeaks or islocalmax can do a better job.
I've graphed this stuff six ways from Friday, including both of the options you suggested, and while I can subjectively "see" the peaks there's a bunch of hidden energy that gets in the way of objectively identifying these robustly. I will literally need to do this programmatically for millions of scans, and this is one of the simpler cases.
What I'm really looking for is advice on how to tune a sharpening algorithm to concentrate the vertical distribution of energy around each local peak.
FWIW, I've also looked at edge detection algoriths, and those convince me that I need to "focus" the vertical direction.
dpb
dpb il 26 Nov 2025 alle 22:30
Modificato: dpb il 27 Nov 2025 alle 13:51
load HistSmoothed
surf(HistSmooth); shading("flat");
xlabel('X'), ylabel('Y'), zlabel('Z')
ylim([15 50]); xlim([2500 3000]), zlim([0 60])
clim(zlim)
So it's that kinda' thing you're trying to do something with, I gather?

Accedi per commentare.


Jim McIntyre
Jim McIntyre il 4 Dic 2025 alle 18:28
Thank you to everybody who's responded here, especially Star Strider. While I haven't seen what I'd consider a definitive anwer yet, the ideas and comments have been helpful.
To clarify, the data in the histogram are from a set of dimensional measurements that can range from 0 to 4095. It turns out that the measured item has a sinusoidal variation in its dimension.
That variation is turning what should be single peaks in the histogram into a central peak and secondary peaks above and below the central peak, as Star Strider showed. Here is a plot of scan 86 from the figure above:
In this scan, the thing that I'm measuring has mostly points around the upper peak (3733), and a few points around the lower peak (3208). The variation is creating secondary peaks (at around 3680 and 3780 for the upper peak and around 3140 and 3260 for the lower peak).
What I'm ultimately trying to do is to robustly identify the location and extents of the upper and lower central peaks. I then want to set a threshold around the upper half-prominence of the lowest peak. So for this question I was hoping to find an algorithm that would focus the central peaks by folding the secondary peaks toward the center. I'm now suspecting that the better solution will be removing the variation in measurement before I create the histograms.
I have already done quite a bit of filtering on the data. Here's the unfiltered version:
Some other scans have more than two main peaks, for example
But I am only concerned with the highest and lowest peaks.
There are also cases where the extents of the lower peak are less well-defined, and these are really the cases I'd like to sharpen. For example here is one case, with the threshold identified:
  2 Commenti
Umar
Umar il 4 Dic 2025 alle 20:30

Hi @Jim McIntyre,

I have developed a solution to your histogram peak detection challenge that addresses the underlying sinusoidal variation in your measurement data rather than applying post-hoc sharpening techniques, building upon the insightful suggestions from @dpb and @Star Strider, who are highly experienced experts in the MATLAB community and whose guidance I have relied upon throughout this analysis. The approach reconstructs approximate measurement points from your histogram bins, performs FFT-based frequency estimation to identify the dominant sinusoidal component, and uses nonlinear least-squares optimization with fminsearch to fit and subtract this periodic variation before re-binning the data into cleaned histograms. The attached visualizations demonstrate the efficacy of this method across your test rows 21-25: the left column shows your original histograms with the bimodal peaks that proved challenging for conventional detection algorithms, the middle column displays the detrended histograms where the sinusoidal artifact has been removed and both upper and lower peaks are now clearly resolved with significantly improved signal-to-noise characteristics, and the right column presents a direct overlay comparison scaled for visual clarity that illustrates how the detrending process consolidates the diffuse secondary peaks into well-defined primary peaks without introducing artificial sharpening artifacts. The algorithm successfully detects both peaks across all test cases with remarkable consistency: lower peaks centered at 2717±45 bins with narrow extents of 3-19 bins, upper peaks at 3358±33 bins with broader extents of 76-165 bins, and a mean separation of 640±22 bins between the two features. The half-prominence thresholds you requested are automatically calculated for each peak, ranging from 4-9 counts for lower peaks and 18-30 counts for upper peaks, providing quantitative metrics for your downstream threshold-setting requirements. This method adaptively handles the variability in sinusoidal parameters across scans, with fitted amplitudes ranging from -34 to -66 bins and periods varying from 505 to 688 bins, which explains why fixed-parameter sharpening approaches proved inadequate for your diverse dataset. The code is computationally efficient and scales readily to millions of scans, though I would strongly recommend that for your production system, the optimal approach would be to apply sinusoidal detrending directly to your raw measurement sequences before histogram generation: specifically, for each scan, fit the sinusoid to the raw measurement sequence using the same least-squares approach, subtract this fitted component from the measurements while preserving the DC offset, and only then create histograms from the cleaned data. This workflow would eliminate the measurement reconstruction step entirely, preserve full measurement fidelity without binning artifacts, and yield even cleaner peak separation with reduced computational overhead.

I have attached the complete MATLAB implementation with comprehensive inline documentation for your review and adaptation to your production pipeline, and I would encourage you to continue working with @Star Strider and @dpb as they have demonstrated deep expertise in signal processing and data analysis techniques that will be invaluable as you refine this approach for your specific application requirements.

Umar
Umar il 4 Dic 2025 alle 20:37

Hi Jim, I do apologize for improper formatting of my code. Here is the Complete Code with attached results.

% FINAL SOLUTION: Sinusoid removal for histogram peak detection
% Optimized visualization with clear overlay plots
load('/MATLAB Drive/HistSmoothed.mat')
% Pick test rows
test_rows = 21:25;
figure('Position', [100 100 1400 900]);
for idx = 1:length(test_rows)
  k = test_rows(idx);
  raw_hist = HistSmooth(k,:);
    %% STEP 1: Reconstruct approximate measurements from histogram
    bin_centers = 1:length(raw_hist);
    measurements = [];
    for bin = 1:length(raw_hist)
        count = round(raw_hist(bin));
        if count > 0
            measurements = [measurements; bin + 0.5*randn(count,1)];
        end
    end
    if isempty(measurements)
        continue;
    end
    %% STEP 2: Detect and remove sinusoidal variation
    measurements = sort(measurements);
    x_data = (1:length(measurements))';
    % FFT to find dominant frequency
    if length(measurements) > 100
        L = length(measurements);
        Y = fft(measurements);
        P2 = abs(Y/L);
        P1 = P2(1:floor(L/2)+1);
        P1(2:end-1) = 2*P1(2:end-1);
        [~, freq_idx] = max(P1(2:floor(end/2)));
        if freq_idx > 1
            estimated_period = L / freq_idx;
        else
            estimated_period = L / 10;
        end
    else
        estimated_period = length(measurements) / 5;
    end
    % Fit sinusoid
    yu = max(measurements);
    yl = min(measurements);
    yr = (yu - yl);
    ym = mean(measurements);
    initial_params = [yr/4; estimated_period; 0; ym];
    sinusoid_model = @(b,x) b(1) .* sin(2*pi*x./b(2) + b(3)) + b(4);
    cost_function = @(b) sum((sinusoid_model(b, x_data) - 
    measurements).^2);
    options = optimset('MaxFunEvals', 5000, 'MaxIter', 1000, 'TolFun', 
    1e-8);
    fitted_params = fminsearch(cost_function, initial_params, options);
    fitted_sinusoid = sinusoid_model(fitted_params, x_data);
    % Remove sinusoid
    detrended_measurements = measurements - (fitted_sinusoid - 
    fitted_params(4));
    %% STEP 3: Create new histogram from cleaned data
    [clean_hist, edges] = histcounts(detrended_measurements, 
    0.5:1:length(raw_hist)+0.5);
    %% STEP 4: Find upper and lower peaks
    % Smooth to reduce noise
    clean_hist_smooth = smoothdata(clean_hist, 'gaussian', 5);
    % Find ALL significant peaks first
    [all_pks, all_locs, all_widths, all_proms] = 
    findpeaks(clean_hist_smooth, ...
        'MinPeakProminence', max(clean_hist_smooth)*0.05, ...
        'MinPeakDistance', 50);
    if length(all_locs) == 0
        fprintf('Warning: No peaks found in row %d\n', k);
        continue;
    end
    % Find HIGHEST peak (upper) and most prominent LOWER peak
    [~, highest_idx] = max(all_pks);
    upper_peak_loc = all_locs(highest_idx);
    upper_peak_height = all_pks(highest_idx);
    upper_peak_width = all_widths(highest_idx);
    upper_peak_prom = all_proms(highest_idx);
    % Find peaks significantly BELOW upper peak (300+ bins away)
    lower_candidates = find(all_locs < (upper_peak_loc - 300));
    if isempty(lower_candidates)
        % Only upper peak detected
        locs = upper_peak_loc;
        pks = upper_peak_height;
        widths = upper_peak_width;
        proms = upper_peak_prom;
        fprintf('Note: Row %d has only upper peak detected\n', k);
    else
        % Find most prominent lower peak
        lower_proms = all_proms(lower_candidates);
        [~, best_lower_idx] = max(lower_proms);
        lower_peak_idx = lower_candidates(best_lower_idx);
        lower_peak_loc = all_locs(lower_peak_idx);
        lower_peak_height = all_pks(lower_peak_idx);
        lower_peak_width = all_widths(lower_peak_idx);
        lower_peak_prom = all_proms(lower_peak_idx);
        % Store in order: lower first, upper second
        locs = [lower_peak_loc; upper_peak_loc];
        pks = [lower_peak_height; upper_peak_height];
        widths = [lower_peak_width; upper_peak_width];
        proms = [lower_peak_prom; upper_peak_prom];
    end
    %% STEP 5: Calculate peak centers and extents
    centers = zeros(size(locs));
    lower_bounds = zeros(size(locs));
    upper_bounds = zeros(size(locs));
    half_prom_thresholds = zeros(size(locs));
    for p = 1:length(locs)
        % Define window around peak
        win_size = max(50, round(2*widths(p)));
        left = max(1, locs(p) - win_size);
        right = min(length(clean_hist_smooth), locs(p) + win_size);
        window = left:right;
        % Calculate weighted center
        data = clean_hist_smooth(window);
        weights = data / sum(data);
        centers(p) = sum(window .* weights);
        % Find extent at half-prominence threshold
        half_prom_thresholds(p) = pks(p) - proms(p)/2;
        above_threshold = clean_hist_smooth(window) > 
        half_prom_thresholds(p);
        extent_indices = window(above_threshold);
        if ~isempty(extent_indices)
            lower_bounds(p) = min(extent_indices);
            upper_bounds(p) = max(extent_indices);
        end
    end
    %% STEP 6: Enhanced Visualization
    colors = {'r', 'g'};
    labels = {'Lower', 'Upper'};
    % Plot 1: Original histogram
    subplot(length(test_rows), 3, 3*idx-2)
    plot(bin_centers, raw_hist, 'k-', 'LineWidth', 1.2);
    hold on;
    % Mark detected peak regions
    for p = 1:length(locs)
        plot(locs(p), raw_hist(locs(p)), [colors{min(p,2)} 'o'], ...
            'MarkerSize', 8, 'LineWidth', 2);
    end
    title(sprintf('Row %d: Original', k), 'FontSize', 10, 'FontWeight', 
    'bold');
    xlabel('Laser Value', 'FontSize', 9);
    ylabel('Count', 'FontSize', 9);
    xlim([2000 4000]);
    grid on;
    box on;
    hold off;
    % Plot 2: Detrended histogram with annotations
    subplot(length(test_rows), 3, 3*idx-1)
    plot(bin_centers, clean_hist_smooth, 'b-', 'LineWidth', 1.5);
    hold on;
    for p = 1:length(locs)
        color = colors{min(p, 2)};
        % Peak location
        plot(locs(p), pks(p), [color 'o'], 'MarkerSize', 8, 'LineWidth',
        2);
        % Center
        plot(centers(p), interp1(bin_centers, clean_hist_smooth, 
        centers(p)), ...
            [color 'x'], 'MarkerSize', 10, 'LineWidth', 3);
        % Extent at half-prominence
        plot([lower_bounds(p), upper_bounds(p)], ...
            [half_prom_thresholds(p), half_prom_thresholds(p)], ...
            [color '-'], 'LineWidth', 2.5);
        % Label
        text(centers(p), pks(p)*1.15, labels{min(p, 2)}, ...
            'HorizontalAlignment', 'center', 'FontWeight', 'bold', ...
            'FontSize', 9, 'Color', color);
    end
    title(sprintf('Row %d: Detrended', k), 'FontSize', 10, 'FontWeight',
    'bold');
    xlabel('Laser Value', 'FontSize', 9);
    ylabel('Count', 'FontSize', 9);
    xlim([2000 4000]);
    grid on;
    box on;
    hold off;
    % Plot 3: Overlay comparison (NO LEGEND - clear view)
    subplot(length(test_rows), 3, 3*idx)
    % Plot original
    plot(bin_centers, raw_hist, 'k-', 'LineWidth', 1.2);
    hold on;
    % Scale detrended to match for comparison
    scale_factor = max(raw_hist) / max(clean_hist_smooth);
    plot(bin_centers, clean_hist_smooth * scale_factor, ...
        'b-', 'LineWidth', 1.8, 'Color', [0 0.4470 0.7410]);
    % Mark peaks on both
    for p = 1:length(locs)
        color = colors{min(p, 2)};
        % On original
        plot(locs(p), raw_hist(locs(p)), [color 'o'], ...
            'MarkerSize', 8, 'LineWidth', 2, 'MarkerFaceColor', color);
        % On detrended (scaled)
        plot(locs(p), pks(p) * scale_factor, [color 's'], ...
            'MarkerSize', 7, 'LineWidth', 2);
    end
    title(sprintf('Row %d: Comparison', k), 'FontSize', 10, 
   'FontWeight', 'bold');
    xlabel('Laser Value', 'FontSize', 9);
    ylabel('Count', 'FontSize', 9);
    xlim([2000 4000]);
    grid on;
    box on;
    % Add text annotations instead of legend
    text(2100, max(raw_hist)*0.95, 'Black = Original', ...
        'FontSize', 8, 'BackgroundColor', 'white', 'EdgeColor', 
    'black');
    text(2100, max(raw_hist)*0.85, 'Blue = Detrended', ...
        'FontSize', 8, 'BackgroundColor', 'white', 'EdgeColor', 'black',
         ...
        'Color', [0 0.4470 0.7410]);
    hold off;
    %% STEP 7: Report results
    fprintf('\n=== Row %d Results ===\n', k);
    fprintf('Sinusoid: Amp=%.1f, Period=%.1f bins, Phase=%.2f rad\n', 
     ...
    fitted_params(1), fitted_params(2), fitted_params(3));
    fprintf('Detected %d peak(s):\n', length(locs));
    for p = 1:length(locs)
        peak_label = labels{min(p, 2)};
        fprintf('  %s Peak: Loc=%.0f, Height=%.1f, Center=%.1f, ', ...
            peak_label, locs(p), pks(p), centers(p));
        fprintf('Extent=[%.0f to %.0f] (width=%.0f bins)\n', ...
            lower_bounds(p), upper_bounds(p), upper_bounds(p)-
        lower_bounds(p));
        fprintf('    Half-prominence threshold: %.1f\n', 
        half_prom_thresholds(p));
    end
    % Store results - FIXED STRUCTURE
    results(idx).row = k;
    results(idx).num_peaks = length(locs);
    results(idx).peak_locations = locs;
    results(idx).peak_heights = pks;
    results(idx).peak_centers = centers;
    results(idx).lower_bounds = lower_bounds;
    results(idx).upper_bounds = upper_bounds;
    results(idx).half_prom_thresholds = half_prom_thresholds;
    results(idx).sinusoid_params = fitted_params;
  end
%% Final Summary with Statistics - FIXED INDEXING
fprintf('\n\n=== FINAL SUMMARY ===\n');
fprintf('Successfully processed %d rows\n\n', length(results));
lower_centers = [];
upper_centers = [];
separations = [];
for idx = 1:length(results)
  fprintf('Row %d: %d peak(s) detected\n', results(idx).row,      results(idx).num_peaks);
    if results(idx).num_peaks >= 2
        % Lower peak (first peak)
        lower_center = results(idx).peak_centers(1);
        lower_extent_width = results(idx).upper_bounds(1) -          results(idx).lower_bounds(1);
        % Upper peak (second peak)
        upper_center = results(idx).peak_centers(2);
        upper_extent_width = results(idx).upper_bounds(2) -         results(idx).lower_bounds(2);
        separation = upper_center - lower_center;
        fprintf('  Lower peak @ %.0f (extent: %.0f bins wide)\n', ...
            lower_center, lower_extent_width);
        fprintf('  Upper peak @ %.0f (extent: %.0f bins wide)\n', ...
            upper_center, upper_extent_width);
        fprintf('  Separation: %.0f bins\n', separation);
        lower_centers = [lower_centers; lower_center];
        upper_centers = [upper_centers; upper_center];
        separations = [separations; separation];
    elseif results(idx).num_peaks == 1
        % Only one peak
        only_center = results(idx).peak_centers(1);
        only_extent_width = results(idx).upper_bounds(1) -         results(idx).lower_bounds(1);
        fprintf('  Single peak @ %.0f (extent: %.0f bins wide)\n', ...
            only_center, only_extent_width);
    end
  end
% Statistics across all rows
if ~isempty(lower_centers)
  fprintf('\n=== STATISTICS ACROSS ROWS ===\n');
  fprintf('Lower peak centers: Mean=%.1f, Std=%.1f, Range=[%.0f to 
  %.0f]\n', ...
      mean(lower_centers), std(lower_centers), min(lower_centers), 
  max(lower_centers));
  fprintf('Upper peak centers: Mean=%.1f, Std=%.1f, Range=[%.0f to 
  %.0f]\n', ...
      mean(upper_centers), std(upper_centers), min(upper_centers), 
      max(upper_centers));
  fprintf('Peak separations: Mean=%.1f, Std=%.1f, Range=[%.0f to 
  %.0f]\n', ...
      mean(separations), std(separations), min(separations), 
      max(separations));
end
fprintf('\n✓ Processing complete!\n');

Accedi per commentare.

Prodotti


Release

R2025a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by