QPSK BER theoretical and simulation doesn't match

I'm trying to plot the QPSK BER curves for 2 cases, which are:

1- Only AWGN exists

2- Rayleigh fading exists with AWGN

but there's a problem, the theoretical and simulation curves do not match. I've seen other people asking the same question about the curves not matching, but all of them were using a different methods and I'm not allowed to use a different method than ours.

I tried many things such as changing the input type of the pskmod and pskdemod, but the below code was the best I could get. I can't figure out what I'm doing wrong.

Here's the current code:

%% Clean Up
clc;
clear;
close all;
%% Parameters Definition
Eb_No = 0:20; % [dB]
Eb_No_linear = db2pow(Eb_No);
N = 64; % number of symbols
N_iterations = 1000; 
M = 4;
Im = log2(M);
Nb = N * Im;
a_msk = 0.85;
a_gmsk = 0.68;
error_awgn = zeros(length(Eb_No), N_iterations);
error_rayleigh_awgn = zeros(length(Eb_No), N_iterations);
BER_awgn = zeros(size(Eb_No));
BER_rayleigh_awgn = zeros(size(Eb_No));
%% Calculating The Theoretical Results
BER_theoretical_awgn = qfunc(sqrt(Eb_No_linear));
BER_theoretical_fading = 1/sqrt(2) * (1 - sqrt( (Eb_No_linear) ./ (2+Eb_No_linear) ));
%% Calculating The Simulation Results
for n = 1:length(Eb_No)
    Es_No = Eb_No(n) + 10*log10(Im);
    for k = 1:N_iterations
        data = randi([0 1], Nb, 1);
        x = pskmod(data, M, pi/M, "bin");
        h = 1/2 * (randn(size(x)) + 1i*randn(size(x)));
        xf = abs(h) .* x;
        y = awgn(x, Es_No, "measured");
        yf = awgn(xf, Es_No, "measured");
        y_est = pskdemod(y, M, pi/M, 'bin');
        yf_est = pskdemod(yf, M, pi/M, 'bin');
        error_awgn(n, k) = nnz(data ~= y_est);
        error_rayleigh_awgn(n, k) = nnz(data ~= yf_est);
    end
    error_awgn_avg = sum(error_awgn(n, :)) / N_iterations;
    error_rayleigh_avg = sum(error_rayleigh_awgn(n, :)) / N_iterations;
    BER_awgn(n) = error_awgn_avg / Nb;
    BER_rayleigh_awgn(n) = error_rayleigh_avg / Nb;
end
%% Plotting The Theoretical and Simulation results
figure(1);
semilogy(Eb_No, BER_theoretical_awgn, 'r-', 'LineWidth', 1.5);
hold on;
semilogy(Eb_No, BER_awgn, 'r-o', 'LineWidth', 1);
semilogy(Eb_No, BER_theoretical_fading, 'b-', 'LineWidth', 1.5);
semilogy(Eb_No, BER_rayleigh_awgn, 'b-o', 'LineWidth', 1);
grid on;
ylim([10e-5 0.5])
hold off;
ylabel('BER');
xlabel('E_b/N_o');
title('QPSK BER curves for AWGN and Rayleigh');
legend('AWGN_{theory}', 'AWGN_{simulation}', 'Rayleigh_{theory}', 'Rayleigh_{simulation}', 'location', 'best');

Answer 1

[Shashi Kiran]

Hi @Haider,

I found a few issues in the code that caused a mismatch between the theoretical and simulated Bit Error Rate values.

Here are the issues and their corresponding fixes:

1. Directly used binary bits (data) with "pskmod", but "pskmod" requires integer symbols.

• Perform bit-to-symbol conversion using "bi2de" to correctly map input bits to QPSK symbols.

2. Used 1/2 instead of the correct factor for Rayleigh fading, resulting in incorrect channel power.

• Normalize the channel coefficient using to ensure an average power of 1.

3. The code only used "abs(h)", ignoring phase shifts, which QPSK is sensitive to. Additionally, no equalization was performed, leading to channel distortion.

• Use the full complex channel coefficient h for both magnitude and phase effects, and apply equalization by dividing the received signal by h.

4. Theoretical BER Calculation:

• For QPSK in AWGN, use .
• For Rayleigh fading, use

5. Errors were averaged per iteration using only Nb, which was incorrect.

• Accumulate total bit errors across all iterations and divide by the total number of transmitted bits (Nb * N_iterations) for accurate BER estimation.

