function baseline = detect_nmr_baseline_plot(frequencies, intensities,sensitivity_factor)
    % Input: 
    % - 'frequencies' is a vector containing the NMR frequencies.
    % - 'intensities' is a vector containing the corresponding signal intensities.
    % - 'sensitivity_factor' is a scalar containing the sensitivity for the
    % threshold calculation 

    % Smooth the intensities to reduce noise (you can adjust the span)
    smoothed_intensities = smooth(intensities, 15);

    % Find peaks in the smoothed intensities using findpeaks
    [peaks, peak_locs] = findpeaks(smoothed_intensities, 'MinPeakProminence', 0.01, 'MinPeakDistance', 1000);

    % Calculate the left and right limits for each peak or group of peaks
    left_limits = zeros(size(peaks));
    right_limits = zeros(size(peaks));

    for i = 1:length(peak_locs)
        [~, left_limit] = min(abs(frequencies - peak_locs(i) - 0.5*(frequencies(2) - frequencies(1))));
        [~, right_limit] = min(abs(frequencies - peak_locs(i) + 0.5*(frequencies(2) - frequencies(1))));
        
        left_limits(i) = left_limit;
        right_limits(i) = right_limit;
    end

    % Calculate sensitivity for each peak and use the average as the sensitivity
    sensitivities = zeros(size(peaks));
    for i = 1:length(peaks)
        sensitivities(i) = calculate_snr(intensities, left_limits(i), right_limits(i));
    end
    sensitivity = sensitivity_factor * mean(sensitivities);

    % Calculate a threshold as a fraction of the maximum intensity
    max_intensity = max(smoothed_intensities);
    threshold = sensitivity * max_intensity;

    % Identify baseline regions by finding areas below the threshold
    baseline = intensities;
    baseline(baseline > threshold) = NaN;

    % Plot the NMR data with baseline points in red
    figure
    plot(frequencies, intensities, 'b', 'LineWidth', 1);
    hold on;
    plot(frequencies, baseline, 'r', 'LineWidth', 1);
    hold off;
    
%     text(gca,peak_locs+.02,peaks,num2str((1:numel(peaks))'))
    xlabel('Frequencies');
    ylabel('Intensities');
    title('NMR Spectrum with Baseline in Red');
    legend('Original Data', 'Baseline');

    set(gca,'Xdir','reverse');
    xlabel('Chemical shift','FontSize',10, 'FontWeight','bold');
%     axis tight    

    %% Create a function to calculate the SNR for each peak or group of peaks
    function snr = calculate_snr(intensity, left_limit, right_limit)
        peak_region = intensity(left_limit:right_limit);
        baseline_region = intensity([1:left_limit-1, right_limit+1:end]);
        snr = max(peak_region) / std(baseline_region);
    end

end

