通信原理与MATLAB(七):2FSK的调制解调_2fsk相干解调原理波形图

作者：Monodyee | 2024-04-25 00:56:24

踩

2fsk相干解调原理波形图

1. 2FSK的调制原理

2FSK调制原理如下图所示，基带码元d(t)中码元为1时，波形为频率为f1的高频载波；基带码元d(t)中码元为0时，波形为频率为f2的高频载波实现2FSK信号的调制，即基带码元和f1的高频正弦波相乘生成2ASK，基带码元的反码和f2的高频正弦波相乘生成第二个2ASK，两个2ASK相加得到2FSK。
在这里插入图片描述
波形图如下图所示

2. 2FSK的解调原理

2FSK的解调原理如下图所示，2FSK信号经过信道传输之后，分为上下两路经过带通滤波器变成两路2ASK信号，再和对应的载波相乘，然后经过低通滤波后抽样判决恢复出原始基带码元信号。
在这里插入图片描述

3. 2FSK的代码

clear all;                  % 清除所有变量
close all;                  % 关闭所有窗口
clc;                        % 清屏
%% 基本参数
M=10;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;                   % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;                    % 码元速率1K
dt=Ts/L;                    % 采样间隔
TotalT=M*Ts;                % 总时间
t=0:dt:TotalT-dt;           % 时间
Fs=1/dt;                    % 采样间隔的倒数即采样频率
%% 产生单极性波形
wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=reshape(x1,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵
%% 基带信号取反
jidai2=ones(1,length(jidai));% 产生和jidai相同长度的jidai2
% 基带信号取反即jidai为1的时候，jidai2为0;jidai为0的时候，jidai2为1
for n=1:length(jidai)
    if jidai(n)==1
        jidai2(n)=0;
    else
        jidai2(n)=1;
    end
end
%% 2FSK调制
fc1=2000;                   % 载波1频率2kHz       
zb1=cos(2*pi*fc1*t);        % 载波1信号
fsk1=jidai.*zb1;            % 基带信号1和载波1相乘 
fc2=10000;                  % 载波2频率10kHz       
zb2=cos(2*pi*fc2*t);        % 载波2信号
fsk2=jidai2.*zb2;           % 基带信号2和载波2相乘  
fsk=fsk1+fsk2;              % 2FSK的调制
figure(1);                  % 绘制第1幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形')% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK的波形 
title('2FSK信号波形')       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-1.1,1.1]);  % 坐标范围限制

