【数字调制】OQPSK调制与解调【附MATLAB代码】

EW Frontier

于 2024-03-13 19:15:24 发布

阅读量906

点赞数 6

文章标签： matlab 信息与通信信号处理算法

本文链接：https://blog.csdn.net/m0_54016576/article/details/136688786

版权

文章来源：微信公众号：EW Frontier 关注可了解更多的雷达、通信、人工智能相关代码。

QPSK信号

QPSK信号与BPSK信号的调制方法相同，仅调制阶数不同，其表达式为：

$s_Q(t)=\sum_n\left\{a_ng(t-nT)\cos(2 \pi f_ct+\theta_0)-b_ng(t-nT)\sin(2 \pi f_ct+\theta_0)\right\}$

上式中 $\begin{matrix}a_{n}&b_{n}\\\end{matrix}$ 分别表示I、Q两路的基带信息序列，取值{-1,1}，调制阶数为 4阶， $\begin{matrix}a_{n}&b_{n}\\\end{matrix}$ 的每一对取值表征两位比特数据信息，其他参数解释同BPSK信号。

OQPSK信号

根据QPSK信号的调制原理知道，它的相邻码元最大相位差为180°。为减小相位突变，将两个正交分量I路与Q路的比特在时间上错开半个码元周期，使两路的比特不会同时改变，保证了相邻码元最大相位差为90°，由此便得到OQPSK 信号，称：偏置正交四相相移键控。其表达式可表示为：

$s_O(t)=\sum_n\left\{a_ng(t-nT)\cos(2\pi f_ct+\theta_0)-b_ng(t-nT-\frac{T}{2})\sin(2\pi f_ct+\theta_0)\right\}$

上式中如果将延迟的T/2时间去掉，即得到QPSK信号的调制。式中an表示基带信息数据序列，取值范围是{-1,1}，即调制阶数为2阶； T表示码元符号周期，fc为载波频率，0为载波的初始相位；g(t) 表示码元成型滤波器，即根升余弦滤波器，其滚降因子取0.5。

MATLAB代码

A=1;         
% 载波幅度
fc=2;      % 载波频率 
Ts=1;                                               % 码元宽度 
fs=1/Ts;
B1=fs;                                              %低通滤波器的宽度
N_sample=32;                                       % 基带码元抽样点数
N=500;                                              % 码元数    
dt=Ts/fc/N_sample;                                  % 抽样时间间隔 
T=N*Ts;                                             % 信号持续时间长度  
t=0:dt:T-dt;                                        % 时间向量 
Lt=length(t);                                       % 时间向量长度 
tx1=0;                                              % 时域波形图横坐标起点 
tx2=10;                                             % 时域波形图横坐标终点  
ty1=-2;                                             % 时域波形图纵坐标起点 
ty2=2;                                              % 时域波形图纵坐
fx1=-10;                                            % 功率谱图横坐标起点  
fx2=10;                                             % 功率谱图横坐标终点  
fy1=-40;                                         % 功率谱图纵坐标起点 
fy2=25;                                          % 功率谱图纵坐标终点
EsN0dB = 3:0.5:10                                %设定EbNo范围 
EsN0 = 10.^(EsN0dB/10);  
rt=zeros(1,length(EsN0dB));                        %初始化误码率向量 
M=4;
%产生二进制信源  
m=randn(1,N);                                    %  产生1到n的随机数
d=sign(m);                            % 将大于0的变为1小于0的变为-1
dd=sigexpand(d,fc*N_sample);          %将序列d的周期变为Ts
gt=ones(1,fc*N_sample);               % 产生宽度为Ts的矩形窗 
d_NRZ=conv(dd,gt);                    % 卷积产生基带信号   
figure(5);  
subplot(2,2,1); 
plot(t,dd(1:Lt));  
axis([tx1,tx2,ty1,ty2]);
figure(1); 
subplot(2,4,1); 
plot(t,d_NRZ(1:Lt));  
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度');  
title('基带信号时域波形图');
grid;  
[f,d_NRZf]=T2F(t,d_NRZ(1:Lt));% 进行傅里叶变换 
figure(1);  
subplot(2,4,5); 
plot(f,10*log10(abs(d_NRZf).^2/T));
axis([fx1,fx2,fy1,fy2]); 
xlabel('频率(Hz)'); 
ylabel('功率谱密度(dB/Hz)'); 
title('基带信号功率谱图'); 
grid; 
figure(4); 
subplot(2,2,1); 
plot(t,d_NRZ(1:Lt));
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('基带信号时域波形图'); 
grid;
figure(4); 
subplot(2,2,2); 
plot(f,10*log10(abs(d_NRZf).^2/T)); 
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)'); 
ylabel('功率谱密度(dB/Hz)');
title('基带信号功率谱图');
grid;
% 串并转换 
d1=[];
d2=[]; 
for i=1:N    
    if rem(i,2)==1
d1((i+1)/2)=d(i);   
    else
        d2(i/2)=d(i);   
    end
