
% Write a Matlab script that creates an EMG peak, then attempt to remove
% the exponential broadening from an EMG by computing its first derivative,
% multiplying it by an adjustable weight, and adding that to the original
% EMG. Start with a positive weight equal to the FWHM of the original EMG
% peak and then adjust the weight iteratively until the resulting
% deconvoluted peak is as narrow as possible without dipping below y=0.
% Plot the original Gaussian peak, the EMG peak, the derivative, and the
% deconvoluted peak all on one plot with different colors and line styles.

% Clear previous figures and workspace
close all;
clear;
clc;

% Parameters
mu = 0;          % Peak center
sigma = 0.5;     % Gaussian standard deviation
lambda = 1.0;    % Exponential decay parameter
x = linspace(-5, 5, 1000);

% Gaussian function
gaussian = @(x, mu, sigma) (1 / (sigma * sqrt(2*pi))) * exp(-0.5 * ((x - mu) / sigma).^2);

% Exponentially Modified Gaussian (EMG) function
emg = @(x, mu, sigma, lambda) ...
    lambda/2 * exp((lambda/2) * (2*mu + lambda*sigma^2 - 2*x)) .* ...
    (1 - erf((mu + lambda*sigma^2 - x) / (sqrt(2)*sigma)));

% Compute Gaussian peak
gauss_peak = gaussian(x, mu, sigma);

% Compute EMG peak
emg_peak = emg(x, mu, sigma, lambda);

% Compute derivative of EMG
emg_derivative = gradient(emg_peak, x);

% Function to compute FWHM
function fwhm = calculate_fwhm(x, y)
    half_max = max(y) / 2;
    [~, idx1] = min(abs(y(1:floor(end/2)) - half_max));
    [~, idx2] = min(abs(y(floor(end/2):end) - half_max));
    idx2 = idx2 + floor(length(y)/2) - 1;
    fwhm = x(idx2) - x(idx1);
end

% Calculate initial FWHM of EMG peak
initial_fwhm = calculate_fwhm(x, emg_peak);

% Iterative deconvolution
best_width = Inf;
best_weight = initial_fwhm;
best_deconvoluted_peak = [];

for weight = [initial_fwhm, initial_fwhm*0.5, initial_fwhm*1.5, initial_fwhm*2]
    % Deconvolution by adding weighted derivative
    deconvoluted_peak = emg_peak + weight * emg_derivative;
    
    % Check if peak remains non-negative
    if all(deconvoluted_peak >= 0)
        % Calculate new FWHM
        current_fwhm = calculate_fwhm(x, deconvoluted_peak);
        
        % Update best solution if narrower
        if current_fwhm < best_width
            best_width = current_fwhm;
            best_weight = weight;
            best_deconvoluted_peak = deconvoluted_peak;
        end
    end
end

% Plotting
figure('Position', [100, 100, 1000, 600]);

% Plot original Gaussian
plot(x, gauss_peak, 'k-', 'LineWidth', 2, 'DisplayName', 'Original Gaussian');
hold on;

% Plot EMG peak
plot(x, emg_peak, 'r--', 'LineWidth', 2, 'DisplayName', 'EMG Peak');

% Plot EMG derivative
plot(x, emg_derivative, 'g:', 'LineWidth', 2, 'DisplayName', 'EMG Derivative');

% Plot deconvoluted peak
plot(x, best_deconvoluted_peak, 'b-.', 'LineWidth', 2, 'DisplayName', 'Deconvoluted Peak');

% Formatting
title('EMG Peak Deconvolution');
xlabel('x');
ylabel('Intensity');
legend('show');
grid on;

% Display best weight
fprintf('Best weight for deconvolution: %f\n', best_weight);
fprintf('Resulting FWHM: %f\n', best_width);