本文为学习所用，严禁转载。

本文参考链接
https://zhuanlan.zhihu.com/p/667382398 QPSK代码及高斯白噪声如何产生
https://ww2.mathworks.cn/help/signal/ref/butter.html 滤波器
https://www.python100.com/html/4LEF79KQK398.html 低通滤波器

本实验使用matlab仿真了五种数字调制方式（OOK、2FSK、2PSK、QPSK、4QAM）在加性高斯白噪声信道中的误码率，与归一化信噪比的关系。其中码元速率为100bps，码元个数为6666，OOK、BPSK载波频率为1KHz，2FSK两个载波频率分别为1KHz，500Hz。OOK、2FSK、2PSK均采用相干解调的方式。下面是各种调制解调方式下的误码率曲线。

clear all; 
close all; 
clc;
M = 6666; % 产生码元数
L = 100; % 每码元复制L次，方便观察
Ts = 0.01; % 每个码元的宽度，即码元的持续时间
Rb = 1/Ts; % 码元速率
dt = Ts/L; % 采样间隔
Fs = 1/dt; % 采样率
TotalT = M*Ts; % 总时间
t = 0:dt:TotalT-dt; % 时间
fc1 = 10*Rb; % 载波频率是码元速率的10倍，即载波的周期是码元周期的1/10
fc2 =  4*Rb; % 2fsk另一个载波的频率
SNR = -10:0.5:10; % 信噪比范围


BER_ask2 = zeros(1,length(SNR)); % 2ASK误码率
BER_psk2 = zeros(1,length(SNR)); % BPSK误码率
BER_fsk2 = zeros(1,length(SNR)); % 2FSK误码率
BER_qpsk = zeros(1,length(SNR)); % QPSK误码率
BER_qam4 = zeros(1,length(SNR)); % 4QAM误码率


    % 产生二进制随机数据
    data_ask2 = randi([0,1],1,M);%产生0和1的二进制随机数
    data_fsk2 = randi([0,1],1,M);%产生0和1的二进制随机数
    data_psk2 = data_ask2*2 -1  ;%借助2ASK的随机数产生-1和1的二进制随机数
    data_qpsk = randi([0 3],M,1);%qpsk的码元范围是0~3
    data_qam4 = randi([0 3],M,1);%qam的码元范围是0~3

    
    % 产生单极性不归零矩形脉冲波形

    
    data_sample_ask2 = repmat(data_ask2,L,1); % 每个码元复制L次
    data_sample_ask2 = reshape(data_sample_ask2,1,M*L); % 产生单极性不归零矩形脉冲波形
    
    
    data_sample_psk2 = repmat(data_psk2,L,1); % 每个码元复制L次
    data_sample_psk2 = reshape(data_sample_psk2,1,M*L); % 产生双极性不归零矩形脉冲波形
   
    
    data_sample_fsk2 = repmat(data_fsk2,L,1); % 每个码元复制L次
    data_sample_fsk2 = reshape(data_sample_fsk2,1,M*L); % 产生单极性不归零矩形脉冲波形

    
    % 产生2ASK已调信号
carrier1 = cos(2*pi*fc1*t); % 载波1的正弦波
carrier2 = cos(2*pi*fc2*t); % 载波2的正弦波

    mod_ask2 = data_sample_ask2.*carrier1; % 2ASK的调制
    mod_psk2 = data_sample_psk2.*carrier1; % 2PSK的调制
    mod_fsk2 = data_sample_fsk2.*carrier1 + (1-data_sample_fsk2).*carrier2;% 2FSK的调制
    mod_qpsk = pskmod(data_qpsk,4,pi/4);% QPSK的调制
    mod_qam4 = qammod(data_qam4,4);% 4QAM的调制

    