end
dd1=sigexpand(d1,2*fc*N_sample);   %功能同上             
gt1=ones(1,2*fc*N_sample);              
d_NRZ1=conv(dd1,gt1);                           
figure(1); 
subplot(2,4,2);
plot(t,d_NRZ1(1:Lt));
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度');  
title('上支路基带信号时域波形图');
grid; 
[f1,d_NRZ1f]=T2F(t,d_NRZ1(1:Lt));
figure(1); 
subplot(2,4,6); 
plot(f1,10*log10(abs(d_NRZ1f).^2/T));
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)');
ylabel('功率谱密度(dB/Hz)'); 
title('上支路基带信号功率谱图');
grid;

figure(3); 
subplot(2,4,1); 
plot(t,d_NRZ1(1:Lt)); 
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度'); 
title('上支路基带信号时域波形图');
grid;
dd2=sigexpand(d2,2*fc*N_sample);         
gt1=ones(1,2*fc*N_sample);                 
d_NRZ2=conv(dd2,gt1);                      
d_NRZ2DLY=[-ones(1,N_sample*fc),d_NRZ2(1:end-N_sample*fc)];  %进行延时，在前面添-1  
figure(1); 
subplot(2,4,3); 
plot(t,d_NRZ2DLY(1:Lt)); 
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('下支路基带信号时域波形图'); 
grid;  
[f2,d_NRZ2f]=T2F(t,d_NRZ2DLY(1:Lt));
figure(1);  
subplot(2,4,7); 
plot(f2,10*log10(abs(d_NRZ2f).^2/T));
axis([fx1,fx2,fy1,fy2]); 
xlabel('频率(Hz)');  
ylabel('功率谱密度(dB/Hz)'); 
title('下支路基带信号功率谱图');
grid;
figure(3); 
subplot(2,4,5); 
plot(t,d_NRZ2DLY(1:Lt));
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度'); 
title('下支路基带信号时域波形图');
grid; % 载波  
h1t=A*cos(2*pi*fc*t);
h2t=A*sin(2*pi*fc*t); 
figure(1); 
subplot(2,4,4);
plot(t,h1t); 
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度');  
title('载波信号时域波形图');
grid;  
[f3,h1tf]=T2F(t,h1t); 
figure(1);  
subplot(2,4,8); 
plot(f3,10*log10(abs(h1tf).^2/T)); 
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)'); 
ylabel('功率谱密度(dB/Hz)'); 
title('载波信号功率谱图');

grid;    % 生成OQPSK信号  
s_qpsk1=d_NRZ1(1:Lt).* h1t;            %上下支路分别调制 
s_qpsk2=d_NRZ2DLY(1:Lt).* h2t; 
figure(2); 
subplot(2,2,1); 
plot(t,s_qpsk1);  
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)'); 
ylabel('幅度');  
title('上支路频带信号时域波形图');
grid; 
[f4,s_qpsk1f]=T2F(t,s_qpsk1);
figure(2);
subplot(2,2,3);
plot(f4,10*log10(abs(s_qpsk1f).^2/T));
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)'); 
ylabel('功率谱密度(dB/Hz)'); 
title('上支路频带信号功率谱图'); 
grid; 
figure(2);
subplot(2,2,2);
plot(t,s_qpsk2);  
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)'); 
ylabel('幅度');  
title('下支路频带信号时域波形图');
grid;

[f5,s_qpsk2f]=T2F(t,s_qpsk2);
figure(2); 
subplot(2,2,4);
plot(f5,10*log10(abs(s_qpsk2f).^2/T));
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)'); 
ylabel('功率谱密度(dB/Hz)'); 
title('下支路频带信号功率谱图');
grid; 

s_oqpsk=s_qpsk1+s_qpsk2;                %两路信号相加得到调制后的信号
figure(8); 
subplot(2,2,1);
plot(t,s_oqpsk);  
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)'); 
ylabel('幅度'); 
title('已调信号时域波形图');
grid;

[f6,s_oqpskf]=T2F(t,s_oqpsk); 
figure(8); 
subplot(2,2,3); 
plot(f6,10*log10(abs(s_oqpskf).^2/T)); 
axis([fx1,fx2,fy1,fy2]); 
xlabel('频率(Hz)');  
ylabel('功率谱密度(dB/Hz)'); 
title('已调信号功率谱图');
grid;

% 信道加入高斯白噪声进行接收解调 
% 产生高斯白噪声
m=1; 
p1=-10; 
noise = wgn(m,Lt,p1); 
% 接收信号
y_qpsk = s_oqpsk + noise; 
figure(8); 
subplot(2,2,2); 
plot(t,y_qpsk); 
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('接收信号时域波形图 '); 
grid; 

[f7,y_qpskf]=T2F(t,y_qpsk);
figure(8); 
subplot(2,2,4);  
plot(f7,10*log10(abs(y_qpskf).^2/T));
axis([fx1,fx2,fy1,fy2]); 
xlabel('频率(Hz)');