%% 信号经过高斯白噪声信道
fsk=awgn(fsk,20);           % 信号fsk中加入白噪声，信噪比为SNR=20dB
figure(2);                  % 绘制第2幅图
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK信号加入白噪声的波形
axis([0,TotalT,-1.5,1.5]);  % 坐标范围设置
title('通过高斯白噪声信道后的信号');% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 带通滤波器设计
%% 第一个带通滤波器设计
fp1=2*(fc1-Rb);             % 第一个带通滤波器频率1
fs1=2*(fc1+Rb);             % 第一个带通滤波器频率2
b1=fir1(30, [fp1/Fs fs1/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h1,w1]=freqz(b1, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band1=fftfilt(b1,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,band1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器1后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 第二个带通滤波器设计
fp2=2*(fc2-Rb);             % 第二个带通滤波器频率1
fs2=2*(fc2+Rb);             % 第二个带通滤波器频率2
b2=fir1(30, [fp2/Fs fs2/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h2,w2]=freqz(b2, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band2=fftfilt(b2,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b2是fir滤波器的系统函数的分子多项式系数
subplot(313)                % 窗口分割成3*1的，当前是第3个子图 
plot(t,band2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器2后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 带通滤波器的频谱
figure(3);                  % 绘制第3幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(w1/pi*Fs/2,20*log(abs(h1)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第一个带通滤波器频谱');    % 标题

subplot(212);               % 窗口分割成2*1的，当前是第2个子图
plot(w2/pi*Fs/2,20*log(abs(h2)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第二个带通滤波器频谱');    % 标题

%% 解调部分
%% fsk信号经过第一个带通滤波后和载波1相乘
figure(4);                  % 绘制第4幅图
jt1=band1.*(-zb1);          % 相干解调，乘以相干载波
subplot(211)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt1,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% fsk信号经过第一个带通滤波后和载波1相乘
jt2=band2.*(-zb2);          % 相干解调，乘以相干载波
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt2,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 加噪信号经过滤波器
%% 低通滤波器设计
fp3=2*Rb;                   % 低通滤波器截止频率
b3=fir1(30, fp3/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗滤波器
[h3,w3]=freqz(b3, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
lvbo1=fftfilt(b3,jt1);      % 对信号进行滤波，jt1是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b3,jt2);      % 对信号进行滤波，jt2是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
figure(5);                  % 绘制第5幅图
subplot(311)                % 窗口分割成3*1的，当前是第1个子图 
plot(t,lvbo1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-0.1,1.1]);  % 设置坐标范围
title("经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

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

subplot(313);               % 窗口分割成3*1的，当前是第3个子图
plot(w3/pi*Fs/2,20*log(abs(h3)),'LineWidth',2); % 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('低通滤波器频谱');    % 标题

%% 抽样判决
pdst=1*(lvbo1>lvbo2);       % 上面支路信号大于下面支路信号判决为1，否则为0 
figure(6);                  % 绘制第6幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(t,pdst,'LineWidth',2) % 画出经过抽样判决后的信号
axis([0,TotalT,-0.1,1.1]); % 设置坐标范用
title("经过抽样判决后的信号")% 标题
xlabel('时间/s');          % x轴标签
ylabel('幅度');            % y轴标签

%% 绘制频谱
%% 2FSK信号频谱
T=t(end);                  % 时间
df=1/T;                    % 频谱分辨率
N=length(fsk);             % 采样长度
f=(-N/2:N/2-1)*df;         % 频率范围
sf=fftshift(abs(fft(fsk)));% 对2FSK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
figure(7)                  % 绘制第7幅图
subplot(211)               % 窗口分割成2*1的，当前是第1个子图 
plot(f,sf,'LineWidth',2)   % 绘制调制信号频谱
title("2FSK信号频谱")      % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 信源频谱
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf,'LineWidth',2)   % 绘制信源频谱波形
title("信源频谱")          % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

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

%% 上面支路和载波1相乘后频谱
figure(9);                 % 绘制第9幅图
mf4=fftshift(abs(fft(jt1)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(211)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf4,'LineWidth',2)  % 绘制信源频谱波形
title("上面支路和载波1相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 下面支路和载波2相乘后频谱
mf5=fftshift(abs(fft(jt2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf5,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路和载波2相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

%% 下面支路经过低通滤波后频谱
mf7=fftshift(abs(fft(lvbo2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf7,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路经过低通滤波后频谱");% 标题
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

4. 结果图

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

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

在这里插入图片描述

5. 特点

抗噪声能力一般，但是可用于多径信道中，抗多径能力强。

6. 代码改进

上述代码中判决不是在码元中间时刻判决的，效果差点意思，现在进行改进，在码元中间时刻判决，代码如下。
BFSK改进代码

clear all;                  % 清除所有变量
close all;                  % 关闭所有窗口
clc;                        % 清屏
%% 基本参数
M=10;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;                   % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;                    % 码元速率1K
dt=Ts/L;                    % 采样间隔
TotalT=M*Ts;                % 总时间
t=0:dt:TotalT-dt;           % 时间
Fs=1/dt;                    % 采样间隔的倒数即采样频率
%% 产生单极性波形
wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=reshape(x1,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵
%% 基带信号取反
jidai2=ones(1,length(jidai));% 产生和jidai相同长度的jidai2
% 基带信号取反即jidai为1的时候，jidai2为0;jidai为0的时候，jidai2为1
for n=1:length(jidai)
    if jidai(n)==1
        jidai2(n)=0;
    else
        jidai2(n)=1;
    end
end
%% 2FSK调制
fc1=2000;                   % 载波1频率2kHz       
zb1=cos(2*pi*fc1*t);        % 载波1信号
fsk1=jidai.*zb1;            % 基带信号1和载波1相乘 
fc2=10000;                  % 载波2频率10kHz       
zb2=cos(2*pi*fc2*t);        % 载波2信号
fsk2=jidai2.*zb2;           % 基带信号2和载波2相乘  
fsk=fsk1+fsk2;              % 2FSK的调制
figure(1);                  % 绘制第1幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形')% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK的波形 
title('2FSK信号波形')       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-1.1,1.1]);  % 坐标范围限制

%% 信号经过高斯白噪声信道
fsk=awgn(fsk,20);           % 信号fsk中加入白噪声，信噪比为SNR=20dB
figure(2);                  % 绘制第2幅图
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK信号加入白噪声的波形
axis([0,TotalT,-1.5,1.5]);  % 坐标范围设置
title('通过高斯白噪声信道后的信号');% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 带通滤波器设计
%% 第一个带通滤波器设计
fp1=2*(fc1-Rb);             % 第一个带通滤波器频率1
fs1=2*(fc1+Rb);             % 第一个带通滤波器频率2
b1=fir1(30, [fp1/Fs fs1/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h1,w1]=freqz(b1, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band1=fftfilt(b1,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,band1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器1后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 第二个带通滤波器设计
fp2=2*(fc2-Rb);             % 第二个带通滤波器频率1
fs2=2*(fc2+Rb);             % 第二个带通滤波器频率2
b2=fir1(30, [fp2/Fs fs2/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h2,w2]=freqz(b2, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band2=fftfilt(b2,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b2是fir滤波器的系统函数的分子多项式系数
subplot(313)                % 窗口分割成3*1的，当前是第3个子图 
plot(t,band2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器2后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 带通滤波器的频谱
figure(3);                  % 绘制第3幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(w1/pi*Fs/2,20*log(abs(h1)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第一个带通滤波器频谱');    % 标题

subplot(212);               % 窗口分割成2*1的，当前是第2个子图
plot(w2/pi*Fs/2,20*log(abs(h2)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第二个带通滤波器频谱');    % 标题

%% 解调部分
%% fsk信号经过第一个带通滤波后和载波1相乘
figure(4);                  % 绘制第4幅图
jt1=band1.*(-zb1);          % 相干解调，乘以相干载波
subplot(211)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt1,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% fsk信号经过第一个带通滤波后和载波1相乘
jt2=band2.*(-zb2);          % 相干解调，乘以相干载波
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt2,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 加噪信号经过滤波器
%% 低通滤波器设计
fp3=2*Rb;                   % 低通滤波器截止频率
b3=fir1(30, fp3/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗滤波器
[h3,w3]=freqz(b3, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
lvbo1=fftfilt(b3,jt1);      % 对信号进行滤波，jt1是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b3,jt2);      % 对信号进行滤波，jt2是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
figure(5);                  % 绘制第5幅图
subplot(311)                % 窗口分割成3*1的，当前是第1个子图 
plot(t,lvbo1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-0.1,1.1]);  % 设置坐标范围
title("经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

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

subplot(313);               % 窗口分割成3*1的，当前是第3个子图
plot(w3/pi*Fs/2,20*log(abs(h3)),'LineWidth',2); % 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('低通滤波器频谱');    % 标题

%% 抽样判决
pdst=1*(lvbo1>lvbo2);       % 上面支路信号大于下面支路信号判决为1，否则为0 
% 取码元中间时刻值为判决值
panjue=[];

for j=(L/2):L:(L*M)
    if pdst(j)>0
        panjue=[panjue,1];
    else
        panjue=[panjue,0];
    end
end
x2=panjue(fz,:);              % 将原来panjue的第一行复制L次，称为L*M的矩阵
panjue_zong=reshape(x2,1,L*M);% 将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵
figure(6);                  % 绘制第6幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(t,panjue_zong,'LineWidth',2) % 画出经过抽样判决后的信号
axis([0,TotalT,-0.1,1.1]); % 设置坐标范用
title("经过抽样判决后的信号")% 标题
xlabel('时间/s');          % x轴标签
ylabel('幅度');            % y轴标签

%% 绘制频谱
%% 2FSK信号频谱
T=t(end);                  % 时间
df=1/T;                    % 频谱分辨率
N=length(fsk);             % 采样长度
f=(-N/2:N/2-1)*df;         % 频率范围
sf=fftshift(abs(fft(fsk)));% 对2FSK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
figure(7)                  % 绘制第7幅图
subplot(211)               % 窗口分割成2*1的，当前是第1个子图 
plot(f,sf,'LineWidth',2)   % 绘制调制信号频谱
title("2FSK信号频谱")      % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 信源频谱
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf,'LineWidth',2)   % 绘制信源频谱波形
title("信源频谱")          % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

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

%% 上面支路和载波1相乘后频谱
figure(9);                 % 绘制第9幅图
mf4=fftshift(abs(fft(jt1)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(211)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf4,'LineWidth',2)  % 绘制信源频谱波形
title("上面支路和载波1相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 下面支路和载波2相乘后频谱
mf5=fftshift(abs(fft(jt2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf5,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路和载波2相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

%% 下面支路经过低通滤波后频谱
mf7=fftshift(abs(fft(lvbo2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf7,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路经过低通滤波后频谱");% 标题
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
249
250
251
252
253
254
255
256
257
258

7. BFSK误码率曲线

BFSK误码率曲线可以通过该链接获取BFSK误码率曲线

8. BFSK改进代码

该代码解决评论中提出的问题。

clear all;                  % 清除所有变量
close all;                  % 关闭所有窗口
clc;                        % 清屏
%% 基本参数
M=10;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;                   % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;                    % 码元速率1K
dt=Ts/L;                    % 采样间隔
TotalT=M*Ts;                % 总时间
t=0:dt:TotalT-dt;           % 时间
Fs=1/dt;                    % 采样间隔的倒数即采样频率
%% 产生单极性波形
wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=reshape(x1,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵
%% 基带信号取反
jidai2=ones(1,length(jidai));% 产生和jidai相同长度的jidai2
% 基带信号取反即jidai为1的时候，jidai2为0;jidai为0的时候，jidai2为1
for n=1:length(jidai)
    if jidai(n)==1
        jidai2(n)=0;
    else
        jidai2(n)=1;
    end
end
%% 2FSK调制
fc1=10000;                   % 载波1频率10kHz       
zb1=cos(2*pi*fc1*t);        % 载波1信号
fsk1=jidai.*zb1;            % 基带信号1和载波1相乘 
fc2=20000;                  % 载波2频率20kHz       
zb2=cos(2*pi*fc2*t);        % 载波2信号
fsk2=jidai2.*zb2;           % 基带信号2和载波2相乘  
fsk=fsk1+fsk2;              % 2FSK的调制
figure(1);                  % 绘制第1幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图 
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形')% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK的波形 
title('2FSK信号波形')       % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-1.1,1.1]);  % 坐标范围限制

%% 信号经过高斯白噪声信道
fsk=awgn(fsk,40);           % 信号fsk中加入白噪声，信噪比为SNR=20dB
figure(2);                  % 绘制第2幅图
subplot(311);               % 窗口分割成3*1的，当前是第1个子图 
plot(t,fsk,'LineWidth',2);  % 绘制2FSK信号加入白噪声的波形
axis([0,TotalT,-1.5,1.5]);  % 坐标范围设置
title('通过高斯白噪声信道后的信号');% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 带通滤波器设计
%% 第一个带通滤波器设计
fp1=2*(fc1-Rb);             % 第一个带通滤波器频率1
fs1=2*(fc1+Rb);             % 第一个带通滤波器频率2
b1=fir1(30, [fp1/Fs fs1/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h1,w1]=freqz(b1, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band1=fftfilt(b1,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
subplot(312)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,band1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器1后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 第二个带通滤波器设计
fp2=2*(fc2-Rb);             % 第二个带通滤波器频率1
fs2=2*(fc2+Rb);             % 第二个带通滤波器频率2
b2=fir1(30, [fp2/Fs fs2/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗带通滤波器
[h2,w2]=freqz(b2, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
band2=fftfilt(b2,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b2是fir滤波器的系统函数的分子多项式系数
subplot(313)                % 窗口分割成3*1的，当前是第3个子图 
plot(t,band2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器2后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
%% 带通滤波器的频谱
figure(3);                  % 绘制第3幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(w1/pi*Fs/2,20*log(abs(h1)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第一个带通滤波器频谱');    % 标题

subplot(212);               % 窗口分割成2*1的，当前是第2个子图
plot(w2/pi*Fs/2,20*log(abs(h2)),'LineWidth',2);% 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('第二个带通滤波器频谱');    % 标题

%% 解调部分
%% fsk信号经过第一个带通滤波后和载波1相乘
figure(4);                  % 绘制第4幅图
jt1=band1.*(-zb1);          % 相干解调，乘以相干载波
subplot(211)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt1,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% fsk信号经过第一个带通滤波后和载波1相乘
jt2=band2.*(zb2);          % 相干解调，乘以相干载波
subplot(212)                % 窗口分割成2*1的，当前是第2个子图 
plot(t,jt2,'LineWidth',1)   % 绘制乘以相干载波后的信号
axis([0,TotalT,-1.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%% 加噪信号经过滤波器
%% 低通滤波器设计
fp3=2*Rb;                   % 低通滤波器截止频率
b3=fir1(30, fp3/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
% fir1函数三个参数分别是阶数，数字截止频率，滤波器类型
% 这里是生成了30阶(31个抽头系数)的矩形窗滤波器
[h3,w3]=freqz(b3, 1,512);   % 生成fir滤波器的频率响应
% freqz函数的三个参数分别是滤波器系统函数的分子多项式的系数，分母多项式的系数(fir滤波器分母系数为1)和采样点数(默认)512
lvbo1=fftfilt(b3,jt1);      % 对信号进行滤波，jt1是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b3,jt2);      % 对信号进行滤波，jt2是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
figure(5);                  % 绘制第5幅图
subplot(311)                % 窗口分割成3*1的，当前是第1个子图 
plot(t,lvbo1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,TotalT,-0.1,1.1]);  % 设置坐标范围
title("经过低通滤波器后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

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

subplot(313);               % 窗口分割成3*1的，当前是第3个子图
plot(w3/pi*Fs/2,20*log(abs(h3)),'LineWidth',2); % 绘制滤波器的幅频响应
xlabel('频率/Hz');          % x轴标签
ylabel('幅度/dB');          % y轴标签
title('低通滤波器频谱');    % 标题

%% 抽样判决
pdst=1*(lvbo1>lvbo2);       % 上面支路信号大于下面支路信号判决为1，否则为0 
% 取码元中间时刻值为判决值
panjue=[];

for j=(L/2):L:(L*M)
    if pdst(j)>0
        panjue=[panjue,1];
    else
        panjue=[panjue,0];
    end
end
x2=panjue(fz,:);              % 将原来panjue的第一行复制L次，称为L*M的矩阵
panjue_zong=reshape(x2,1,L*M);% 将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵
figure(6);                  % 绘制第6幅图
subplot(211);               % 窗口分割成2*1的，当前是第1个子图
plot(t,jidai,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形');      % 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(t,panjue_zong,'LineWidth',2) % 画出经过抽样判决后的信号
axis([0,TotalT,-0.1,1.1]); % 设置坐标范用
title("经过抽样判决后的信号")% 标题
xlabel('时间/s');          % x轴标签
ylabel('幅度');            % y轴标签

%% 绘制频谱
%% 2FSK信号频谱
T=t(end);                  % 时间
df=1/T;                    % 频谱分辨率
N=length(fsk);             % 采样长度
f=(-N/2:N/2-1)*df;         % 频率范围
sf=fftshift(abs(fft(fsk)));% 对2FSK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
figure(7)                  % 绘制第7幅图
subplot(211)               % 窗口分割成2*1的，当前是第1个子图 
plot(f,sf,'LineWidth',2)   % 绘制调制信号频谱
title("2FSK信号频谱")      % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 信源频谱
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf,'LineWidth',2)   % 绘制信源频谱波形
title("信源频谱")          % 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

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

%% 上面支路和载波1相乘后频谱
figure(9);                 % 绘制第9幅图
mf4=fftshift(abs(fft(jt1)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(211)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf4,'LineWidth',2)  % 绘制信源频谱波形
title("上面支路和载波1相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

%% 下面支路和载波2相乘后频谱
mf5=fftshift(abs(fft(jt2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf5,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路和载波2相乘后频谱");% 标题
xlabel('频率/Hz');         % x轴标签
ylabel('幅度');            % y轴标签

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

%% 下面支路经过低通滤波后频谱
mf7=fftshift(abs(fft(lvbo2)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(212)               % 窗口分割成2*1的，当前是第2个子图 
plot(f,mf7,'LineWidth',2)  % 绘制信源频谱波形
title("下面支路经过低通滤波后频谱");% 标题
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
249
250
251
252
253
254
255
256
257
258
259
260

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Monodyee/article/detail/482461