代码
%% 基本参数M=240; % 产生码元数L=100; % 每个码元采样次数fc=50e3; % 载波频率50kHz% flocal = 50010; % 接收端的本地载波频率flocal = 50100; % 模拟接收端载波频率不同步的情况Rb =10e3; % 码元速率Ts=1/Rb; % 码元的持续时间dt=Ts/L; % 采样间隔TotalT=M*Ts; % 总时间t=0:dt:TotalT-dt; % 时间Fs=L*Rb; % 采样频率C1 = 2^(-4); % costas环滤波器系数c1C2 = C1 * 2^(-3); % costas环滤波器系数c2%% 产生信号源wave=randi([0,1],1,M); % 随机产生信号%帧头oxcc,23时24分25秒的一个数据包，最后一字节为校验和%wave=[1 1 0 0 1 1 0 0 0 0 0 1 0 1 1 1 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0 1 0 0];wave=2*wave-1; % 单极性变双极性fz=ones(1,L); % 定义复制的次数L,L为每码元的采样点数x1=wave(fz,:); % 将原来wave的第一行复制L次，称为L*M的矩阵baseband=reshape(x1,1,L*M); % 产生双极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵%% I、Q路码元% I路码元是基带码元的奇数位置码元，Q路码元是基带码元的偶数位置码元I=[]; Q=[];for i=1:M if mod(i, 2)~=0 I = [I, wave(i)]; else Q = [Q, wave(i)]; endendfz2 = ones(1,2*L);x2 = I(fz2,:); % 将原来I的第一行复制2L次，成为2L*(M/2)的矩阵I_signal = reshape(x2,1,L*M);% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*M)的矩阵x3 = Q(fz2,:); % 将原来Q的第一行复制2L次，称为2L*(M/2)的矩阵Q_signal = reshape(x3,1,L*M);% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*M)的矩阵%% 成形滤波% 通过Filter Designer生成了40阶(41个抽头系数)的升余弦平方根滤波器rcosfilter% 采样频率为Fs,截止频率为Rb/2Q_filtered = filter(rcosfilter,Q_signal); %Q路成形滤波I_filtered = filter(rcosfilter,I_signal); %I路成形滤波Q_filtered = double(Q_filtered);I_filtered = double(I_filtered);%% QPSK调制carry_cos=cos(2*pi*fc*t); % 载波1psk1=I_filtered.*carry_cos; % PSK1的调制carry_sin=sin(2*pi*fc*t); % 载波2psk2=Q_filtered.*carry_sin; % PSK1的调制qpsk=psk1+psk2; % QPSK的实现%% 信号经过高斯白噪声信道%qpsk_n = qpsk; %不加噪qpsk_n=awgn(qpsk,20); % 信号qpsk中加入白噪声，信噪比为SNR=20dB%% 解调部分err_phase = zeros(1,length(t));phase_ctrl= zeros(1,length(t));carry_cos_local = zeros(1, length(t));carry_sin_local = zeros(1, length(t));demo_I = zeros(1, length(t));demo_Q = zeros(1, length(t));filtered_I = zeros(1, length(t));filtered_Q = zeros(1, length(t));pd_I = zeros(1, length(t));pd_Q = zeros(1, length(t));inv_Q = zeros(1, length(t));inv_I = zeros(1, length(t));%% 载波同步与下变频for i = 1:length(t)
carry_cos_local(i) = cos(2*pi*flocal*t(i)-err_phase(i)); % 本地载波受相位控制字err_phase调整 carry_sin_local(i) = sin(2*pi*flocal*t(i)-err_phase(i)); % 本地载波受相位控制字err_phase调整 % 利用可调整频率的本地载波与QPSK信号相乘 demo_I(i)=qpsk_n(i)*carry_cos_local(i); % 相干解调，乘以本地相干载波 demo_Q(i)=qpsk_n(i)*carry_sin_local(i);
%低通滤波 filtered_Q = double(filter(demo_lowpass,demo_Q)); %Q路低通滤波 filtered_I = double(filter(demo_lowpass,demo_I)); %I路低通滤波 % 低通滤波后进行载波同步鉴相，模拟costas环鉴相器 inv_Q(i) = -1*filtered_Q(i); inv_I(i) = -1*filtered_I(i); % 依据I路正负选择I路待相乘的鉴相值 if filtered_I(i)>=0
pd_I(i) = filtered_Q(i); else
pd_I(i) = inv_Q(i); end% ind = find(filtered_I>=0);% pd_I(ind) = filtered_Q(ind);% ind = find(filtered_I<0);% pd_I(ind) = inv_Q(ind); %依据Q路正负选择Q路待相乘的鉴相值 if filtered_Q(i)>=0
pd_Q(i) = filtered_I(i); else
pd_Q(i) = inv_I(i); end % 鉴相器原始输出(未滤波) pd(i) = pd_I(i) - pd_Q(i); %鉴相器输出环路滤波 if i==1 err_phase(i+1) = C1*pd(i);
elseif i ~= length(t)
err_phase(i+1) = err_phase(i) + C1*pd(i)+(C2-C1)*pd(i-1);
end% ind = find(filtered_Q>=0);% pd_Q(ind) = filtered_I(ind);% ind = find(filtered_Q<0);% pd_Q(ind) = inv_I(ind);end% err_phase(1) = C1*pd(1); % 滤波器输入输出第一个数据是一样的% for i=2:length(t)% err_phase(i) = err_phase(i-1) + C1*pd(i)+(C2-C1)*pd(i-1);% end%% 抽样判决k=0; % 设置抽样限值sample_d_I=1*(filtered_I>k); % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0sample_d_Q=1*(filtered_Q>k); % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0%% I、Q路合并I_comb = [];Q_comb = [];% 取码元的中间位置上的值进行判决for j=L:2*L:(L*M)
if sample_d_I(j)>0 I_comb=[I_comb,1]; else I_comb=[I_comb,-1]; endend% 取码元的中间位置上的值进行判决for k=L:2*L:(L*M)
if sample_d_Q(k)>0 Q_comb=[Q_comb,1]; else Q_comb=[Q_comb,-1]; endendcode = [];% 将I路码元为最终输出的奇数位置码元，将Q路码元为最终输出的偶数位置码元for n=1:M if mod(n, 2)~=0 code = [code, I_comb((n+1)/2)]; else code = [code, Q_comb(n/2)]; endendx4=code(fz,:); % 将原来code的第一行复制L次，称为L*M的矩阵dout=reshape(x4,1,L*M); % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵%% 绘制原始信号figure();
subplot(311); % 窗口分割成3*1的，当前是第1个子图plot(t,baseband,'LineWidth',2); % 绘制基带码元波形，线宽为2title('基带信号波形');
xlabel('时间/s');
ylabel('幅度');
subplot(312); % 窗口分割成3*1的，当前是第2个子图plot(t,I_signal,'LineWidth',2); % 绘制基带码元波形，线宽为2title('I路信号波形');
xlabel('时间/s');
ylabel('幅度');
subplot(313); % 窗口分割成3*1的，当前是第3个子图plot(t,Q_signal,'LineWidth',2); % 绘制基带码元波形，线宽为2title('Q路信号波形'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1.1,1.1]) % 坐标范围限制%% 绘制成形滤波后信号figure();
subplot(211);
plot(t,Q_filtered,'LineWidth',2);% 绘制成形滤波后Q路信号title('成形滤波后Q路波形'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1,1]); % 设置坐标范围subplot(212);
plot(t,I_filtered,'LineWidth',2);% 绘制成形滤波后I路信号title('成形滤波后I路波形'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1,1]); % 设置坐标范围%% 绘制QPSK调制信号以及加噪后信号figure();
subplot(211);
plot(t,qpsk,'LineWidth',2); % 绘制基带码元波形，线宽为2title('QPSK信号波形'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1,1]); % 设置坐标范围subplot(212); % 窗口分割成2*1的，当前是第2个子图plot(t,qpsk_n,'LineWidth',2); % 绘制QPSK信号加入白噪声的波形axis([0,TotalT,-1,1]); % 设置坐标范围title('通过高斯白噪声信道后的信号');% 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签%% 绘制绘制IQ两路乘以本地相干载波后的信号figure();
subplot(211) % 窗口分割成2*1的，当前是第1个子图plot(t,demo_I,'LineWidth',2) % 绘制I路乘以相干载波后的信号axis([0,TotalT,-1,1]); % 设置坐标范围title("I路乘以相干载波后的信号") % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签subplot(212) % 窗口分割成2*1的，当前是第2个子图plot(t,demo_Q,'LineWidth',2) % 绘制Q路乘以相干载波后的信号axis([0,TotalT,-1,1]); % 设置坐标范围title("Q路乘以相干载波后的信号") % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签%% 绘制鉴相器输出figure();
subplot(311)
plot(t,pd,'LineWidth',2)
title("鉴相器计算结果"); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签subplot(312)
plot(t,pd_I,'LineWidth',2)
title("I路鉴相器输入");
xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签subplot(313)
plot(t,pd_Q,'LineWidth',2)
title("I路鉴相器输入");
xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签% %% 载波同步鉴相器结果展示% figure();% subplot(411)% plot(t,filtered_I,'LineWidth',2) % 绘制I路乘以相干载波后的信号% title("I路"); % 标题% xlabel('时间/s'); % x轴标签% ylabel('幅度'); % y轴标签%% subplot(413)% plot(t,filtered_Q,'LineWidth',2) % 绘制载波同步环路滤波器输出% title("Q路");% xlabel('时间/s'); % x轴标签% ylabel('幅度'); % y轴标签%% subplot(414)% plot(t,pd_Q,'LineWidth',2) % 绘制载波同步环路滤波器输出% title("鉴相器Q路");% xlabel('时间/s'); % x轴标签% ylabel('幅度'); % y轴标签%% subplot(412)% plot(t,pd_I,'LineWidth',2) % 绘制载波同步环路滤波器输出% title("鉴相器I路");% xlabel('时间/s'); % x轴标签% ylabel('幅度'); % y轴标签%% 载波同步鉴相器结果展示figure();
subplot(211)
plot(t,pd,'LineWidth',2) % 绘制I路乘以相干载波后的信号title("鉴相器计算结果"); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签subplot(212)
plot(t,err_phase,'LineWidth',2) % 绘制载波同步环路滤波器输出title("载波同步环路滤波器输出");
xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签%% 绘图比较本地载波和发送端载波figure()
nop=300; %由于数据很多，为了便于观察选取前nop点进行绘图start=1000; %开始观察的点的索引subplot(211) % 窗口分割成2*1的，当前是第1个子图% 绘制正弦载波plot(t(start+1:start+nop),carry_sin(start+1:start+nop),'LineWidth',2)
hold on% 绘制接收端正弦载波plot(t(start+1:start+nop),carry_sin_local(start+1:start+nop),'LineWidth',2)
hold onlegend('调制端正弦载波','接收端本地正弦载波');title("正弦载波") % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签subplot(212) % 窗口分割成2*1的，当前是第1个子图% 绘制余弦载波plot(t(start+1:start+nop),carry_cos(start+1:start+nop),'LineWidth',2)
hold on% 绘制本地余弦载波plot(t(start+1:start+nop),carry_cos_local(start+1:start+nop),'LineWidth',2)
hold onlegend('调制端余弦载波','接收端本地余弦载波');title("余弦载波") % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签%% 绘制加噪信号经过低通滤波器后的信号figure();
subplot(211)
plot(t,filtered_I,'LineWidth',2); % 绘制I路经过低通滤波器后的信号axis([0,TotalT,-1.1,1.1]); % 设置坐标范围title("I路经过低通滤波器后的信号");xlabel('时间/s');
ylabel('幅度');
subplot(212)
plot(t,filtered_Q,'LineWidth',2); % 绘制Q路经过低通滤波器后的信号axis([0,TotalT,-1.1,1.1]);
title("Q路经过低通滤波器后的信号");xlabel('时间/s');
ylabel('幅度');
%% 绘制抽样判决结果figure();subplot(311) % 窗口分割成3*1的，当前是第1个子图plot(t,sample_d_I,'LineWidth',2)% 画出经过抽样判决后的信号axis([0,TotalT,-0.1,1.1]); % 设置坐标范围title("I路经过抽样判决后的信号");xlabel('时间/s');
ylabel('幅度');
subplot(312) % 窗口分割成3*1的，当前是第2个子图plot(t,sample_d_Q,'LineWidth',2)% 画出经过抽样判决后的信号axis([0,TotalT,-0.1,1.1]); % 设置坐标范用title("Q路经过抽样判决后的信号")% 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签%% 绘图比较调制解调的信号figure()
subplot(411)
plot(t,I_signal,'LineWidth',2);% 绘制基带码元波形，线宽为2title('I路信号波形');
xlabel('时间/s');
ylabel('幅度');
subplot(412)
plot(t,sample_d_I,'LineWidth',2);title("I路经过抽样判决后的信号");subplot(413)
plot(t,Q_signal,'LineWidth',2);title('Q路信号波形');
xlabel('时间/s');
ylabel('幅度');
subplot(414)
plot(t,sample_d_Q,'LineWidth',2);title("Q路经过抽样判决后的信号");figure()
subplot(211)
plot(t,baseband,'LineWidth',2);% 绘制基带码元波形，线宽为2title('基带信号波形');
xlabel('时间/s');
ylabel('幅度');
subplot(212)
plot(t,dout,'LineWidth',2);% 绘制基带码元波形title('QPSK解调判决后信号'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1.1,1.1]) % 坐标范围限制subplot(313); % 窗口分割成3*1的，当前是第3个子图plot(t,dout,'LineWidth',2);% 绘制基带码元波形title('QPSK解调判决后信号'); % 标题xlabel('时间/s'); % x轴标签ylabel('幅度'); % y轴标签axis([0,TotalT,-1.1,1.1]) % 坐标范围限制%% 将仿真波形输出为txt文本作为testbench数据输入Width = 15; %数据位宽I_n=round(filtered_I*(2^(Width)-1));Q_n=round(filtered_Q*(2^(Width)-1));fid=fopen('dataI.txt','w');
for k=1:length(I_n)
B_s=dec2bin(I_n(k)+((I_n(k))<0)*2^Width,Width); for j=1:Width if B_s(j)=='1' tb=1; else tb=0; end fprintf(fid,'%d',tb); end fprintf(fid,'\r\n');endfprintf(fid,';');fclose(fid);fid=fopen('dataQ.txt','w');
for k=1:length(Q_n)
B_s=dec2bin(Q_n(k)+((Q_n(k))<0)*2^Width,Width); for j=1:Width if B_s(j)=='1' tb=1; else tb=0; end fprintf(fid,'%d',tb); end fprintf(fid,'\r\n');endfprintf(fid,';');fclose(fid);
波形
