通信原理与MATLAB(十):QPSK的调制解调_qpsk调制

作者：weixin_40725706 | 2024-08-20 18:19:18

踩

qpsk调制

1.QPSK的调制原理

QPSK调制原理如下图所示，QPSK相当于两个正交的BPSK相加而成。其调制原理是将基带码元分成I、Q两路，I路是原始基带码元的奇数位置码元，Q路是原始基带码元的偶数位置码元，然后两条支路分别和对应的载波相乘实现BPSK的调制，然后将两条支路相加实现QPSK的调制。
在这里插入图片描述

2.QPSK的解调原理

QPSK的解调原理如下图所示，DPSK信号再分为I、Q两路和对应的载波相乘，然后经过低通滤波器后进行抽样判决，相当于作两路的BPSK解调。判决之后的I、Q路码元进行合并，I路为最终码元序列的奇数位置码元，Q路为最终码元序列的偶数位置码元，恢复出原始的码元序列。
在这里插入图片描述

3.QPSK代码

clear all;                  % 清除所有变量
close all;                  % 关闭所有窗口
clc;                        % 清屏
%% 基本参数
M=20;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;                   % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;                    % 码元速率1K
dt=Ts/L;                    % 采样间隔
TotalT=M*Ts;                % 总时间
t=0:dt:TotalT-dt;           % 时间
TotalT2=(M/2)*Ts;           % 总时间2
t2=0:dt:TotalT2-dt;         % 时间2
Fs=1/dt;                    % 采样间隔的倒数即采样频率

%% 产生双极性波形
wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
wave=2*wave-1;              % 单极性变双极性
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=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)];
    end
end
x2 = I(fz,:);               % 将原来I的第一行复制L次，称为L*(M/2)的矩阵
I_lu = reshape(x2,1,L*(M/2));% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*(M/2))的矩阵
x3 = Q(fz,:);               % 将原来Q的第一行复制L次，称为L*(M/2)的矩阵
Q_lu = reshape(x3,1,L*(M/2));% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*(M/2))的矩阵
figure(1);                  % 绘制第1幅图
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(t,jidai,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('基带信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-1.1,1.1])   % 坐标范围限制

subplot(312);               % 窗口分割成3*1的，当前是第2个子图 
plot(t2,I_lu,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('I路信号波形');       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-1.1,1.1])  % 坐标范围限制

subplot(313);               % 窗口分割成3*1的，当前是第3个子图 
plot(t2,Q_lu,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('Q路信号波形');       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-1.1,1.1])  % 坐标范围限制
%% QPSK调制
fc=2000;                    % 载波频率2kHz       
zb1=cos(2*pi*fc*t2);        % 载波1
psk1=I_lu.*zb1;             % PSK1的调制 
zb2=sin(2*pi*fc*t2);        % 载波2
psk2=Q_lu.*zb2;             % PSK2的调制 
qpsk=psk1+psk2;             % QPSK的实现

