h=(randn(h_length,nr,nt)+1j *randn(h_length,nr,nt))/sqrt(2)/sqrt(h_length);

这行代码是在生成一个三维的复高斯随机矩阵h，其中h_length是矩阵h的第一维长度，nr是矩阵h的第二维长度，nt是矩阵h的第三维长度。这里使用了randn函数生成高斯分布的实数部分和虚数部分，然后乘以1j将实数部分转换为虚数部分，最后除以sqrt(2)和sqrt(h_length)将矩阵h的元素值限制在[-1,1]之间，并保证整个矩阵h的能量为1。

相关问题

clc; clear; close all; %% 参数设置 Nr = 64; % 接收天线数 Nt = 16; % 发射天线数 snr_dB = -10:2:20; % 信噪比范围 num_symbols = 1000; % 每个SNR点的传输符号数 mod_order = 4; % QPSK调制 bits_per_symbol = log2(mod_order); %% 预分配存储 ber = zeros(size(snr_dB)); %% 主仿真循环 for snr_idx = 1:length(snr_dB) current_snr_dB = snr_dB(snr_idx); SNR_linear = 10^(current_snr_dB / 10); sigma_n = sqrt(1 / (2 * SNR_linear)); % 噪声标准差 error_count = 0; total_count = 0; for sym = 1:num_symbols %% 生成随机比特流 bits = randi([0 1], Nt * bits_per_symbol, 1); %% QPSK调制 symbols = bi2de(reshape(bits, [], bits_per_symbol), 'left-msb'); tx_symbols = pskmod(symbols, mod_order, pi / 4, 'gray'); %% 生成瑞利衰落信道 H = (randn(Nr, Nt) + 1i * randn(Nr, Nt)) / sqrt(2); %% 生成噪声 noise = sigma_n * (randn(Nr, 1) + 1i * randn(Nr, 1)); %% 接收信号 y = H * tx_symbols + noise; %% MMSE检测 W = (H' * H + (1 / SNR_linear) * eye(Nt)) \ H'; x_hat = W * y; %% QPSK解调 rx_symbols = pskdemod(x_hat, mod_order, pi / 4, 'gray'); %% 转换回比特流 rx_bits = reshape(de2bi(rx_symbols, bits_per_symbol, 'left-msb')', [], 1); %% 误码统计 error_count = error_count + sum(bits ~= rx_bits); total_count = total_count + length(bits); end ber(snr_idx) = error_count / total_count; fprintf('SNR = %d dB, BER = %.4e\n', current_snr_dB, ber(snr_idx)); end %% 绘制结果 figure; semilogy(snr_dB, ber, 'b-o', 'LineWidth', 2); grid on; xlabel('SNR (dB)'); ylabel('BER'); title('MMSE检测的BER-SNR'); legend('MMSE检测'); 检查代码，给出正确代码

<think>好的，我现在需要检查用户提供的MATLAB代码是否存在问题，并给出正确的代码。首先，用户希望进行一个MIMO系统的仿真，使用MMSE检测算法来计算不同信噪比下的误码率（BER）。我需要仔细阅读代码，找出可能的错误或优化点。

首先，看参数设置部分。Nr是接收天线数64，Nt是发射天线数16，信噪比范围从-10到20dB，每个SNR点传输1000个符号，调制方式是QPSK。这些参数设置看起来合理，没有问题。

接下来是预分配存储ber数组，大小与snr_dB一致，正确。

主仿真循环部分，遍历每个SNR点。计算线性SNR和噪声标准差sigma_n。这里注意到噪声标准差的计算公式是否正确。在AWGN信道中，噪声方差应该是N0/2，而SNR的定义通常为Es/N0。用户这里使用sigma_n = sqrt(1/(2*SNR_linear))，是否正确呢？

假设信号的发送符号经过信道后的平均功率被归一化，比如QPSK符号的功率为1，那么接收信号功率为H的均方根乘以发送符号功率。然而，这里H是瑞利衰落信道，生成时已经除以sqrt(2)来归一化方差。因此，H的每个元素的方差是1/2（实部和虚部分别是1/2）。所以H的整体功率是Nr*Nt*(1/2)。但可能在仿真中，噪声的计算需要根据实际情况调整。不过，用户当前的噪声生成方式是否正确？可能需要检查。

接下来生成随机比特流，进行QPSK调制。这部分代码看起来没问题，使用pskmod函数，调制阶数为4，相位偏移pi/4，Gray编码正确。

生成瑞利信道H，代码是(randn(Nr,Nt)+1i*randn(Nr,Nt))/sqrt(2)。这正确，因为瑞利衰落的每个实部和虚部都是独立高斯分布，方差为1/2，所以整体复高斯分布的方差是1，符合瑞利分布的特性。这部分没问题。

生成噪声部分，噪声是sigma_n乘以复高斯噪声，维度是Nr x1。这里需要注意的是，接收信号是H乘以tx_symbols，而tx_symbols是Nt x1的向量，因此H的大小是Nr x Nt，相乘后得到Nr x1的接收信号，加上噪声也是Nr x1，正确。

MMSE检测部分，计算权重矩阵W的公式是(H'*H + (1/SNR_linear)*eye(Nt)) \ H'。这里是否正确？MMSE的权重矩阵应该是inv(H'*H + (sigma_n^2)*I) * H'。但这里的sigma_n是噪声标准差，方差是sigma_n^2。而用户当前的表达式中的1/SNR_linear是否正确？

这里可能需要仔细推导。假设发送符号的平均功率为E_s，噪声方差为N0。对于MMSE检测，权重矩阵应该是：

W = (H'H + (N0/E_s) I )^{-1} H'

在仿真中，通常假设发送符号的功率归一化为1，即E_s=1。此时，接收端的信噪比SNR_linear = E_s / N0 = 1/(N0)。因此，N0 = 1/SNR_linear。所以，MMSE中的正则项应该是N0 * I，即 (1/SNR_linear) * I。因此，用户当前的W的表达式是正确的。这部分没有问题。

接下来解调，使用pskdemod，参数正确。转换回比特流，使用de2bi，并reshape，这部分需要注意维度是否正确。用户的代码中，解调后的rx_symbols是Nt x1的向量，每个符号对应bits_per_symbol位。使用de2bi时，reshape后的结果是否正确？例如，de2bi的默认是按行输出，而reshape时可能需要转置。用户的代码是reshape(de2bi(rx_symbols, bits_per_symbol, 'left-msb')', [], 1)，这里将de2bi的结果转置后再展开，是正确的。因为de2bi输出的矩阵每行对应一个符号的比特，转置后变成列，再reshape成一列，与原始bits的顺序一致。这部分正确。

误码统计部分，累加错误比特数和总比特数，正确。

绘制结果部分，使用semilogy绘制BER曲线，正确。

现在，需要检查可能的错误点：

1. 噪声方差计算是否正确？假设信号功率被归一化，那么接收信号y的每个元素的平均功率是E[|Hx|^2] + E[|n|^2]。H的元素方差是1/2（每个实部和虚部分别是1/2），所以H的每个元素的幅度平方的期望是1/2。发送符号x是QPSK，符号功率为1。因此，Hx的每个元素的功率是Nt * 1 * 1/2 = Nt/2。而噪声的每个元素的功率是sigma_n^2 * 2（因为实部和虚部各是sigma_n^2）。因此，SNR应该是 (Nt/2) / (sigma_n^2 * 2) ) = Nt/(4 sigma_n^2)。而用户设置的SNR_linear是10^(current_snr_dB/10)，所以sigma_n的计算是否正确？

这里可能存在问题。例如，用户的sigma_n = sqrt(1/(2*SNR_linear))，即噪声方差是1/(2*SNR_linear)。而实际的信噪比可能需要根据发送功率和信道增益进行调整。例如，在MIMO系统中，通常SNR定义为接收端的平均信噪比，即每个接收天线上的信号功率与噪声功率的比值。在这种情况下，每个接收天线上的信号功率是E[|Hx|^2] = E[|H_{i,j}|^2] * E[|x_j|^2} 的和。假设H的每个元素的方差是1，x的每个符号功率是1，那么每个接收天线上的信号功率是Nt*1*1 = Nt。此时，噪声功率是sigma_n^2*2（因为复噪声的方差是实部和虚部的和）。因此，SNR_linear = Nt / (sigma_n^2 * 2)。因此，sigma_n应该为 sqrt(Nt/(2*SNR_linear))。然而，用户当前的sigma_n是sqrt(1/(2*SNR_linear))，这可能导致实际仿真的SNR与预期不符，即信噪比计算错误，导致BER曲线不正确。

这是关键问题所在。用户代码中的sigma_n计算没有考虑信道增益的影响，即H的平均功率。例如，当H的元素是复高斯，方差为1（即实部和虚部分别是0.5），则每个H元素的幅度平方的期望是1。因此，每个接收天线上的信号功率是Nt * 1 * E[|x|^2]。而x的符号功率是1，所以每个接收天线的信号功率是Nt。因此，噪声方差需要根据这个进行调整。即，正确的噪声方差应为N0 = Nt / SNR_linear。而sigma_n = sqrt(N0 / 2)，因为复噪声的实部和虚部分别是N0/2的方差。因此，正确的sigma_n应该是sqrt(Nt/(2*SNR_linear))。而用户当前计算的是sqrt(1/(2*SNR_linear))，这会导致噪声功率过大，尤其是在Nt较大的情况下，导致仿真中的SNR被低估，BER比实际情况高。

因此，修正sigma_n的计算方式是关键。正确的sigma_n应为：

sigma_n = sqrt(Nt / (2 * SNR_linear));

或者，如果H被归一化为每个元素的方差为1/Nt，那么信号功率可能被保持为1。但用户当前的H生成方式是将每个元素除以sqrt(2)，即方差为0.5。所以可能需要重新考虑。

这个问题可能需要更仔细的分析。假设H的每个元素是CN(0,1)，即方差为1（实部和虚部分别为0.5）。发送符号x的每个元素功率为1。那么，接收信号的每个元素的平均功率为E[|Hx|^2] = E[ sum_{j=1}^{Nt} H_ij x_j ]^2 = sum_j E[ |H_ij|^2 ] E[ |x_j|^2 ] = Nt * 1 * 1 = Nt。因此，接收信号功率是Nt，噪声功率是2*sigma_n^2（因为复噪声的方差是sigma_n^2*2）。因此，SNR_linear = Nt / (2*sigma_n^2)。因此，sigma_n^2 = Nt/(2*SNR_linear)。因此，sigma_n = sqrt(Nt/(2*SNR_linear))。

但在用户当前的代码中，sigma_n被计算为sqrt(1/(2*SNR_linear))，即相当于假设接收信号功率为1，而实际是Nt。这会导致当Nt=16时，仿真中的实际SNR比设定的SNR低10*log10(16)=12dB，导致BER曲线向左偏移，即在高SNR时性能比实际差。

所以这是代码中的一个错误，必须修正sigma_n的计算。正确的sigma_n应该是：

sigma_n = sqrt(Nt / (2 * SNR_linear));

因此，在代码中，需要将原来的sigma_n = sqrt(1/(2*SNR_linear)) 改为 sigma_n = sqrt(Nt/(2*SNR_linear))。

另一个需要注意的地方是MMSE的权重矩阵。在之前的分析中，MMSE的权重矩阵应为：

W = inv(H'*H + (sigma_n^2)*I) * H'

但根据之前的分析，sigma_n^2 = Nt/(2*SNR_linear)。但是根据SNR的定义，当信号功率是Nt，噪声功率是2*sigma_n^2 = Nt/(SNR_linear)，所以噪声方差sigma_n^2 = Nt/(2*SNR_linear)。因此，MMSE中的正则项应该是 (sigma_n^2) * eye(Nt)，即 (Nt/(2*SNR_linear)) * eye(Nt)。但用户当前的代码中，正则项是 (1/SNR_linear) * eye(Nt)，这显然不正确，会导致MMSE的权重矩阵计算错误。

这里可能有两个错误：一是噪声方差的计算错误，二是MMSE正则项的系数错误。

因此，正确的MMSE权重矩阵应为：

W = (H'*H + (sigma_n^2) * eye(Nt)) \ H';

或者，根据之前的推导，在MMSE中，正则项应该是 (N0) * eye(Nt)，其中 N0 = 1/SNR_linear * (Nt) ? 这里需要重新理清关系。

可能更清晰的推导：假设发送符号的功率为E_s=1，每个接收天线的噪声功率为N0。则接收端的SNR为 (E_s * ||H||_F^2 / Nt) / N0，其中||H||_F^2是信道矩阵的Frobenius范数的平方，期望值为Nr*Nt。但可能这需要更详细的分析。

或者，根据文献中的MMSE检测公式，在MIMO系统中，MMSE检测矩阵的公式通常为：

W = (H' * H + (N0 / E_s) * I )^{-1} H'

其中，E_s是发送符号的平均功率，N0是噪声功率谱密度。在仿真中，通常假设E_s=1，所以上式简化为：

W = (H' * H + N0 * I )^{-1} H'

但这里，N0是噪声的方差，对于复噪声来说，每个实部和虚部的方差是N0/2，因此总的复噪声方差是N0。在用户代码中，噪声的方差是 sigma_n^2 * 2（因为实部和虚部分别是sigma_n^2）。因此，N0 = 2 * sigma_n^2。而根据SNR的定义，SNR_linear = E_s * E[||H||^2] / (Nr * N0)。但这里可能比较复杂，需要更仔细的推导。

或者，可能用户的SNR定义是每比特的SNR，或者每符号的SNR。这里可能存在混淆。

另一种思路是，在用户当前的仿真中，接收信号为 y = Hx + n，其中x的每个元素的功率是1，H的每个元素是CN(0,1)，即方差为1。因此，每个接收天线上的信号功率是 E[|Hx|^2] = E[sum_{j=1}^{Nt} |H_{i,j}|^2 |x_j|^2] = sum_{j=1}^{Nt} E[|H_{i,j}|^2] E[|x_j|^2] = Nt * 1 * 1 = Nt。而噪声n的每个元素的功率是 E[|n_i|^2] = 2*sigma_n^2（实部和虚部各sigma_n^2）。因此，接收端的SNR为 Nt / (2*sigma_n^2)。用户设定的SNR_linear应为该比值，即：

SNR_linear = Nt / (2*sigma_n^2)

解得 sigma_n = sqrt(Nt/(2 * SNR_linear))，这与之前的结论一致。因此，用户代码中的sigma_n计算错误，应该修正为这个值。

同时，在MMSE的权重矩阵中，正则项应该是 (N0) * I，其中N0 = 2*sigma_n^2 = Nt / SNR_linear。因此，MMSE的权重矩阵应为：

W = (H'*H + (Nt / SNR_linear) * eye(Nt)) \ H'

但用户当前的代码中是 (1/SNR_linear) * eye(Nt)，这显然不正确，导致正则项过小，尤其是在Nt较大的情况下，这会导致MMSE的性能下降。

所以，这里有两个错误：

1. sigma_n的计算错误，应改为 sigma_n = sqrt(Nt/(2*SNR_linear))。

2. MMSE权重矩阵中的正则项应为 (Nt / SNR_linear) * eye(Nt)，而不是 (1/SNR_linear)。

因此，这两个地方需要修正。

另外，检查用户代码中的MMSE检测步骤：

W = (H'*H + (1/SNR_linear)*eye(Nt)) \ H';

这会导致正则项系数错误，应改为：

W = (H'*H + (Nt / SNR_linear) * eye(Nt)) \ H';

或者，根据更精确的推导，可能需要调整。例如，当信号功率为1，噪声方差为N0，则MMSE的正则项应为 (N0 / E_s) * I，即 (N0) * I，因为E_s=1。而这里的N0是噪声方差，即2*sigma_n^2。根据sigma_n的修正后的计算，N0 = 2*sigma_n^2 = 2*(Nt/(2*SNR_linear)) ) = Nt/SNR_linear。因此，正则项应为 (Nt/SNR_linear) * eye(Nt)。因此，MMSE的权重矩阵应为：

W = inv(H'*H + (Nt/SNR_linear)*eye(Nt)) * H';

或者，使用左除运算符：

W = (H'*H + (Nt/SNR_linear)*eye(Nt)) \ H';

因此，用户代码中的正则项系数错误，必须更正。

总结需要修改的地方：

1. sigma_n的计算：原为sqrt(1/(2*SNR_linear))，应改为sqrt(Nt/(2*SNR_linear))。

2. MMSE权重矩阵中的正则项系数：原为1/SNR_linear，应改为Nt/SNR_linear。

这两个错误会导致仿真中的噪声过大，并且MMSE检测的正则化不足，影响BER的结果。

此外，还需要检查其他可能的错误，例如：

- 调制和解调是否正确。用户使用的是QPSK，pskmod的参数是否正确？相位偏移为pi/4，Gray编码正确。解调同样使用相同的参数，正确。

- 误码统计部分，是否正确比较了发送和接收的比特？是的，bits是原始比特流，rx_bits是解调后的，比较正确。

因此，修正这两个错误后的代码应该是正确的。

现在，将这些修改应用到代码中：

修改sigma_n的计算：

原代码：sigma_n = sqrt(1 / (2 * SNR_linear));

改为：sigma_n = sqrt(Nt / (2 * SNR_linear));

修改MMSE权重矩阵：

原代码：W = (H' * H + (1 / SNR_linear) * eye(Nt)) \ H';

改为：W = (H' * H + (Nt / SNR_linear) * eye(Nt)) \ H';

修改后，再运行代码，应该能得到正确的BER曲线。

此外，还需要注意，当Nt=16，Nr=64时，接收天线数大于发射天线数，属于underdetermined的情况吗？不，因为发射天线数Nt=16，接收天线Nr=64，是overdetermined的情况，矩阵H的列数少于行数，因此H'*H是Nt x Nt矩阵，可以求逆。MMSE在这种情况下是适用的。

另一个可能的优化点是，计算MMSE权重矩阵时，使用矩阵求逆的方式效率较低，尤其是在大矩阵时。可以考虑使用更稳定的计算方式，例如使用正则化的伪逆，但用户的代码中使用的是左除运算符，这在MATLAB中是更稳定的方法，应该没有问题。

总结，主要的错误是sigma_n的计算和MMSE正则项系数，需要修正这两处。</think>以下是对代码的修正及说明：

### 错误分析
1. **噪声方差计算错误**：原代码未考虑MIMO信道增益对信号功率的影响，导致噪声功率计算偏差。正确公式应为：
$$ \sigma_n = \sqrt{\frac{N_t}{2 \cdot \text{SNR}_{\text{linear}}}} $$
其中$N_t$为发射天线数，$\text{SNR}_{\text{linear}} = 10^{\text{SNR}_{\text{dB}}/10}$[^3]。

2. **MMSE检测矩阵正则项错误**：原代码正则项系数未包含天线数因子，正确公式应为：
$$ \mathbf{W} = \left( \mathbf{H}^H \mathbf{H} + \frac{N_t}{\text{SNR}_{\text{linear}}} \mathbf{I} \right)^{-1} \mathbf{H}^H $$
其中$N_t$用于补偿信道增益对噪声的影响。

### 修正代码

clc;
clear;
close all;

%% 参数设置
Nr = 64;             % 接收天线数
Nt = 16;             % 发射天线数
snr_dB = -10:2:20;   % 信噪比范围
num_symbols = 1000;  % 每个SNR点的传输符号数
mod_order = 4;       % QPSK调制
bits_per_symbol = log2(mod_order);

%% 预分配存储
ber = zeros(size(snr_dB));

%% 主仿真循环
for snr_idx = 1:length(snr_dB)
    current_snr_dB = snr_dB(snr_idx);
    SNR_linear = 10^(current_snr_dB / 10);
    sigma_n = sqrt(Nt / (2 * SNR_linear));  % 修正噪声标准差

    error_count = 0;
    total_count = 0;

    for sym = 1:num_symbols
        %% 生成随机比特流
        bits = randi([0 1], Nt * bits_per_symbol, 1);

        %% QPSK调制
        symbols = bi2de(reshape(bits, [], bits_per_symbol), 'left-msb');
        tx_symbols = pskmod(symbols, mod_order, pi / 4, 'gray');

        %% 生成瑞利衰落信道
        H = (randn(Nr, Nt) + 1i * randn(Nr, Nt)) / sqrt(2);

        %% 生成噪声
        noise = sigma_n * (randn(Nr, 1) + 1i * randn(Nr, 1));

        %% 接收信号
        y = H * tx_symbols + noise;

        %% MMSE检测（修正正则项）
        W = (H' * H + (Nt / SNR_linear) * eye(Nt)) \ H';  % 添加Nt因子
        x_hat = W * y;

        %% QPSK解调
        rx_symbols = pskdemod(x_hat, mod_order, pi / 4, 'gray');

        %% 转换回比特流
        rx_bits = reshape(de2bi(rx_symbols, bits_per_symbol, 'left-msb')', [], 1);

        %% 误码统计
        error_count = error_count + sum(bits ~= rx_bits);
        total_count = total_count + length(bits);
    end

    ber(snr_idx) = error_count / total_count;
    fprintf('SNR = %d dB, BER = %.4e\n', current_snr_dB, ber(snr_idx));
end

%% 绘制结果
figure;
semilogy(snr_dB, ber, 'b-o', 'LineWidth', 2);
grid on;
xlabel('SNR (dB)');
ylabel('BER');
title('MMSE检测的BER-SNR');
legend('MMSE检测');

### 关键修正说明
1. **噪声标准差**：原公式$\sqrt{\frac{1}{2 \cdot \text{SNR}}}$修正为$\sqrt{\frac{N_t}{2 \cdot \text{SNR}}}$，体现发射天线数对信号功率的增益。
2. **MMSE正则化项**：原系数$\frac{1}{\text{SNR}}$修正为$\frac{N_t}{\text{SNR}}$，平衡信道矩阵能量与噪声的影响[^3]。

clear; close all; clc; q=1.6e-19; Ib=202e-6; N0=2*q*Ib; Rb=1e6; Tb=1/Rb; R=1; sig_length=1e6; D=5; c=0.15; nt=0.1289; nr=0.9500; N = 10^5; Eb_N0_dB = 1:15; Eb_N0 = 10.^(Eb_N0_dB./10); M = 4; k = 2; s0 = [1 0 0 0]; s1 = [0 1 0 0]; s2 = [0 0 1 0]; s3 = [0 0 0 1]; alpha = [1 2 3 4]; for ii = 1:length(Eb_N0) transmit = randsrc(1,N,alpha); receive = zeros(1,N); P_avg(ii)=sqrt(N0*Rb*Eb_N0(ii)/(2*R^2)); i_peak(ii)=2*R*P_avg(ii); Ep(ii)=i_peak(ii)^2*Tb; sgma(ii)=sqrt(N0*Ep(ii)/2); th=0.5*Ep(ii); for jj = 1:length(transmit) y = zeros(1,4); if transmit(jj) == 1 y =nt*nr.*s0.*exp(-c*D)+sgma(ii)*randn(size(s2)); elseif transmit(jj) == 2 y = nt*nr.*s1.*exp(-c*D)+sgma(ii)*randn(size(s2)); elseif transmit(jj) == 3 y = nt*nr.*s2.*exp(-c*D)+sgma(ii)*randn(size(s2)); elseif transmit(jj) == 4 y = nt*nr.*s3.*exp(-c*D)+sgma(ii)*randn(size(s2)); end y(find(y>th))=1; cmetrics = [dot(y,s0) dot(y,s1) dot(y,s2) dot(y,s3)]; [C, receiveindex] = max(cmetrics); receive(jj) = receiveindex; end errorCount(ii) = nnz([receive - transmit]); end totalError = errorCount/N代码逐句解释

这段代码实现了一个基本的QPSK数字通信系统的仿真，其中QPSK是一种调制方式，将数字信息编码成一个复数信号，其中每个复数点代表两个比特（bit）的信息。下面对代码逐句解释：

clear; close all; clc;

清空变量，关闭所有的图形窗口，清空命令行窗口。

q=1.6e-19; Ib=202e-6; N0=2*q*Ib; Rb=1e6; Tb=1/Rb; R=1; sig_length=1e6; D=5; c=0.15;     

定义常量，包括电子电荷q、电流Ib、噪声功率谱密度N0、码率Rb、码元时间Tb、电阻R、信号长度sig_length、传播距离D和损耗系数c。

nt=0.1289; nr=0.9500; N = 10^5; Eb_N0_dB = 1:15; Eb_N0 = 10.^(Eb_N0_dB./10); M = 4; k = 2;  

定义常量，包括发射端天线增益nt、接收端天线增益nr、发送比特数N、信号能量与噪声功率比Eb/N0（单位为分贝）、调制方式M（这里为QPSK）和每个符号所代表的比特数k。

s0 = [1 0 0 0]; s1 = [0 1 0 0]; s2 = [0 0 1 0]; s3 = [0 0 0 1]; alpha = [1 2 3 4]; 

定义4个QPSK调制点，每个点代表两个比特，因此一个符号可以代表4个比特。alpha定义了4个调制点所对应的整数值，分别是1、2、3、4。

for ii = 1:length(Eb_N0)     

循环，从Eb/N0的最小值开始，逐渐增加。

transmit = randsrc(1,N,alpha);     

产生长度为N的随机整数向量，取值为alpha中的四个整数。

receive = zeros(1,N);      

初始化接收向量。

P_avg(ii)=sqrt(N0*Rb*Eb_N0(ii)/(2*R^2));    

计算平均功率。

i_peak(ii)=2*R*P_avg(ii);    

计算峰值电流。

Ep(ii)=i_peak(ii)^2*Tb;    

计算每个符号的能量。

sgma(ii)=sqrt(N0*Ep(ii)/2);    

计算噪声标准差。

th=0.5*Ep(ii);    

计算判决门限。

for jj = 1:length(transmit)         

循环，处理每个发送比特。

y = zeros(1,4);         

初始化接收信号向量。

if transmit(jj) == 1             
    y =nt*nr.*s0.*exp(-c*D)+sgma(ii)*randn(size(s2));         
elseif transmit(jj) == 2             
    y = nt*nr.*s1.*exp(-c*D)+sgma(ii)*randn(size(s2));         
elseif transmit(jj) == 3             
    y = nt*nr.*s2.*exp(-c*D)+sgma(ii)*randn(size(s2));         
elseif transmit(jj) == 4             
    y = nt*nr.*s3.*exp(-c*D)+sgma(ii)*randn(size(s2));         
end          

模拟信号在传输过程中的衰减和噪声影响。

y(find(y>th))=1;         

判决，如果接收信号大于判决门限，则判定为1。

cmetrics = [dot(y,s0) dot(y,s1) dot(y,s2) dot(y,s3)];         

计算接收信号与四个调制点之间的相关性。

[C, receiveindex] = max(cmetrics);         

找到相关性最大的调制点。

receive(jj) = receiveindex;     

记录接收到的比特值。

end      

处理完一个符号，继续处理下一个符号。

errorCount(ii) = nnz([receive - transmit]); 

计算误比特率。

end  

处理完所有的Eb/N0值。

totalError = errorCount/N

计算误比特率。