function classify_signal_with_additional_features()
    % Parameters
    fs = 1000;  % Sampling frequency
    fc = 100;   % Carrier frequency
    dataLength = fs;  % Length of data (1 second of data)
    numTests = 1000;  % Number of signals to test per modulation type

    modTypes = {'ASK', 'FSK', 'PSK', 'AM', 'FM', 'MSK'};

    correctCount = 0;
    totalCount = 0;
    actualModTypes = {};
    predictedModTypes = {};
    
    % Store feature values for logging
    featureValues = struct();
    for modType = modTypes
        featureValues.(modType{1}) = [];
    end

    % Loop over each modulation type for tests
    for i = 1:length(modTypes)
        modType = modTypes{i};

        for j = 1:numTests
            data = randi([0 1], 1, dataLength);
            signal = generateModulatedSignal(data, fs, fc, modType);
            [sigma_ap, sigma_af, sigma_dp] = extract_features(signal, fs);
            renyi_entropy = compute_renyi_entropy_stft(signal, fs);
            phase_variance = compute_phase_variance(signal, fs);
            spectral_flatness = compute_spectral_flatness(signal);
            frequency_stability = compute_frequency_stability(signal, fs);

            % Add the nonlinear phase peak factor
            nonlinear_phase_peak_factor = compute_differential_nonlinear_phase_peak_factor(signal, fs, fc);
            
            detectedModType = classify_with_additional_features(sigma_ap, sigma_af, sigma_dp, renyi_entropy, phase_variance, spectral_flatness, frequency_stability, nonlinear_phase_peak_factor);
            
            % Append to classification results
            actualModTypes{end+1} = modType;
            predictedModTypes{end+1} = detectedModType;
            
            % Store the feature values
            featureValues.(modType){end+1} = [sigma_ap, sigma_af, sigma_dp, renyi_entropy, phase_variance, spectral_flatness, frequency_stability, nonlinear_phase_peak_factor];

            if strcmp(detectedModType, modType)
                correctCount = correctCount + 1;
            end
            totalCount = totalCount + 1;
        end
    end
    
    accuracy = (correctCount / totalCount) * 100;
    fprintf('Total Accuracy: %.2f%%\n', accuracy);
    
    % Print feature values for all modulation types
    for i = 1:length(modTypes)
        modType = modTypes{i};
        fprintf('\n%s Feature Averages (over %d tests):\n', modType, numTests);
        featuresMatrix = cell2mat(featureValues.(modType)');
        avgFeatures = mean(featuresMatrix, 1);
        fprintf('Sigma_ap (Amplitude STD): %.4f\n', avgFeatures(1));
        fprintf('Sigma_af (Frequency STD): %.4f\n', avgFeatures(2));
        fprintf('Sigma_dp (Phase STD): %.4f\n', avgFeatures(3));
        fprintf('Renyi Entropy: %.4f\n', avgFeatures(4));
        fprintf('Phase Variance: %.4f\n', avgFeatures(5));
        fprintf('Spectral Flatness: %.4f\n', avgFeatures(6));
        fprintf('Frequency Stability: %.4f\n', avgFeatures(7));
        fprintf('Nonlinear Phase Peak Factor: %.4f\n', avgFeatures(8));
    end
    
    % Plot the confusion matrix
    figure;
    confusionchart(actualModTypes, predictedModTypes, 'RowSummary', 'row-normalized', 'ColumnSummary', 'column-normalized');
    title('Confusion Matrix for Modulation Classification');
end

function modType = classify_with_additional_features(sigma_ap, sigma_af, sigma_dp, renyi_entropy, phase_variance, spectral_flatness, frequency_stability, nonlinear_phase_peak_factor)
    % Updated classification rules based on expanded features

    if sigma_ap >= 0.38 && sigma_ap <= 0.39 && sigma_af >= 41 && sigma_af <= 45 && sigma_dp >= 181 && sigma_dp <= 182 && renyi_entropy >= 1.9 && renyi_entropy <= 2.2
        modType = 'ASK';
    elseif sigma_ap >= 0.30 && sigma_ap <= 0.33 && sigma_af >= 141 && sigma_af <= 170 && sigma_dp >= 249 && sigma_dp <= 274 && renyi_entropy > 2.3
        modType = 'FSK';
    elseif sigma_ap < 0.13 && sigma_af < 2.1 && sigma_dp < 182 && renyi_entropy > 2.2 && renyi_entropy < 2.4
        modType = 'AM';
    elseif sigma_ap < 0.08 && sigma_af >= 16 && sigma_af <= 24 && sigma_dp >= 223 && sigma_dp <= 230 && renyi_entropy >= 1.7 && renyi_entropy <= 2.6
        modType = 'FM';
    elseif sigma_ap >= 0.3351 && sigma_ap <= 0.3556 && sigma_af >= 204.28 && sigma_af <= 227.49 && sigma_dp >= 319.30 && sigma_dp <= 364.10 && renyi_entropy >= 2.29 && renyi_entropy <= 2.31
        % Use nonlinear phase peak factor to differentiate
        if nonlinear_phase_peak_factor < 100
            modType = 'MSK';
        else
            modType = 'PSK';
        end
    else
        modType = 'Unknown';
    end
end

function frequency_stability = compute_frequency_stability(signal, fs)
    % Compute the instantaneous frequency
    phase = unwrap(angle(hilbert(signal)));
    instantaneous_frequency = diff(phase) * fs / (2 * pi);
    frequency_stability = std(instantaneous_frequency);
end

function [stft_coeff, f, t] = stft(x, win_len, overlap, nfft, fs)
    % Implementing a basic Short Time Fourier Transform (STFT)
    window = hamming(win_len);
    step = win_len - overlap;
    idx = 0:step:(length(x)-win_len);
    num_frames = length(idx);
    stft_coeff = zeros(nfft/2 + 1, num_frames);
    
    for i = 1:num_frames
        frame = x(idx(i) + (1:win_len)) .* window;
        X = fft(frame, nfft);
        stft_coeff(:, i) = X(1:nfft/2 + 1);
    end
    
    % Frequency and time vectors
    f = (0:nfft/2) * fs / nfft;
    t = idx / fs;
end

function signal = generateModulatedSignal(data, fs, fc, modType)
    dataLength = length(data);
    t = (0:dataLength-1) / fs;
    
    switch modType
        case 'ASK'
            signal = (1 + data) .* cos(2 * pi * fc * t);
        case 'FSK'
            f1 = fc;  
            f2 = fc + 50;  
            signal = cos(2 * pi * (f1 + (f2-f1) * data) .* t);
        case 'PSK'
            signal = cos(2 * pi * fc * t + pi * data);
        case 'AM'
            dataMod = movmean(data, 20); % Low-pass filter for amplitude modulation characteristic
            signal = (1 + dataMod) .* cos(2 * pi * fc * t);
        case 'FM'
            kf = 50; % Frequency sensitivity
            signal = cos(2 * pi * fc * t + 2*pi*kf*cumsum(data)/fs);
        case 'MSK'
            % Generate MSK: Continuous-phase frequency shift keying
            kf = pi / (2*fs); % Frequency deviation for MSK
            data_diff = [data(1), diff(data)];
            phase = pi * cumsum(data_diff);
            signal = cos(2 * pi * fc * t + phase);
        otherwise
            error('Unknown modulation type');
    end
end

function [sigma_ap, sigma_af, sigma_dp] = extract_features(x, fs)
    % Amplitude statistics
    x_envelope = abs(hilbert(x));  
    sigma_ap = std(x_envelope);    
    
    % Phase statistics
    phase = unwrap(angle(hilbert(x)));  
    sigma_dp = std(phase);   
    
    % Frequency statistics
    freqDiff = fs * diff(phase) / (2*pi); % Frequency deviation based on phase difference
    sigma_af = std(freqDiff);  
end

function re = compute_renyi_entropy_stft(signal, fs)
    % Parameters for STFT
    windowLength = 256;
    overlap = windowLength / 2;
    nfft = 512; % Number of FFT points
    
    % Perform STFT (Short-Time Fourier Transform)
    [stft_coefficients, ~, ~] = stft(signal, windowLength, overlap, nfft, fs);
    
    % Compute the magnitude spectrogram (i.e., energy)
    energy = abs(stft_coefficients) .^ 2;
    total_energy = sum(energy(:));
    
    % Avoid division by zero in subsequent calculations
    if total_energy == 0
        total_energy = eps; % Small constant to avoid division by zero
    end
    
    % Renyi entropy calculation (order 2)
    p = energy / total_energy;
    re = -log(sum(p(:) .^ 2));
end

function phase_variance = compute_phase_variance(signal, fs)
    phase = unwrap(angle(hilbert(signal)));
    phase_variance = var(diff(phase));
end

function spec_flatness = compute_spectral_flatness(signal)
    % Compute the power spectral density
    [pxx,~] = periodogram(signal);
    % Spectral flatness calculation
    geometric_mean = exp(mean(log(pxx + eps)));
    arithmetic_mean = mean(pxx + eps);
    spec_flatness = geometric_mean / arithmetic_mean;
end

function pf = compute_differential_nonlinear_phase_peak_factor(signal, fs, fc)
    phase = unwrap(angle(hilbert(signal)));
    
    % Calculate the corrected phase
    corrected_phase = phase - (2 * pi * (0:length(phase)-1) * fc / fs);

    % Normalize and limit range to [0,1]
    normalized_phase = corrected_phase - min(corrected_phase);
    normalized_phase = normalized_phase / max(normalized_phase);
    
    % Calculate the difference
    phase_diff = diff(normalized_phase);
    
    % Calculate PF as described
    threshold = 0.01;  % Select an appropriate threshold for significant peaks
    peaks = phase_diff > threshold;
    sum_peaks = sum(peaks);
    mean_phase_diff = mean(phase_diff);
    
    % To avoid division by zero
    if mean_phase_diff == 0
        mean_phase_diff = eps;
    end
    
    pf = sum_peaks / mean_phase_diff;
end
Also plot a confusion matrix for this