figure(2);                  % 绘制第2幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(t2,qpsk,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('QPSK信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-2.1,2.1])  % 坐标范围限制

%% 信号经过高斯白噪声信道
tz=awgn(qpsk,20);           % 信号qpsk中加入白噪声，信噪比为SNR=20dB
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(t2,tz,'LineWidth',2);  % 绘制2ASK信号加入白噪声的波形
axis([0,TotalT2,-2.5,2.5]); % 坐标范围设置
title('通过高斯白噪声信道后的信号');% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 解调部分
figure(3);
tz1=tz.*zb1;                % 相干解调，乘以相干载波
subplot(211)                % 窗口分割成2*1的，当前是第1个子图 
plot(t2,tz1,'LineWidth',2)  % 绘制I路乘以相干载波后的信号
axis([0,TotalT2,-2.5,2.5]); % 设置坐标范围
title("I路乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

tz2=tz.*zb2;                % 相干解调，乘以相干载波
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t2,tz2,'LineWidth',2)  % 绘制Q路乘以相干载波后的信号
axis([0,TotalT2,-2.5,2.5]); % 设置坐标范围
title("Q路乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 加噪信号经过滤波器
% 低通滤波器设计
fp=2*Rb;                    % 低通滤波器截止频率，乘以2是因为下面要将模拟频率转换成数字频率wp=Rb/(Fs/2)
b=fir1(30, fp/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗滤波器
[h,w]=freqz(b, 1,512);      % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
lvbo1=fftfilt(b,tz1);       % 对信号进行滤波，tz1是等待滤波的信号，b是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b,tz2);       % 对信号进行滤波，tz2是等待滤波的信号，b是fir滤波器的系统函数的分子多项式系数
figure(4);                  % 绘制第4幅图  
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(w/pi*Fs/2,20*log(abs(h)),'LineWidth',2); % 绘制滤波器的幅频响应
title('低通滤波器的频谱');  % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签

subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t2,lvbo1,'LineWidth',2); % 绘制I路经过低通滤波器后的信号
axis([0,TotalT2,-2.1,2.1]);  % 设置坐标范围
title("I路经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(313)                % 窗口分割成3*1的，当前是第3个子图 
plot(t2,lvbo2,'LineWidth',2); % 绘制Q路经过低通滤波器后的信号
axis([0,TotalT2,-2.1,2.1]);  % 设置坐标范围
title("Q路经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 抽样判决
figure(5);
k=0;                        % 设置抽样限值
pdst1=1*(lvbo1>0);          % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0
pdst2=1*(lvbo2>0);          % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0
subplot(311)                % 窗口分割成3*1的，当前是第1个子图 
plot(t2,pdst1,'LineWidth',2)% 画出经过抽样判决后的信号
axis([0,TotalT2,-0.1,1.1]); % 设置坐标范用
title("I路经过抽样判决后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t2,pdst2,'LineWidth',2)% 画出经过抽样判决后的信号
axis([0,TotalT2,-0.1,1.1]); % 设置坐标范用
title("Q路经过抽样判决后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% I、Q路合并
I_zong = [];
Q_zong = [];
% 取码元的中间位置上的值进行判决
for j=L/2:L:(L*M/2)
    if pdst1(j)>0
        I_zong=[I_zong,1];
    else
        I_zong=[I_zong,-1];
    end
end
% 取码元的中间位置上的值进行判决
for k=L/2:L:(L*M/2)
    if pdst2(k)>0
        Q_zong=[Q_zong,1];
    else
        Q_zong=[Q_zong,-1];
    end
end
code = [];
% 将I路码元为最终输出的奇数位置码元，将Q路码元为最终输出的偶数位置码元
for n=1:M
    if mod(n, 2)~=0
        code = [code, I_zong((n+1)/2)];
    else
        code = [code, Q_zong(n/2)];
    end
end

x4=code(fz,:);             % 将原来code的第一行复制L次，称为L*M的矩阵
dout=reshape(x4,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵

subplot(313);              % 窗口分割成3*1的，当前是第3个子图 
plot(t,dout,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('I、Q路合并信号波形'); % 标题
xlabel('时间/s');          % x轴标签
ylabel('幅度');            % y轴标签
axis([0,TotalT,-1.1,1.1])  % 坐标范围限制
%% 绘制频谱
%% 信源频谱
figure(6)                   % 绘制第6幅图
T=t(end);                   % 时间
df=1/T;                     % 频谱分辨率
N=length(jidai);            % 采样长度
f=(-N/2:N/2-1)*df;          % 频率范围
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f,mf,'LineWidth',2);   % 绘制信源频谱波形
title("基带信号频谱");      % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% QPSK信号频谱
T2=t2(end);                 % 时间2
df2=1/T2;                   % 频谱分辨率2
N2=length(qpsk);            % 采样长度2
f2=(-N2/2:N2/2-1)*df2;      % 频率范围
sf=fftshift(abs(fft(qpsk)));% 对QPSK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(f2,sf,'LineWidth',2)   % 绘制QPSK调制信号频谱
title("QPSK信号频谱")       % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])   % 坐标范围限制

%% I路乘以相干载波后的频谱
mmf=fftshift(abs(fft(tz1)));% 对相干载波信号采用快速傅里叶变换并移到矩阵中心
figure(7)                   % 绘制第7幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f2,mmf,'LineWidth',2)  % 画出乘以相干载波后的频谱
title("I路乘以相干载波后的频谱")
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% Q路乘以相干载波后的频谱
mmf2=fftshift(abs(fft(tz2)));% 对相干载波信号采用快速傅里叶变换并移到矩阵中心
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(f2,mmf2,'LineWidth',2) % 画出乘以相干载波后的频谱
title("Q路乘以相干载波后的频谱")
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% 经过低通滤波后的频谱
figure(8);
dmf=fftshift(abs(fft(lvbo1)));%对低通滤波信号采用快速傅里叶变换并移到矩阵中心
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f2,dmf,'LineWidth',2)  % 画出经过低通滤波后的频谱
title("I路经过低通滤波后的频谱");
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签

dmf2=fftshift(abs(fft(lvbo2)));%对低通滤波信号采用快速傅里叶变换并移到矩阵中心
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(f2,dmf2,'LineWidth',2) % 画出经过低通滤波后的频谱
title("Q路经过低通滤波后的频谱");
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248

4.结果图

结果图中QPSK信号是经过信道，加了高斯白噪声的。
如果不想加噪声，把下面这行代码注释即可。

tz=awgn(qpsk,20);           % 信号qpsk中加入白噪声，信噪比为SNR=20dB
1

在这里插入图片描述

5.特点

随着进制数增加，抗噪声性能下降，但是传输速率提高，牺牲通信系统的可靠性，提高有效性。

6.加星座图的QPSK代码

clear all;                  % 清除所有变量
close all;                  % 关闭所有窗口
clc;                        % 清屏
%% 基本参数
M=20;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;                   % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;                    % 码元速率1K
dt=Ts/L;                    % 采样间隔
TotalT=M*Ts;                % 总时间
t=0:dt:TotalT-dt;           % 时间
TotalT2=(M/2)*Ts;           % 总时间2
t2=0:dt:TotalT2-dt;         % 时间2
Fs=1/dt;                    % 采样间隔的倒数即采样频率

%% 产生双极性波形
wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
wave=2*wave-1;              % 单极性变双极性
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=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)];
    end
end
x2 = I(fz,:);               % 将原来I的第一行复制L次，称为L*(M/2)的矩阵
I_lu = reshape(x2,1,L*(M/2));% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*(M/2))的矩阵
x3 = Q(fz,:);               % 将原来Q的第一行复制L次，称为L*(M/2)的矩阵
Q_lu = reshape(x3,1,L*(M/2));% 将刚得到的L*(M/2)矩阵，按列重新排列形成1*(L*(M/2))的矩阵
figure(1);                  % 绘制第1幅图
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(t,jidai,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('基带信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-1.1,1.1])   % 坐标范围限制

subplot(312);               % 窗口分割成3*1的，当前是第2个子图 
plot(t2,I_lu,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('I路信号波形');       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-1.1,1.1])  % 坐标范围限制

subplot(313);               % 窗口分割成3*1的，当前是第3个子图 
plot(t2,Q_lu,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('Q路信号波形');       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-1.1,1.1])  % 坐标范围限制
%% QPSK调制
fc=2000;                    % 载波频率2kHz       
zb1=cos(2*pi*fc*t2);        % 载波1
psk1=I_lu.*zb1;             % PSK1的调制 
zb2=sin(2*pi*fc*t2);        % 载波2
psk2=Q_lu.*zb2;             % PSK2的调制 
qpsk=psk1+psk2;             % QPSK的实现

