本文為學習所用，嚴禁轉載。

本文參考鏈接
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~3data_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('誤碼率曲線');

經過上述仿真得到誤碼率曲線如下。