%%%%%%%%%%%%%%%%%%%%% QPSK调制解调仿真 %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% QPSK_modula_sim2.m %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 时间：2022年4月20日 %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% 程序说明 %%%%%%%%%%%%%%%%%%%%%%%
%完成了QPSK调制解调仿真，比较不同信噪比下的误码率性能
%采用复数形式完成正交调制IQ
%调制方式：QPSK 编码方式：无 滚降因子：0.5
%解调方式：相干解调 译码方式：无 噪声：线性高斯白噪声
%对于复信号：SNR=Eb/N0*br/sr/sps，
%dB：SNR=EbN0_db+log10(ml/sps), ml=br/sr = 2(QPSK)
%%%%%%%%%%%%%%%%%%%%%%%% 仿真环境 %%%%%%%%%%%%%%%%%%%%%%%
%软件版本： MATLAB R2020b
clear all;
close all;
format long;
tic;
%%%%%%%%%%%%%%%%%%%%%%%% 程序主体 %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 参数设定 %%%%%%%%%%%%%%%%
bit_rate = 2000; %比特率
symbol_rate = 1000; %波特率
ml = bit_rate/symbol_rate;
sps = 16; %每个符号的采样点数
fc = 2000; %载波频率
fs = 16000; %采样频率
num_source = 2e6; %信号长度
%%%%%%%%%%%%%%%%% 信源 %%%%%%%%%%%%%%%%
%%%%%%%%随机信号
msg_source = randi([0, 1], 1, num_source);
%%%%%%%%%%%%%%%%% 发射机 %%%%%%%%%%%%%%%%
%%%%%%%%%%调制器
%%%%%串并转换
msg_source_i = msg_source(1:ml:(num_source-1));
msg_source_q = msg_source(2:ml:num_source);
%%%%QPSK星座图映射
bipolar_msg_source_i = zeros(1, num_source/2);
bipolar_msg_source_q = zeros(1, num_source/2);
%pi/4QPSK映射
for i = 1:num_source/2
if msg_source_i(i)==1 && msg_source_q(i)==1
bipolar_msg_source_i(i) = 1;
bipolar_msg_source_q(i) = 1;
elseif msg_source_i(i)==0 && msg_source_q(i)==1
bipolar_msg_source_i(i) = -1;
bipolar_msg_source_q(i) = 1;
elseif msg_source_i(i)==0 && msg_source_q(i)==0
bipolar_msg_source_i(i) = -1;
bipolar_msg_source_q(i) = -1;
elseif msg_source_i(i)==1 && msg_source_q(i)==0
bipolar_msg_source_i(i) = 1;
bipolar_msg_source_q(i) = -1;
end
end
%%%%波形观察
% figure;
% subplot(1,2,1);
% plot(bipolar_msg_source_i);
% title('I路时域波形');
% subplot(1,2,2);
% plot(bipolar_msg_source_q);
% title('Q路时域波形');
%%%%%插值（过采样）
up16_bipolar_msg_source_i = upsample(bipolar_msg_source_i, sps);
up16_bipolar_msg_source_q = upsample(bipolar_msg_source_q, sps);
%%%%%滚降滤波器，fir滤波器，线性相位
rollof_factor = 0.5; %滚降因子
span = 6;
rcos_fir = rcosdesign(rollof_factor, span, sps, 'sqrt');
% fvtool(rcos_fir, 'Analysis', 'impulse');
n1_fir = span*sps+1; %滤波器阶数
filter_delay1 = (n1_fir-1)/2; %滤波器的延迟点数
%%%%%滚降滤波
rcos_msg_source_i = filter(rcos_fir, 1, up16_bipolar_msg_source_i);
rcos_msg_source_q = filter(rcos_fir, 1, up16_bipolar_msg_source_q);
%%%%载波调制
%%%生成基带信号
base_signal = rcos_msg_source_i + 1i*rcos_msg_source_q;
%%%生成中频信号
time = 1:length(rcos_msg_source_i);
rcos_msg_source_carrier = (rcos_msg_source_i + 1j*rcos_msg_source_q).*exp(1j*2*pi*fc.*time/fs);
%%%%%%%%%%%%%%%%% 信道 %%%%%%%%%%%%%%%%
%%%%设置信噪比，单位dB
ebn0_db = -2:11;
snr = ebn0_db - 10*log10(sps/ml);
err_number = zeros(1, length(snr));
bit_err_rate = zeros(1, length(snr));
for i = 1:length(snr)
%高斯信号
rcos_msg_source_carrier_addnoise = awgn(rcos_msg_source_carrier, snr(i), 'measured');
%%%%%%%%%%%%%%%%% 接收机 %%%%%%%%%%%%%%%%
%%%%相干解调
rcos_msg_source_addnoise = rcos_msg_source_carrier_addnoise.*exp(-1j*2*pi*fc.*time/fs);
rcos_msg_source_receive_i = real(rcos_msg_source_addnoise);
rcos_msg_source_receive_q = imag(rcos_msg_source_addnoise);
%%%%匹配滤波
rcos_msg_source_match_i = filter(rcos_fir, 1, rcos_msg_source_receive_i);
rcos_msg_source_match_q = filter(rcos_fir, 1, rcos_msg_source_receive_q);
filter_delay2 = filter_delay1; %滤波器的延迟点数
%%%%最佳采样
delay_point = filter_delay1 + filter_delay2;
rcos_msg_source_mathch_i_ = rcos_msg_source_match_i(delay_point+1:sps:end);
rcos_msg_source_mathch_q_ = rcos_msg_source_match_q(delay_point+1:sps:end);
%%%%QPSK星座图
if(ebn0_db(i)==8)
scatterplot(base_signal(48+1:sps:end));
title(['Eb/N0=',num2str(ebn0_db(1,i)),'时发射端星座图']);
receive_signal = rcos_msg_source_match_i + 1j*rcos_msg_source_match_q;
scatterplot(receive_signal(96+1:sps:end));
title(['Eb/N0=',num2str(ebn0_db(1,i)),'时接收端星座图']);
end
%%%%判决
msg_source_match_judge_i = sign(rcos_msg_source_mathch_i_);
msg_source_match_judge_q = sign(rcos_msg_source_mathch_q_);
%%%%反映射
num_receive = num_source-delay_point/sps*2;
msg_receive_i = zeros(1, num_receive/2);
msg_receive_q = zeros(1, num_receive/2);
for k = 1:num_receive/2
if msg_source_match_judge_i(k)==1 && msg_source_match_judge_q(k)==1
msg_receive_i(k) = 1;
msg_receive_q(k) = 1;
elseif msg_source_match_judge_i(k)==-1 && msg_source_match_judge_q(k)==1
msg_receive_i(k) = 0;
msg_receive_q(k) = 1;
elseif msg_source_match_judge_i(k)==-1 && msg_source_match_judge_q(k)==-1
msg_receive_i(k) = 0;
msg_receive_q(k) = 0;
elseif msg_source_match_judge_i(k)==1 && msg_source_match_judge_q(k)==-1
msg_receive_i(k) = 1;
msg_receive_q(k) = 0;
end
end
%%%%%并串转换
%两路信号因为滤波都损失了6个点，所以最终信号损失了12个点
msg_receive = zeros(1,num_receive);
msg_receive(1:ml:(num_receive-1)) = msg_receive_i;
msg_receive(2:ml:num_receive) = msg_receive_q;
%%%%%%%%%%%%%%%%% 信宿 %%%%%%%%%%%%%%%%
%%%%误码/误比特率性能对比
[err_number(i), bit_err_rate(i)] = biterr(msg_source(1:num_receive), msg_receive);
end
toc
%%%%%%%%%%%%%%%%% 仿真结果 %%%%%%%%%%%%%%%%
berbpsk_th = berawgn(ebn0_db, 'psk', 2, 'nondiff');
berqpsk_th = berawgn(ebn0_db, 'psk', 4, 'nondiff');
figure;
semilogy(ebn0_db, berbpsk_th, '-o',ebn0_db, bit_err_rate,'-*', ebn0_db, berqpsk_th,'-+');
xlabel('比特信噪比Eb/N0 /dB');
ylabel('误比特率Pb');
title('不同信噪比下误比特率曲线');
legend('BPSK理论曲线','QPSK仿真曲线', 'QPSK理论曲线');