ylabel('功率谱密度(dB/Hz)'); 
title('接收信号功率谱图 Pn=-10dB'); 
grid; 
%[t,ny_qpsk]=bpf(f7,y_qpskf,1,8);

% 相干解调 
% 通过乘法器1 
r_qpsk1 = y_qpsk .* h1t;   
figure(3);  
subplot(2,4,2);
plot(t,r_qpsk1); 
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('上支路通过乘法器信号时域波形图'); 
grid; 

%通过低通滤波器 
[f8,r_qpsk1f]=T2F(t,r_qpsk1);   
B1=1;
[t1,r_qpsk11]=lpf(f8,r_qpsk1f,B1); 
figure(3);  
subplot(2,4,3); 
plot(t1,r_qpsk11)  
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)'); 
ylabel('幅度');  
title('上支路通过低通滤波器信号时域波形图');
grid;

%抽样判决  
dd11=r_qpsk11(fc*N_sample:2*fc*N_sample:end);   
dd22=sign(dd11);
%%判决 
dd222=sigexpand(dd22,2*fc*N_sample); 
d_NRZ11=conv(dd222,gt1);          %
得到上支路信号  
d_NRZ11DLY=[-ones(1,N_sample*fc),d_NRZ11(1:end-N_sample*fc)];   %上支路信号延时  
figure(3);

subplot(2,4,4);
plot(t,d_NRZ11DLY(1:Lt)); 
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('上支路抽样判决后信号时域波形图'); 
grid; 
% 通过乘法器2  
r_qpsk2 = y_qpsk .* h2t;
figure(3);  
subplot(2,4,6); 
plot(t,r_qpsk2);  
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度'); 
title('下支路通过乘法器信号时域波形图');
grid; 

%通过低通滤波器
[f9,r_qpsk2f]=T2F(t,r_qpsk2);  
[t2,r_qpsk21]=lpf(f9,r_qpsk2f,B1);
figure(3); 
subplot(2,4,7); 
plot(t2,r_qpsk21);
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)'); 
ylabel('幅度'); 
title('下支路通过低通滤波器信号时域波形图');
grid; 
%抽样判决  
dd33=r_qpsk21(fc*N_sample:2*fc*N_sample:end);  
dd44=sign(dd33);     %判决
dd444=sigexpand(dd44,2*fc*N_sample); 
d_NRZ21=conv(dd444,gt1)               
figure(3);
subplot(2,4,8);
plot(t,d_NRZ21(1:Lt));
axis([tx1,tx2,ty1,ty2]);
xlabel('时间(S)');
ylabel('幅度'); 
title('下支路抽样判决后信号时域波形图'); 

grid;  
dd22DLY=d_NRZ11DLY(fc*N_sample:2*fc*N_sample:end);
% 并串转换 
ddd=[]; 
for s=1:N/2  
    ddd(2*s-1)=dd22DLY(s);  
    ddd(2*s)=dd44(s); 
end
rddd=sigexpand(ddd,fc*N_sample); 
%解调信号输出  
r_qpsk=conv(rddd,gt); 
qr_qpsk=[r_qpsk(2*N_sample*fc+1:N*N_sample*fc),ones(1,2*N_sample*fc)];
%去除延时  
figure(4);

subplot(2,2,3); 

plot(t,qr_qpsk(1:Lt)); 
axis([tx1,tx2,ty1,ty2]); 
xlabel('时间(S)');
ylabel('幅度');
title('解调信号时域波形图 Pn=-10dB'); 
grid;  
[f10,qr_qpskf]=T2F(t,r_qpsk(1:Lt));
figure(4);
subplot(2,2,4); 
plot(f10,10*log10(abs(qr_qpskf).^2/T));
axis([fx1,fx2,fy1,fy2]);
xlabel('频率(Hz)');
ylabel('功率谱密度(dB/Hz)');
title('解调信号功率谱图 Pn=-10dB');
grid;

%眼图  
eyediagram(y_qpsk,32,2,8);
x=d_NRZ1(1:fc*N*N_sample)+j*d_NRZ2DLY(1:fc*N*N_sample); 
%星座图  
scatterplot(x); 
axis([-4,4,-4,4]); 
grid; 
xlabel('实部');
ylabel('虚部'); 
title('星座图');  
d_NRZJ=d_NRZ(1:fc*N*N_sample); 
p=find(d_NRZJ<0);
d_NRZJ(p)=0; 
q=find(qr_qpsk<0); 
qr_qpsk(q)=0; 
r=find(d_NRZ>0); 
d_NRZJ(r)=1; 
s=find(qr_qpsk>0); 
qr_qpsk(s)=1;
figure 
[num,rt] = biterr(d_NRZJ,qr_qpsk);
Ps = erfc(sqrt(EsN0)*sin(pi/M));

semilogy(EsN0dB,Ps,'rd -'); 
xlabel('Es/N0(dB)');
ylabel('误码率');
grid on;  
title('OQPSK系统的误码率');