1．试编写程序，仿真4PAM调制信号在高斯信道下的性能，并与理论分析结果相比。(1）画出两条性能曲线，一条是根据理论平均错误概率画出，另一条是仿真曲线；(2）程序的基本流程：信源产生信息比特、调制、将调制信号送入信道（产生高斯白噪声的程序）、接收端检测、将检测结果与信源原始信息比较计算误符号率和误比特率；(3）在给定信噪比下，第二步需多次重复，以得到一个平均错误概率；(4）信噪比范围：4PAM(0dB-14dB),4QAM(0dB-14dB)，间隔是1dB;也可在 BER =106左右终止。(5)信噪比计算 SNR =10log( Es /N0)=10log( REb /N0)。注意调制不能运用MATLAB内置函数pammod。注意标注中文注释。注意检测方法使用多进制调制信号软输出检测。注意一定用MATLAB编写。编写程序时注意矩阵维度要一致不要出错。注意索引值不能超出数组边界，不要出错。注意仿真曲线不能是直线。

以下是一个可能的程序实现，其中包括了注释以方便理解，但可能不是最优的实现方式：

clear all;
close all;

% 仿真参数
SNRdB = 0:1:14; % 信噪比范围
numBits = 1e6; % 仿真比特数
numSymbols = numBits/2; % 仿真符号数
PAMorder = 4; % PAM调制阶数
QAMorder = 4; % QAM调制阶数

% 生成随机比特序列
bits = randi([0 1],1,numBits);

%% PAM调制
% 将比特序列分组
dataInMatrix = reshape(bits,length(bits)/log2(PAMorder),log2(PAMorder));
% 将二进制数据转换为十进制数据
dataInMatrix = bi2de(dataInMatrix);
% 将十进制数据映射到符号
dataMod = pammod(dataInMatrix,PAMorder);
% 对调制后的信号进行归一化
dataMod = dataMod./sqrt(mean(abs(dataMod).^2));

%% QAM调制
% 将比特序列分组
dataInMatrixQAM = reshape(bits,length(bits)/log2(QAMorder),log2(QAMorder));
% 将二进制数据转换为十进制数据
dataInMatrixQAM = bi2de(dataInMatrixQAM);
% 将十进制数据映射到符号
dataModQAM = qammod(dataInMatrixQAM,QAMorder);
% 对调制后的信号进行归一化
dataModQAM = dataModQAM./sqrt(mean(abs(dataModQAM).^2));

%% 发送信号经过高斯白噪声信道
for k=1:length(SNRdB)
    % 计算噪声方差
    SNR = 10.^(SNRdB(k)/10);
    Es = mean(abs(dataMod).^2);
    N0 = Es/SNR;
    % 生成高斯白噪声
    noise = sqrt(N0/2)*(randn(1,numSymbols)+1j*randn(1,numSymbols));
    % 发送信号经过信道并加上噪声
    dataRx = dataMod + noise;
    
    %% PAM调制信号软输出检测
    % 计算标准差
    sigma = sqrt(N0/2);
    % 接收端对接收信号进行软判决
    dataSym = repmat(dataRx',1,PAMorder);
    decisionSoft = zeros(size(dataSym));
    for ii = 1:PAMorder
        decisionSoft(:,ii) = exp(-(dataSym-real(pammod(ii-1,PAMorder)).^2)/(2*sigma^2));
    end
    % 找出最大似然符号
    [~,decisionIndexes] = max(decisionSoft,[],2);
    % 将符号转换为比特
    receivedBits = de2bi(decisionIndexes-1,log2(PAMorder));
    % 计算比特错误率
    bitErrs(k) = sum(bits ~= receivedBits,'all');
    BER(k) = bitErrs(k)/numBits;

    %% QAM调制信号软输出检测
    % 接收端对接收信号进行软判决
    dataSymQAM = repmat(dataRx',1,QAMorder);
    decisionSoftQAM = zeros(size(dataSymQAM));
    for ii = 1:QAMorder
        decisionSoftQAM(:,ii) = exp(-abs(dataSymQAM-complex(qammod(ii-1,QAMorder))).^2/(2*sigma^2));
    end
    % 找出最大似然符号
    [~,decisionIndexesQAM] = max(decisionSoftQAM,[],2);
    % 将符号转换为比特
    receivedBitsQAM = de2bi(decisionIndexesQAM-1,log2(QAMorder));
    % 计算比特错误率
    bitErrsQAM(k) = sum(bits ~= receivedBitsQAM,'all');
    BERQAM(k) = bitErrsQAM(k)/numBits;
end

%% 理论误比特率
EbN0 = 10.^(SNRdB/10);
Pb = qfunc(sqrt(2*EbN0));
% 4PAM理论值
PbPAM = (2*(sqrt(PAMorder)-1)/sqrt(PAMorder))*Pb;
% 4QAM理论值
PbQAM = (4*(sqrt(QAMorder)-1)/sqrt(QAMorder))*Pb.*(1-Pb);

%% 画图
figure
semilogy(SNRdB,BER,'or-','LineWidth',1.5,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',8);
hold on;
semilogy(SNRdB,PbPAM,'b-','LineWidth',1.5);
semilogy(SNRdB,BERQAM,'sg-','LineWidth',1.5,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',8);
semilogy(SNRdB,PbQAM,'m-','LineWidth',1.5);
hold off;
grid on;
legend('PAM仿真','PAM理论','QAM仿真','QAM理论','Location','SouthWest');
xlabel('信噪比（dB）');
ylabel('比特错误率');
title('4PAM与4QAM在高斯信道下的性能比较')