figure(2);                  % 绘制第2幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(t2,qpsk,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('QPSK信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT2,-2.1,2.1])  % 坐标范围限制

%% 信号经过高斯白噪声信道
tz=awgn(qpsk,20);           % 信号qpsk中加入白噪声，信噪比为SNR=20dB
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(t2,tz,'LineWidth',2);  % 绘制2ASK信号加入白噪声的波形
axis([0,TotalT2,-2.5,2.5]); % 坐标范围设置
title('通过高斯白噪声信道后的信号');% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 解调部分
figure(3);
tz1=tz.*zb1;                % 相干解调，乘以相干载波
subplot(211)                % 窗口分割成2*1的，当前是第1个子图 
plot(t2,tz1,'LineWidth',2)  % 绘制I路乘以相干载波后的信号
axis([0,TotalT2,-2.5,2.5]); % 设置坐标范围
title("I路乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

tz2=tz.*zb2;                % 相干解调，乘以相干载波
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t2,tz2,'LineWidth',2)  % 绘制Q路乘以相干载波后的信号
axis([0,TotalT2,-2.5,2.5]); % 设置坐标范围
title("Q路乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 加噪信号经过滤波器
% 低通滤波器设计
fp=2*Rb;                    % 低通滤波器截止频率，乘以2是因为下面要将模拟频率转换成数字频率wp=Rb/(Fs/2)
b=fir1(30, fp/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗滤波器
[h,w]=freqz(b, 1,512);      % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
lvbo1=fftfilt(b,tz1);       % 对信号进行滤波，tz1是等待滤波的信号，b是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b,tz2);       % 对信号进行滤波，tz2是等待滤波的信号，b是fir滤波器的系统函数的分子多项式系数
figure(4);                  % 绘制第4幅图  
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(w/pi*Fs/2,20*log(abs(h)),'LineWidth',2); % 绘制滤波器的幅频响应
title('低通滤波器的频谱');  % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签

subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t2,lvbo1,'LineWidth',2); % 绘制I路经过低通滤波器后的信号
axis([0,TotalT2,-2.1,2.1]);  % 设置坐标范围
title("I路经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(313)                % 窗口分割成3*1的，当前是第3个子图 
plot(t2,lvbo2,'LineWidth',2); % 绘制Q路经过低通滤波器后的信号
axis([0,TotalT2,-2.1,2.1]);  % 设置坐标范围
title("Q路经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 抽样判决
figure(5);
k=0;                        % 设置抽样限值
pdst1=1*(lvbo1>0);          % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0
pdst2=1*(lvbo2>0);          % 滤波后的向量的每个元素和0进行比较，大于0为1，否则为0
subplot(311)                % 窗口分割成3*1的，当前是第1个子图 
plot(t2,pdst1,'LineWidth',2)% 画出经过抽样判决后的信号
axis([0,TotalT2,-0.1,1.1]); % 设置坐标范用
title("I路经过抽样判决后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t2,pdst2,'LineWidth',2)% 画出经过抽样判决后的信号
axis([0,TotalT2,-0.1,1.1]); % 设置坐标范用
title("Q路经过抽样判决后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% I、Q路合并
I_zong = [];
Q_zong = [];
% 取码元的中间位置上的值进行判决
for j=L/2:L:(L*M/2)
    if pdst1(j)>0
        I_zong=[I_zong,1];
    else
        I_zong=[I_zong,-1];
    end
end
% 取码元的中间位置上的值进行判决
for k=L/2:L:(L*M/2)
    if pdst2(k)>0
        Q_zong=[Q_zong,1];
    else
        Q_zong=[Q_zong,-1];
    end
end
code = [];
% 将I路码元为最终输出的奇数位置码元，将Q路码元为最终输出的偶数位置码元
for n=1:M
    if mod(n, 2)~=0
        code = [code, I_zong((n+1)/2)];
    else
        code = [code, Q_zong(n/2)];
    end
end

x4=code(fz,:);             % 将原来code的第一行复制L次，称为L*M的矩阵
dout=reshape(x4,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵

subplot(313);              % 窗口分割成3*1的，当前是第3个子图 
plot(t,dout,'LineWidth',2);% 绘制基带码元波形，线宽为2
title('I、Q路合并信号波形'); % 标题
xlabel('时间/s');          % x轴标签
ylabel('幅度');            % y轴标签
axis([0,TotalT,-1.1,1.1])  % 坐标范围限制
%% 绘制频谱
%% 信源频谱
figure(6)                   % 绘制第6幅图
T=t(end);                   % 时间
df=1/T;                     % 频谱分辨率
N=length(jidai);            % 采样长度
f=(-N/2:N/2-1)*df;          % 频率范围
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f,mf,'LineWidth',2);   % 绘制信源频谱波形
title("基带信号频谱");      % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% QPSK信号频谱
T2=t2(end);                 % 时间2
df2=1/T2;                   % 频谱分辨率2
N2=length(qpsk);            % 采样长度2
f2=(-N2/2:N2/2-1)*df2;      % 频率范围
sf=fftshift(abs(fft(qpsk)));% 对QPSK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(f2,sf,'LineWidth',2)   % 绘制QPSK调制信号频谱
title("QPSK信号频谱")       % 标题
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])   % 坐标范围限制

%% I路乘以相干载波后的频谱
mmf=fftshift(abs(fft(tz1)));% 对相干载波信号采用快速傅里叶变换并移到矩阵中心
figure(7)                   % 绘制第7幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f2,mmf,'LineWidth',2)  % 画出乘以相干载波后的频谱
title("I路乘以相干载波后的频谱")
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% Q路乘以相干载波后的频谱
mmf2=fftshift(abs(fft(tz2)));% 对相干载波信号采用快速傅里叶变换并移到矩阵中心
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(f2,mmf2,'LineWidth',2) % 画出乘以相干载波后的频谱
title("Q路乘以相干载波后的频谱")
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签
axis([-20000,20000,-inf,inf])% 坐标范围限制

%% 经过低通滤波后的频谱
figure(8);
dmf=fftshift(abs(fft(lvbo1)));%对低通滤波信号采用快速傅里叶变换并移到矩阵中心
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(f2,dmf,'LineWidth',2)  % 画出经过低通滤波后的频谱
title("I路经过低通滤波后的频谱");
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签

dmf2=fftshift(abs(fft(lvbo2)));%对低通滤波信号采用快速傅里叶变换并移到矩阵中心
subplot(212);               % 窗口分割成2*1的，当前是第2个子图 
plot(f2,dmf2,'LineWidth',2) % 画出经过低通滤波后的频谱
title("Q路经过低通滤波后的频谱");
xlabel('频率/Hz');          % x轴标签
ylabel('幅度');             % y轴标签

figure(9);                  % 绘制第9幅图
scatter(I,Q);               % 绘制星座图，如果没有4个点，是因为码元数量不够的原因
% 让坐标轴放在坐标原点
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254