Here is the full corrected code:

%% Clean Up
clc;
clear;
close all;

%% Parameters Definition
Eb_No = 0:20; % [dB]
Eb_No_linear = db2pow(Eb_No); % Linear scale
N = 64; % number of symbols
N_iterations = 1000; 
M = 4; % QPSK
Im = log2(M); % Bits per symbol
Nb = N * Im; % Total number of bits per iteration

% Initialize error arrays
error_awgn = zeros(length(Eb_No), N_iterations);
error_rayleigh_awgn = zeros(length(Eb_No), N_iterations);
BER_awgn = zeros(size(Eb_No));
BER_rayleigh_awgn = zeros(size(Eb_No));

%% Calculating The Theoretical Results
BER_theoretical_awgn = qfunc(sqrt(2 * Eb_No_linear)); % Corrected QPSK BER for AWGN
BER_theoretical_fading = 0.5 * (1 - sqrt(Eb_No_linear ./ (1 + Eb_No_linear))); % Corrected BER for Rayleigh fading

%% Calculating The Simulation Results
for n = 1:length(Eb_No)
    Es_No = Eb_No(n) + 10*log10(Im); % Es/No in dB
    total_bit_errors_awgn = 0;
    total_bit_errors_rayleigh = 0;
    
    for k = 1:N_iterations
        % Generate random binary data
        data = randi([0 1], Nb, 1);
        
        % Map bits to symbols (0, 1, 2, 3 for QPSK)
        data_symbols = bi2de(reshape(data, Im, []).', 'left-msb');
        
        % QPSK Modulation with Gray coding
        x = pskmod(data_symbols, M, pi/M, 'gray');
        
        % Rayleigh Fading Channel Coefficient
        h = 1/sqrt(2) * (randn(N, 1) + 1i * randn(N, 1)); % Normalized
        
        % AWGN Channel
        y = awgn(x, Es_No, "measured");
        
        % Rayleigh + AWGN Channel
        xf = h .* x; % Apply full Rayleigh fading (magnitude and phase)
        yf = awgn(xf, Es_No, "measured");
        
        % Equalize Rayleigh channel effects
        yf_equalized = yf ./ h; % Equalize by dividing with channel coefficient
        
        % Demodulation
        y_est = pskdemod(y, M, pi/M, 'gray');
        yf_est = pskdemod(yf_equalized, M, pi/M, 'gray');
        
        % Convert received symbols back to bits
        y_est_bits = de2bi(y_est, Im, 'left-msb').';
        yf_est_bits = de2bi(yf_est, Im, 'left-msb').';
        
        % Reshape bits for comparison
        y_est_bits = y_est_bits(:);
        yf_est_bits = yf_est_bits(:);
        
        % Count bit errors
        total_bit_errors_awgn = total_bit_errors_awgn + nnz(data ~= y_est_bits);
        total_bit_errors_rayleigh = total_bit_errors_rayleigh + nnz(data ~= yf_est_bits);
    end
    
    % Calculate BER for current Eb/No
    BER_awgn(n) = total_bit_errors_awgn / (Nb * N_iterations);
    BER_rayleigh_awgn(n) = total_bit_errors_rayleigh / (Nb * N_iterations);
end

%% Plotting The Theoretical and Simulation Results
figure(1);
semilogy(Eb_No, BER_theoretical_awgn, 'r-', 'LineWidth', 1.5); hold on;
semilogy(Eb_No, BER_awgn, 'r-o', 'LineWidth', 1);
semilogy(Eb_No, BER_theoretical_fading, 'b-', 'LineWidth', 1.5);
semilogy(Eb_No, BER_rayleigh_awgn, 'b-o', 'LineWidth', 1);
grid on;
ylim([10e-5 0.5]);
hold off;
ylabel('BER');
xlabel('E_b/N_o');
title('QPSK BER curves for AWGN and Rayleigh');
legend('AWGN_{theory}', 'AWGN_{simulation}', 'Rayleigh_{theory}', 'Rayleigh_{simulation}', 'location', 'best');

Refer to the following documentation for more details about the function "pskomd": https://www.mathworks.com/help/comm/ref/pskmod.html

Hope this solves your query.

[Haider]

Thanks a lot for explaining what I did wrong, but can you give me a reference where I can find the BER equations for both AWGN and Rayleigh? and what is the meaning of this h vector that represents Rayleigh fading? I'd really appreciate it if you gave me sources I can read to understand the topic even further.