微信公众号：智能电磁频谱算法 QQ交流群：949444104
主要内容
MATLAB代码
% Parameters
M = 16;
N = 4; % Number of circles for CQAM
SNR_dB = 0:2:25; % Extended SNR range to reach higher values
num_symbols = 1e5; % Total number of symbols for simulation
% Define the constellation points for 16-PAM, 16-PSK, and 16-QAM
constellation_pam = pammod(0:(M-1), M, 0, 'gray'); % 16-PAM
constellation_psk = pskmod(0:(M-1), M, 0, 'gray'); % 16-PSK
constellation_qam = qammod(0:(M-1), M, 'gray'); % 16-QAM
constellation_cqam = generate_16CQAM(M, N);
% Calculate SEP vs SNR for each modulation scheme
SEP_vs_SNR = zeros(length(SNR_dB), 4);
for i = 1:length(SNR_dB)
SNR_linear = 10^(SNR_dB(i)/10);
SEP_vs_SNR(i, 1) = sep_pam(SNR_linear);
SEP_vs_SNR(i, 2) = sep_psk(SNR_linear);
SEP_vs_SNR(i, 3) = sep_qam(SNR_linear);
SEP_vs_SNR(i, 4) = sep_cqam(SNR_linear);
end
% Plot SEP vs SNR
figure;
semilogy(SNR_dB, SEP_vs_SNR(:, 1), 'r-', 'LineWidth', 2);
hold on;
semilogy(SNR_dB, SEP_vs_SNR(:, 2), 'g--', 'LineWidth', 2);
semilogy(SNR_dB, SEP_vs_SNR(:, 3), 'b-.', 'LineWidth', 2);
semilogy(SNR_dB, SEP_vs_SNR(:, 4), 'm:', 'LineWidth', 2);
xlabel('E_b/N_0 (dB)');
ylabel('SEP');
legend('16-PAM', '16-PSK', '16-QAM', '16-CQAM');
title('SEP vs SNR for 16-PAM, 16-PSK, 16-QAM, and 16-CQAM');
grid on;
% Function to generate 16-CQAM constellation
function constellation = generate_16CQAM(M, N)
n = M / N;
R1 = 1; % Smallest radius
constellation = [];
for i = 1:N
Ri = R1 * i; % Increasing radius for each circle
base_angle = pi/4 * (i-1); % Rotate each circle by 45 degrees more than the previous one
angles = (0:n-1) * (2 * pi / n) + base_angle;
constellation = [constellation, Ri * exp(1i * angles)]; % Append points
end
end
% Function to calculate theoretical SEP for 16-PAM
function sep = sep_pam(SNR)
M = 16;
sep = 2 * (M - 1) / M * qfunc(sqrt(6 * SNR / (M^2 - 1)));
end
% Function to calculate theoretical SEP for 16-PSK
function sep = sep_psk(SNR)
M = 16;
sep = 2 * qfunc(sqrt(2 * SNR) * sin(pi / M));
end
% Function to calculate theoretical SEP for 16-QAM
function sep = sep_qam(SNR)
M = 16;
Ps_M = 2 * (1 - 1/sqrt(M)) * qfunc(sqrt(3 / (M - 1) * SNR));
sep = 1 - (1 - Ps_M)^2;
end
% Function to calculate theoretical SEP for 16-CQAM
function sep = sep_cqam(SNR)
M = 16;
R = 1; % Smallest radius
N = 4;
sum_sq = sum((1:N).^2);
Es_avg = (4 / M) * R^2 * sum_sq;
SNR_eff = SNR * Es_avg;
sep = 2 * (M - 1) / M * qfunc(sqrt(6 * SNR_eff / (M^2 - 1)));
end% Define the constellation points for 16-PAM, 16-PSK, and 16-QAM
M = 16;
constellation_pam = pammod(0:(M-1), M, 0, 'gray'); % 16-PAM
constellation_psk = pskmod(0:(M-1), M, 0, 'gray'); % 16-PSK
constellation_qam = qammod(0:(M-1), M, 'gray'); % 16-QAM
% Normalize constellations so that average symbol energy Es = 1
constellation_pam = normalize_energy(constellation_pam);
constellation_psk = normalize_energy(constellation_psk);
constellation_qam = normalize_energy(constellation_qam);
% Generate 16-CQAM constellation
N = 4; % Change N to other values as needed
constellation_cqam = generate_16CQAM(M, N);
constellation_cqam = normalize_energy(constellation_cqam);
% Calculate Euclidean distances for each modulation
dmin_pam = calculate_dmin(constellation_pam);
dmin_psk = calculate_dmin(constellation_psk); % Calculate dmin for 16-PSK
dmin_qam = calculate_dmin(constellation_qam);
dmin_cqam = calculate_dmin(constellation_cqam);
% Calculate PAPR for each modulation
calculatePAPR = @(signal) max(abs(signal).^2) / mean(abs(signal).^2);
papr_pam = calculatePAPR(constellation_pam);
papr_psk = calculatePAPR(constellation_psk);
papr_qam = calculatePAPR(constellation_qam);
papr_cqam = calculatePAPR(constellation_cqam);
% Plot PAPR vs dmin
figure;
hold on;
scatter(dmin_pam, papr_pam, 'filled', 'DisplayName', '16-PAM');
scatter(dmin_psk, papr_psk, 'filled', 'DisplayName', '16-PSK');
scatter(dmin_qam, papr_qam, 'filled', 'DisplayName', '16-QAM');
scatter(dmin_cqam, papr_cqam, 'filled', 'DisplayName', '16-CQAM');
xlabel('d_{min}');
ylabel('PAPR');
title('Comparison of PAPR and Euclidean Distance / d_{min} in 16-PAM, 16-PSK, 16-QAM, and 16-CQAM');
legend('Location', 'best');
grid on;
hold off;
% Plot constellations
figure;
scatter(real(constellation_pam), imag(constellation_pam), 'filled', 'o');
title('16-PAM Constellation');
grid on;
axis equal;
for i = 1:length(constellation_pam)
text(real(constellation_pam(i)), imag(constellation_pam(i)), ['s', num2str(i)], 'VerticalAlignment','bottom', 'HorizontalAlignment','right')
end
figure;
scatter(real(constellation_psk), imag(constellation_psk), 'filled', 'o');
title('16-PSK Constellation');
grid on;
axis equal;
for i = 1:length(constellation_psk)
text(real(constellation_psk(i)), imag(constellation_psk(i)), ['s', num2str(i)], 'VerticalAlignment','bottom', 'HorizontalAlignment','right')
end
figure;
scatter(real(constellation_qam), imag(constellation_qam), 'filled', 'o');
title('16-QAM Constellation');
grid on;
axis equal;
for i = 1:length(constellation_qam)
text(real(constellation_qam(i)), imag(constellation_qam(i)), ['s', num2str(i)], 'VerticalAlignment','bottom', 'HorizontalAlignment','right')
end
figure;
scatter(real(constellation_cqam), imag(constellation_cqam), 'filled', 'o');
title('16-CQAM Constellation');
grid on;
axis equal;
for i = 1:length(constellation_cqam)
text(real(constellation_cqam(i)), imag(constellation_cqam(i)), ['s', num2str(i)], 'VerticalAlignment','bottom', 'HorizontalAlignment','right')
end
% Function to generate 16-CQAM constellation
function constellation = generate_16CQAM(M, N)
n = M / N;
R1 = 1; % Smallest radius
constellation = [];
for i = 1:N
Ri = R1 * i; % Increasing radius for each circle
base_angle = pi/4 * (i-1); % Rotate each circle by 45 degrees more than the previous one
angles = (0:n-1) * (2 * pi / n) + base_angle;
constellation = [constellation, Ri * exp(1i * angles)]; % Append points
end
end
% Function to calculate the minimum distance (dmin) in a constellation
function dmin = calculate_dmin(constellation)
num_points = length(constellation);
distances = inf(num_points*(num_points-1)/2, 1);
k = 1;
for i = 1:num_points
for j = i+1:num_points
distances(k) = abs(constellation(i) - constellation(j));
k = k + 1;
end
end
dmin = min(distances);
end
% Function to normalize the average symbol energy to 1
function normalized_constellation = normalize_energy(constellation)
avg_energy = mean(abs(constellation).^2);
normalized_constellation = constellation / sqrt(avg_energy);
end
MATLAB仿真结果