for i = 1:length(SNR)
    % 通过实时测量已调信号的功率，对已调信号加入高斯白噪声
     noise_ask2 = awgn (mod_ask2,SNR(i),'measured');
     noise_psk2 = awgn (mod_psk2,SNR(i),'measured');
     noise_fsk2 = awgn (mod_fsk2,SNR(i),'measured');
     noise_qpsk = awgn (mod_qpsk,SNR(i),'measured');
     noise_qam4 = awgn (mod_qam4,SNR(i),'measured');
     % 对接受信号进行带通滤波，滤除通频带外的噪声
    [b1,a1] = butter(4,[0.8*fc1 1.2*fc1]/(Fs/2),'bandpass'); % 设计带通滤波器
    [b2,a2] = butter(4,[0.8*fc2 1.2*fc2]/(Fs/2),'bandpass'); % 设计带通滤波器
    r_ask2 = filter(b1,a1,noise_ask2); % 对接收的已调信号进行带通滤波，,滤除带外噪声
    r_psk2 = filter(b1,a1,noise_psk2); % 对接收的已调信号进行带通滤波，,滤除带外噪声
    r_fsk2_fc1 =  filter(b1,a1,noise_fsk2); % 对接收的已调信号进行带通滤波，,滤除带外噪声
    r_fsk2_fc2 =  filter(b2,a2,noise_fsk2); % 对接收的已调信号进行带通滤波，,滤除带外噪声

      % 对接收信号进行相干解调
      dem_ask2 = r_ask2.*carrier1;
      dem_psk2 = r_psk2.*carrier1;
      dem_fsk2_fc1 = 2*r_fsk2_fc1.*carrier1;
      dem_fsk2_fc2 = 2*r_fsk2_fc2.*carrier2;

      % 对解调信号进行低通滤波，滤除载波倍频分量，以获得码元
     coe_lowpass = fir1(6,2*fc1/Fs); % 低通滤波，将两倍的载波频率分量滤除
     lowpass_ask2 = filter(coe_lowpass,1,dem_ask2); 
     lowpass_psk2 = filter(coe_lowpass,1,dem_psk2);
     lowpass_fsk2_fc1 = filter(coe_lowpass,1,dem_fsk2_fc1);
     lowpass_fsk2_fc2 = filter(coe_lowpass,1,dem_fsk2_fc2); 

    
    % 对解调信号进行抽样判决
    sample_ask2 = lowpass_ask2(L/2:L:end); % 码元中点时间抽样
    sample_psk2 = lowpass_psk2(L/2:L:end); % 码元中点时间抽样
    sample_fsk2_fc1 = lowpass_fsk2_fc1(L/2:L:end); % 码元中点时间抽样
    sample_fsk2_fc2 = lowpass_fsk2_fc2(L/2:L:end); % 码元中点时间抽样
    decision_ask2 = (sample_ask2>0.5);
    decision_psk2 = (sample_psk2>0);
    decision_fsk2 = (sample_fsk2_fc1>sample_fsk2_fc2); 
    decision_qpsk = pskdemod(noise_qpsk,4,pi/4);
    decision_qam4 = qamdemod(noise_qam4,4);
    % 计算误码个数和误码率
    error_ask2 = sum(xor(data_ask2,decision_ask2));
    error_psk2 = sum(xor(data_ask2,decision_psk2)); %这里因为PSK是ASK的随机序列产生的，所以和ASK比较
    error_fsk2 = sum(xor(data_fsk2,decision_fsk2)); 
    error_qpsk = sum(xor(data_qpsk,decision_qpsk)); 
    error_qam4 = sum(xor(data_qam4,decision_qam4)); 
    BER_ask2(i) = error_ask2/M; % 2ASK误码率
    BER_psk2(i) = error_psk2/M; % 2PSK误码率
    BER_fsk2(i) = error_fsk2/M; % 2FSK误码率
    BER_qpsk(i) = error_qpsk/M; % QPSK误码率
    BER_qam4(i) = error_qam4/M; % 4QAM误码率
end
% 计算理论误码率
% Pe = zeros(1,length(SNR));
% for i = 1:length(SNR)
% r = 10^(SNR(i)/10);
% Pe(i) = qfunc(sqrt(r)); % 2ASK理论误码率公式
% end
% 绘制误码率曲线

% semilogy横轴是线性，纵轴10倍一格
semilogy(SNR,BER_ask2,'b-d','LineWidth',2); % 2ASK仿真曲线
hold on;
semilogy(SNR,BER_psk2,'y-h','LineWidth',2); % BPSK仿真曲线
hold on;
semilogy(SNR,BER_fsk2,'g-o','LineWidth',2); % 2FSK仿真曲线
hold on;
semilogy(SNR,BER_qpsk,'k-*','LineWidth',2); % QPSK仿真曲线
hold on;
semilogy(SNR,BER_qam4,'r-s','LineWidth',2); % 4QAM仿真曲线
hold on;
% semilogy(SNR,Pe,'m-+'); % 理论曲线
% grid on;

xlabel('归一化信噪比Eb/N0 (dB)');
ylabel('误码率BER');
legend('OOK','BPSK','2FSK','QPSK','4QAM');
title('误码率曲线');

经过上述仿真得到误码率曲线如下。