function classify_signal_with_entropy_expanded()
    % 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 = {};
    
    % 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);
            detectedModType = classify_with_entropy_expanded(sigma_ap, sigma_af, sigma_dp, renyi_entropy);
            
            % Append to classification results
            actualModTypes{end+1} = modType;
            predictedModTypes{end+1} = detectedModType;

            if strcmp(detectedModType, modType)
                correctCount = correctCount + 1;
            end
            totalCount = totalCount + 1;
        end
    end
    
    accuracy = (correctCount / totalCount) * 100;
    fprintf('Total Accuracy: %.2f%%\n', accuracy);
    
    % Plot the confusion matrix
    figure;
    confusionchart(actualModTypes, predictedModTypes, 'RowSummary', 'row-normalized', 'ColumnSummary', 'column-normalized');
    title('Confusion Matrix for Modulation Classification');
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 modType = classify_with_entropy_expanded(sigma_ap, sigma_af, sigma_dp, renyi_entropy)
    % Updated classification rules based on provided average 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.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
        modType = 'PSK';
    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.327 && sigma_ap <= 0.3568 && sigma_af >= 190.55 && sigma_af <= 223.19 && sigma_dp >= 328.80 && sigma_dp <= 376.86 && renyi_entropy >= 2.29 && renyi_entropy <= 2.31
        modType = 'MSK';
    else
        modType = 'Unknown';  % Assign a default label when there is no match
    end
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: A form of 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
The MSK and PSK signals are not being differentiated correctly. Find something new parameter to differentiate and classifiy it.