using the fact that CWT = CWT(a, τ) =
 +∞
−∞
s(t)ψ∗
a,τ (t)dt,
and The |CWT|of PSK signals is constant while |CWT|
of ASK, FSK, MSK, and QAM signals is multi-step
function.
(ii) The |CWT|of the normalized ASK, PSK, and QAM
signals is constant while the |CWT| of normalized
FSK and MSK signals is multi-step function.
(iii) The statistical properties including the mean, the
variance, and higher order moments (HOM) of
wavelet transforms are different from modulation
scheme to another. These statistical properties also
differ depending on the order of modulation, since
the frequency, amplitude, and other signal properties
may change depending on the modulation order.
(iv) There are distinct peaks in wavelet transforms of different modulated signals and their normalized ones
when the Haar wavelet covers a symbol change. Note
that the median filtering helps in removing these
peaks which will affect |CWT|statistical properties.
Modify the code 
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);

            detectedModType = classify_with_additional_features(sigma_ap, sigma_af, sigma_dp, renyi_entropy, phase_variance, spectral_flatness, frequency_stability);
            
            % 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];

            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));
    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)
    % 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';
        
    % Differentiation between PSK and MSK using phase variance, spectral flatness, and frequency stability
    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
        if (phase_variance > 1.735 && spectral_flatness > 0.5595) || (frequency_stability > 209)
            modType = 'PSK';
        else
            modType = 'MSK';
        end
    else
        modType = 'Unknown';  % Assign a default label when there is no match
    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);  % Standard deviation as stability measure
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