通信原理实验_通信原理课内实验代码-优快云博客

本文链接：https://blog.youkuaiyun.com/m0_51077616/article/details/125400170

本文详细记录了一次通信原理实验，涵盖了AMI和HDB3编码的原理与实现，以及2ASK、2FSK、2PSK和2DPSK的调制解调过程。实验结果显示，HDB3编码有效地减少了直流分量，增强了抗扰能力。在不同信噪比下，2PSK和2DPSK的误码率曲线展示了其在高信噪比时的优越性能。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

通信原理实验报告

有具体实验需求可私聊定制

这篇文章记录了一下我的通信原理实验，文末有完整代码，但是最后的绘制误码率曲线部分我是直接进行散点图绘制，因此会跑大量数据，可以根据自己的需求减少点数通过插值等方法也可以获得平滑的曲线。中间的图片排版顺序因为导入原因有点顺序问题，也懒得改了

仿真报告1

实验目的

$M A T L A B$ 生成不同占空比单双极性方波，分析其频谱；
$M A T L A B$ 仿真进行 $A M I$ ， $HD B 3$ 编码；

实验原理

单极性方波

特点：极性单一、有直流分量和低频分量

双极性方波

优点：无直流分量(等概)、抗扰能力较强

$A M I$ 编码

不含直流分量，低频成分少
·三电平
编译码电路简单，有宏观自检能力
AMI编码缺点：信码有长连0串时，难以获取定时信息
AMI编码应用： $PCM 24$ 路基群(北美系列) $1.544 M b / s$ 的线路码型

$HD B 3$ 编译码

三阶高密度双（ $HD B 3$ 码）是一种适用于基带传输的编码方式，它是为了克服 $A M I$ 码的缺点而出现的，具有能量分散，抗破坏性强等特点。

编码规则

连0的个数不超过3时，规则与 $A M I$ 相同，即0不变，1变为-1、+1交替；
若连0的个数超过3，则将每4个0看作一小节，定义为 $B 00 V$ ， $B$ 可以是-1、0、+1， $V$ 可以是-1、+1；
B和V具体值满足以下条件： $V$ 和前面相邻非0符号极性相同；不看 $V$ 时极性交替； $V$ 与 $V$ 之间极性交替；
一般第一个 $B$ 取0，第一个非0符取-1。在 $V$ 与 $V$ 之间如果出现偶数个B时，应在后一个V字节补一个 $B^{'}$ ，定义为 $B^{'} 00 V$ ， $B^{'}$ 与前面相邻的 $B$ 之间符号极性相反，这个字节的 $V$ 和 $B^{'}$ 符号极性相同。由于 $V$ 会破坏极性交替的规律、 $B$ 有3种变化以满足规则，所以 $V$ 称为破坏脉冲、 $B$ 称为调节脉冲， $B 00 V$ 称为取代节、破坏节。

编码特点

由 $HD B 3$ 码确定的基带信号无直流分量，且只有很小的低频分量;
$HD B 3$ 中连0串的数目至多为3个，易于提取定时信号。
编码规则复杂，但译码较简单。

译码规则

译码的本质即把原来的取代节（4个连零）找到即可，若3连“0”前后非零脉冲同极性，则将最后一个非零元素译为零，如+1000+1 就应该译成“10000”，否则不用改动；若2连 “0”前后非零脉冲极性相同，则两零前后都译为零，如-100-1，就应该译为0000，否则也不用改动.
将所有的-1变换成+1后，就可以得到原消息码。

应用

数字基带信号的传输是数字通信系统的重要组成部分。在数字通信中，有些场合可不经过载波调制和解调过程，而对基带信号进行直接传输。采用AMI码的信号交替反转，有可能出现四连零现象，这不利于接收端的定时信号提取。而 $HD B 3$ 码因其无直流成份、低频成份少和连0个数最多不超过三个等特点，而对定时信号的恢复十分有利。

实验结果

单双极性归零码功率谱密度

在这里插入图片描述

可见占空比越大，主瓣带宽越宽，旁瓣越少，符合理想状态

$A M I$ 码 $HD B 3$ 码编码译码

在这里插入图片描述

$HD B 3$ 编码不会产生连续相同码元，减少直流分量

仿真报告二

实验目的

仿真 $2 A S K$ / $2 FS K$ 调制波形，分析其频谱
仿真对比 $A W GN$ 信道下 $2 PS K$ 和 $2 D PS K$ 的误码率曲线

实验原理

$2 A S K$

$A S K$ 即“振幅键控”， $2 A S K$ 是二进制振幅键控。这是正弦载波的幅度随数字基带信号而变化的数字调制。当数字基带信号为二进制时，则为二进制振幅键控。设发送的二进制符号序列由 $0$ ， $1$ 序列组模拟信号源调制器信道解调器受信者噪声源调制器信道解调器基带信号输入噪声源基带信号输出成，发送 $0$ 符号的概率为 $P$ ，发送 $1$ 符号的概率为 $1 - P$ ，且相互独立。
该二进制符号序列可表示为：

$s(t)=\sum_{n}a_{n}g(t-nT_{s})$

则二进制振幅键控信号可表示为:

$e{2ASK}(t)=\sum_{n}a_{n}g(t-nT_{s})cos\omega_{n}t$

$2 A S K$ 信号的产生方式通常有两种：模拟调制法和键控法。
$2 A S K$ 信号与模拟调制中的 $A M$ 信号类似。所以，对 $2 A S K$
信号也能够采用非相干解调（包络检波法）和相干解调（同步检测法）

$2 FS K$

$2 FS K$ 为二进制数字频率调制即二进制频移键控，用载波的频率来传送数字信息，即用所传送的数字信息控制载波的频率。 $2 FS K$ 信号便是符号“ $0$ ”对应于载频$ f1$，而符号“ $1$ ”对应于载频 $f 2$ 即与 $f 1$ 不同的另一载频的已调波形，而且 $f 1$ 与 $f 2$ 之间的改变是瞬间的。传“ $0$ ”信号时，发送频率为 $f 1$ 的载波；传“ $1$ ”信号时，发送频率为 $f 2$ 的载波。可见， $FS K$ 是用不同频率的载波来传递数字消息的。
一般来说，其信号产生有两种方法，即频率键控法和直接调频法。
频率键控法：两个分别产生正弦振荡的独立振荡器经由数字基带信号控制的电子开关后，选出的高频振荡信号就是 $FS K$ 调制信号。
直接调频法是利用数字基带信号直接控制载频振荡器的振荡频率。与键控法调频相比较，它产生的信号频率稳定性比键控法产生的信号差，且存在过渡频率。
在接收端，信号的解调方法有两种，一种为相干解调法，另一种叫非相干解调法也叫包络检波法。
相干载波与原调制的载波信号必须同频同相，理论上来说虽然在信号中确实存在着载波分量，但是由于提取载波分量的过程需要加上额外的电路，会给设备增加复杂度，因此，一般情况下均采用非相干解调的方式还原信号。幅移键控调制的方式出现较早，实现虽然容易，但相对于其他方式来说抗干扰能力不强，因此在实际中不常使用。
非相干解调首先将得到的信号进行带通滤波后滤除载波频率以外的噪声以及干扰，使得信号可以完整的通过，再经过全波整流器输出正极端的包络曲线，然后经过低通滤波器或者整流模块输出基带包络信号，再经过抽样判决器输出基带二进制信号。包络检波各个部分最终输出的波形在时间上相对于原基带二进制信号有一定的延时，这是硬件部分进行信号处理时无法避免的，在信号速率不大的情况下这种延时可以忽略。

$2 PS K$

$2 PS K$ 即为二进制相移键控，它用两个初相相隔为180的载波来传递二进制信息。所以也被称为BPSK。

相移键控是利用载波的相位变化来传递数字信息，而振幅和频率保持不变。

$2 PS K$ 信号的时域表达式为：

$e_{2psk(t)=Acos(\omega_{c}t+\phi_{n})}$

由于表示信号的两种码元的1波形相同，极性相反，故 $2 PS K$ 信号一般可以表述为一个双极性全占空矩形脉冲序列与一个正弦载波相乘，即

$e_{2psk(t)}=\sum_{n}a_{n}g(t-nT^{s}cos\omega_{c}t)$

发送为二进制符号“0”时（取+1,），取0相位；发送为二进制符号“1”时（取-1,），取Π相位。这种以载波的不同相位直接去表示相应二进制数字信号的调制方式，称为二进制绝对相移方式。
$2 PS K$ 的功率谱密度为：

$P_{2psk(f)}=\frac14[P_{s(f+f_{c})}P_{s}(f-f_{c})]$

$2 PS K$ 信号的调制方法分为模拟调制方法和数字调制方法。
$2 PS K$ 信号的解调通常采用相干解调法。假设相干载波的基准相位与 $2 PS K$ 信号的调制载波的基准相位一致（通常默认为0相位）。但是，由于在 $2 PS K$ 的载波恢复过程中存在着的相位模糊（，即恢复的本地载波与所需的相干载波可能同相，也有可能相反，这种相位关系的不确定性将会造成解调出的数字基带信号与发送的数字基带信号正好相反，即“1”变为“0”，“0”变为“1”，判决器输出数字信号全部出错。这种现象称为 $2 PS K$ 方式的“倒 ”现象或“反相工作”。这也是 $2 PS K$ 方式在实际中很少采用的主要原因。另外，在随机信号码元序列中，信号波形有可能出现长时间连续正弦波形，致使在接收端无法辨识信号码元的起止时刻。为了解决上述问题，可以采用差分相移键控（DPSK）体制。

$2 D PS K$

$2 PS K$ 信号中，相位的变化时以未调载波的相位作为参考基准的。又称为绝对相移。因此具有 $0$ 、 $\pi$ 模糊性，导致借条过程中出现“反相工作”现象，恢复出的数字信号“ $1$ ”，“ $0$ ”倒置，从而致使 $2 PS K$ 难以实际应用。为了克服此缺点，提出了二进制差分相移键控( $2 D PS K$ )方式。
$2 D PS K$ 是利用前后相邻码元的载波相对相位变化传递数字信息，所以又称为相对相移键控。假设 $\delta\phi$ 为当前码元的载波相位差，可定义一种数字信息与 $\delta\phi$ 之间的关系为：

$\Delta\phi = \left\{ \begin{array}{lr} 0 & : 表示数字信息“0”\\ \pi & : 表示数字信息“1” \end{array} \right.$

二进制数字信息	1	1	0	1	0	0	1	1	0
$2 D PS K$ 信号相位（ $0$ ）	$\pi$	$0$	$0$	$\pi$	$\pi$	$\pi$	$0$	$\pi$	$\pi$
或（ $\pi$ ）	$0$	$\pi$	$\pi$	$0$	$0$	$0$	$\pi$	$0$	$0$

(a)相对码		1	1	0	1
(b)绝对码	0	1	0	0	1

在接收端只能采用相干解调。 $2 D PS K$ 信号先经过带通滤波器，滤除调制信号频带以外的在信道中混入的噪声，此后该信号分为两路，一路延时一个码元的时间后与另一路的信号相乘，再经过低通滤波器去除高频成分，得到包含基带信号的低频信号，将其送入抽样判决器中进行抽样判决，抽样判决器的输出即为原基带信号。

实验结果

$2 A S K$ 调制解调

在这里插入图片描述

$2 FS K$ 调制

在这里插入图片描述

$2 PS K$ 调制解调

在这里插入图片描述

上图为信噪比分别在 $- 15 d B$ ， $0 d B$ ， $15 d B$ 时的调制波形，解调波形，信源以及译码序列

$2 D PS K$ 调制解调

在这里插入图片描述

上图分别为信噪比在 $- 15 d B$ ， $0 d B$ ， $15 d B$ 时的信源序列，译码序列，调制波形以及解调波形

误码率

在这里插入图片描述

第一个图为误码率在对数坐标系下的曲线，第二个图为误码率在直角坐标系下的曲线，可以看出在误码率水平较高的时候并不能看出两种编码的误码率关系，但当误码率逐渐降低的时候，两种调制方式的误码率逐渐趋近于 $2 PS K$ 误码率为 $2 D PS K$ 误码率的2倍，同样，误码率也和滤波器性能有关。

附录

双极性归零码

Ts = 1;             % 码元周期
fs = 1000;          % 单个码元抽样点数
dt = Ts/ fs;        % 抽样时间间隔
N = 100;            % 码元数
t=(0: N*fs-1)/fs;   % 时间刻度

nfft = fs*N         % fft长度
RanNum = randi([0, 1], 1, N) % 随机01序列

duty = 25;          % 占空比
gt1 = [ones(1, fs*duty/100), zeros(1, fs*(1-duty/100))]; % 定义每个码元
gt2 = -gt1;

se = []; % 初始化
for i = 1 : N
    if RanNum(i) == 1
        se = [se,gt1];
    else
        se = [se,gt2];
end
end

fft_se = fftshift(fft(se)); 
PE = 10 * log10(abs(fft_se) .^ 2 / (N * Ts));
PEL = (-length(fft_se) / 2 : length(fft_se) / 2 - 1) / 100*(1/Ts);
subplot(221)
% plot(t, se); grid on; axis([0 20 -1.5 1.5]);
plot(PEL, PE); grid on; axis([-20 20 -50 75]); hold on;
subtitle("占空比25%")

$A M I$ , $HD B 3$ 编码译码

% 初始化消息码
N = 22;% N = 100;        % 码元数
num_AMI = 0;        % 初始化计数器
num_HDB3 = 0;
ran_num = [1,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1]
% ran_num = randi([0, 1], 1, N); % 随机01串
ran_num_ori = ran_num;
ran_num_AMI = ran_num;
ran_num_HDB3 = ran_num;
num_0 = 0; % 0计数器
num_V = 0; % V计数器

% AMI编码
for i = 1:N
    if  ran_num_AMI(i) == 1
        num_AMI = num_AMI+1;
        ran_num_AMI(i) = (-1)^num_AMI*ran_num_AMI(i);
    end
end

% AMI解码
decode_AMI = ran_num_AMI;
for j = 1:N
    if  abs(decode_AMI(j)) == 1
        decode_AMI(j) = abs(decode_AMI(j));
    end
end

% HDB3编码
for j = 1:N
    if  ran_num_HDB3(j) == 1
        num_HDB3 = num_HDB3+1;
        ran_num_HDB3(j) = (-1)^num_HDB3*ran_num_HDB3(j);
        flag = ran_num_HDB3(j); % 上一个非0数作为标志
        num_0 = 0;        
    else if ran_num_HDB3(j) == 0            
        num_0 = num_0+1;
        if num_0 == 4
            num_V = num_V+1;
            num_0 = 0; % 连续四个0,0计数器清零
            if num_V == 1
                ran_num_HDB3(j) = 1*flag;
                flag_V = ran_num_HDB3(j);
                flag = flag_V;
            else
                ran_num_HDB3(j) = -flag_V; % 破坏脉冲
                flag_V = ran_num_HDB3(j);
                if ran_num_HDB3(j) ~= flag
                    ran_num_HDB3(j-3) = ran_num_HDB3(j); % 调节脉冲
                end
                flag = flag_V;
            end            
        end
    end
    end
end

% HDB3解码
decode_HDB3 = ran_num_HDB3;
count = 0;
for j = 1:N  
    if ran_num_HDB3(j) == 0 
        count  = count + 1;
        if  count == 3  %若3连“0”前后非零脉冲同极性，则将最后一个非零元素译为零
            if decode_HDB3(j + 1) * decode_HDB3(j - 3) == 1
                decode_HDB3(j + 1) = 0;
            end
        end
        if count == 2 %若2连 “0”前后非零脉冲极性相同，则两零前后都译为零
            if decode_HDB3(j + 1) * decode_HDB3(j  - 2) == 1
                decode_HDB3(j  - 2) = 0;
                decode_HDB3(j + 1) = 0;
            end
        end
    else
         count = 0;   
    end
end
decode_HDB3 = abs(decode_HDB3);

% 绘图区
disp('消息码：'); disp(ran_num_ori);
disp('AMI码：'); disp(ran_num_AMI);
disp('HDB3码：'); disp(ran_num_HDB3);
figure;
subplot(3, 1, 1); stairs(0:(N-1), ran_num_ori); grid on; axis([0 21 -2 2]); title('消息码');
subplot(3, 1, 2); stairs(0:(N-1), ran_num_AMI); grid on; axis([0 21 -2 2]); title('AMI编码');
subplot(3, 1, 3); stairs(0:(N-1), ran_num_HDB3); grid on; axis([0 21 -2 2]); title('HDB3编码');
figure;
subplot(3, 1, 1); stairs(0:(N-1), ran_num_ori); grid on; axis([0 21 -2 2]); title('消息码');
subplot(3, 1, 2); stairs(0:(N-1), decode_AMI); grid on; axis([0 21 -2 2]); title('AMI解码');
subplot(3, 1, 3); stairs(0:(N-1), decode_HDB3); grid on; axis([0 21 -2 2]); title('HDB3解码');

$2 A S K$

clc
clear
N = 1000                        % 码元数
fs = 100                        % 采样频率
source = randi([0, 1], 1, N);   % 信源
T = 1;                          % 单个码元持续时间
rate = 4;                       % 载波频率与码元频率之比
fc = rate*(1/T);                % 载波频率
t = (0:N*fs*T-1)/fs;            % 时间刻度
L = fs*T;
t_decode = (0:(N*fs*T+L/2-1))/fs;
decode_2ASK = zeros(1, N);

carrier = sin(2*pi*fc*t);       % 载波
s = reshape((repmat(source', 1, (fs*T)))', 1, N*fs*T);
ASK_res = carrier.*s;


% 计算信号的傅里叶变换
Y = fft(ASK_res);
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/(N*fs*T)); 
P1 = P2(1:N*fs*T/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% 定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，频率和幅值与理论相同
f = fs*(0:(N*fs*T/2))/(N*fs*T);


s_awgn = awgn(ASK_res, 15);
x = abs(s_awgn);    % 加噪，全波整流
Rp=5;               % 通带纹波系数
Rs=6;               % 阻带衰减系数
Wp=4;               % 通带截止频率
Fs=3;               % 阻带起始频率


Wp=Wp/(fs/2);    % 通带截止频率
Fs=Fs/(fs/2);    % 阻带起始频率
 
% 计算滤波器的最低阶数n和3dB截止频率Wn
[n,Wn]=buttord(Wp,Fs,Rp,Rs);
% 计算系统函数H(z)的分子、分母多项式系数
[b,a]=butter(n,Wn);
y=filter(b,a,x);      % filter函数实现IIR滤波器的直接形式


for i = 1:N
    if y(L*i-L/2)>=0.5
        decode_2ASK(i)=1;
    else
        decode_2ASK(i)=0;
    end
end
source;
decode_2ASK;
de_2ASK = [zeros(1, L/2), reshape(repmat(decode_2ASK', 1, L)', 1, N*L)];
k = xor(source, decode_2ASK);
error = length(find(k(1,:)==1))/N;

disp('误码率：')
disp(error)

df = 0.01;
ts = T/fs;
[s_2askf,~,~,f_]=T2F(ASK_res,ts,df,fs);
plot(f_,10*log10(abs(fftshift(s_2askf).^2/length(f_))));

figure;
subplot(4, 1, 1); plot(t, s); axis([0, 10, -0.5, 1.5]); title("信源")
subplot(4, 1, 2); plot(t, ASK_res); axis([0,10,-1.5,1.5]); title("调制信号")
subplot(4, 1, 3); plot(f,P1); xlabel('f (Hz)'); title('X(t)的单边幅值谱')
subplot(4, 1, 4); plot(f_,10*log10(abs(fftshift(s_2askf).^2/length(f_)))); axis([-50,50,-200,100]); title('功率谱密度')

figure;
subplot(4,1,1);plot(t, s_awgn);axis([0,10,-1.5,1.5]);
subplot(4,1,2);plot(t, x);axis([0,10,-0.5,1.5]);
subplot(4,1,3);plot(t, y);axis([0,10,-0.5,1.5]); 
subplot(4,1,4);plot(t_decode, de_2ASK);axis([0,10,-0.5,1.5]);

$2 FS K$

clc
clear
N = 10000                       % 码元数
fs = 900000                      % 采样频率
source = randi([0, 1], 1, N);   % 信源
T = 1/9000;                      % 单个码元持续时间
rate1 = 2;                      % 载波频率与码元频率之比
rate2 = 4;
fc1 = rate1*(1/T);              % 载波频率
fc2 = rate2*(1/T);
t = (0:N*fs*T-1)/fs;            % 时间刻度
L = fs*T;
t_decode = (0:(N*L+L/2-1))/fs;
decode_2FSK = zeros(1, N);

carrier1 = sin(2*pi*fc1*t); % 载波
carrier2 = sin(2*pi*fc2*t); 
s = reshape((repmat(source', 1, L))', 1, N*L);
ASK_res = carrier1.*s+carrier2.*(~s);


%计算信号的傅里叶变换
Y = fft(ASK_res);
%计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/(N*L)); 
P1 = P2(1:N*L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
%定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，频率和幅值与理论相同
f = fs*(0:(N*L/2))/(N*L);


s_awgn = awgn(ASK_res, 15);
x = abs(s_awgn);    % 加噪，全波整流
Rp=8;               % 通带纹波系数
Rs=10;              % 阻带衰减系数
% Wp=16000;               % 通带截止频率
% Fs=12000;               % 阻带起始频率
Wp=10000;             
Fs=8000;             
Wp1=[12000 24000];               
Fs1=[16000 20000];               
Wp2=[30000 42000];               
Fs2=[34000 38000];               

Wp=Wp/(fs/2);    % 通带截止频率
Fs=Fs/(fs/2);    % 阻带起始频率
Wp1=Wp1/(fs/2);
Fs1=Fs1/(fs/2);
Wp2=Wp2/(fs/2);
Fs2=Fs2/(fs/2);
 
% 计算滤波器的最低阶数n和3dB截止频率Wn
[n,Wn]=buttord(Wp,Fs,Rp,Rs);
[n1,Wn1]=buttord(Wp1,Fs1,Rp,Rs);
[n2,Wn2]=buttord(Wp2,Fs2,Rp,Rs);
% 计算系统函数H(z)的分子、分母多项式系数
[b,a]=butter(n,Wn);
[b1,a1]=butter(n1,Wn1);
[b2,a2]=butter(n2,Wn2);
y=filter(b,a,x);      % filter函数实现IIR滤波器的直接形式
de_car1=filter(b1,a1,s_awgn);
de_car2=filter(b2,a2,s_awgn);
de_car1_abs = abs(de_car1);
de_car2_abs = abs(de_car2);
y1 = filter(b, a, de_car1_abs);
y2 = filter(b, a, de_car2_abs);


figure;
subplot(4,1,1);plot(t, de_car1);axis([0,10*T,-1.5,1.5]);
subplot(4,1,2);plot(t, de_car2);axis([0,10*T,-1.5,1.5]);
subplot(4,1,3);plot(t, de_car1_abs);axis([0,10*T,-1.5,1.5]);
subplot(4,1,4);plot(t, de_car2_abs);axis([0,10*T,-1.5,1.5]);

for i = 1:N
    if y1(L*i-L/2)>y2(L*i-L/2)
        decode_2FSK(i)=1;
    else
        decode_2FSK(i)=0;
    end
end
source;
decode_2FSK;
de_2FSK = [zeros(1, L/2), reshape(repmat(decode_2FSK', 1, L)', 1, N*L)];
k = xor(source, decode_2FSK);
length(find(k(1,:)==1))
error = length(find(k(1,:)==1))/N;

disp('误码率：')
disp(error)

figure;
subplot(3, 1, 1); plot(t, s); axis([0, 10*T, -0.5, 1.5])
subplot(3, 1, 2); plot(t, ASK_res); axis([0,10*T,-1.5,1.5])
subplot(3,1,3); plot(f,P1); xlabel('f (Hz)'); title('X(t)的单边幅值谱')

figure;
subplot(4,1,1);plot(t, s_awgn);axis([0,10*T,-1.5,1.5]);
subplot(4,1,2);plot(t, x);axis([0,10*T,-0.5,1.5]);
subplot(4,1,3);plot(t, y1);axis([0,10*T,-0.5,1.5]);
subplot(4,1,4);plot(t, y2);axis([0,10*T,-0.5,1.5]);

figure;
subplot(2, 1, 1); plot(t, s); axis([0, 10*T, -0.5, 1.5])
subplot(2,1,2);plot(t_decode, de_2FSK);axis([0,10*T,-0.5,1.5]);

$2 PS K$

clc
clear
N = 10000                       % 码元数
fs = 900000                     % 采样频率
source = randi([0, 1], 1, N);   % 信源
AMI = source
AMI(AMI==0) = -1

T = 1/9000;                     % 单个码元持续时间
rate = 2;                       % 载波频率与码元频率之比
fc = rate*(1/T);                % 载波频率
t = (0:N*fs*T-1)/fs;            % 时间刻度
L = fs*T;
t_decode = (0:(N*fs*T+L/2-1))/fs;
decode_2ASK = zeros(1, N);

carrier = sin(2*pi*fc*t);       % 载波
Ori = reshape((repmat(source', 1, (fs*T)))', 1, N*fs*T);
s = reshape((repmat(AMI', 1, (fs*T)))', 1, N*fs*T);
ASK_res = carrier.*s;

% 计算信号的傅里叶变换
Y = fft(ASK_res);
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/(N*fs*T)); 
P1 = P2(1:N*fs*T/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% 定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，频率和幅值与理论相同
f = fs*(0:(N*fs*T/2))/(N*fs*T);

s_awgn = awgn(ASK_res, -15);
Rp=5;                   % 通带纹波系数
Rs=6;                   % 阻带衰减系数
Wp=14000;               % 通带截止频率
Fs=12000;               % 阻带起始频率

Wp1=[12000 24000];               
Fs1=[16000 20000];
Wp1=Wp1/(fs/2);
Fs1=Fs1/(fs/2);
[n1,Wn1]=buttord(Wp1,Fs1,Rp,Rs);
[b1,a1]=butter(n1,Wn1);
de_car1=filter(b1,a1,s_awgn);

m = de_car1.*carrier

figure
subplot(3,1,1);plot(t, de_car1);axis([0,10*T,-1.5,1.5]);
subplot(3,1,2);plot(t, m);axis([0,10*T,-1.5,1.5]);

Wp=Wp/(fs/2);    % 通带截止频率
Fs=Fs/(fs/2);    % 阻带起始频率
 
% 计算滤波器的最低阶数n和3dB截止频率Wn
[n,Wn]=buttord(Wp,Fs,Rp,Rs);
% 计算系统函数H(z)的分子、分母多项式系数
[b,a]=butter(n,Wn);
y=filter(b,a,m);      % filter函数实现IIR滤波器的直接形式

subplot(3,1,3);plot(t, y);axis([0,10*T,-1.5,1.5]);

for i = 1:N
    if y(L*i-L/2)>=0
        decode_2ASK(i)=1;
    else
        decode_2ASK(i)=0;
    end
end
source;
decode_2ASK;
de_2ASK = [zeros(1, L/2), reshape(repmat(decode_2ASK', 1, L)', 1, N*L)];
k = xor(source, decode_2ASK);
error = length(find(k(1,:)==1))/N;

disp('误码率：')
disp(error)

figure;
subplot(3, 1, 1); plot(t, s); axis([0, 10*T, -1.5, 1.5])
subplot(3, 1, 2); plot(t, ASK_res); axis([0,10*T,-1.5,1.5])
subplot(3, 1, 3); plot(t, s_awgn); axis([0,10*T,-3,3])
% subplot(3,1,3); plot(f,P1); xlabel('f (Hz)'); title('X(t)的单边幅值谱')

figure;
subplot(3,1,1);plot(t, s_awgn);axis([0,10*T,-5,7]);
% subplot(4,1,2);plot(t, x);axis([0,10*T,-2,5]);
subplot(3,1,2);plot(t, y);axis([0,10*T,-1.5,1.5]); 
subplot(3,1,3);plot(t_decode, de_2ASK);axis([0,10*T,-0.5,1.5]); 

figure;
ax.TitleHorizontalAlignment = 'center';
subplot(4, 1, 1); plot(t, ASK_res); axis([0,10*T,-1.5,1.5]);
subplot(4, 1, 2);plot(t, y, 'r');axis([0,10*T,-1.5,1.5]); 
subplot(4, 1, 3); plot(t, s); axis([0, 10*T, -1.5, 1.5]);
subplot(4, 1, 4);plot(t_decode, de_2ASK, 'r');axis([0,10*T,-0.5,1.5]); 
sgtitle('-15dB')

$2 D SP K$

2DPSK
clc
clear
N = 10000                          % 码元数
fs = 900000                        % 采样频率

error_all = zeros(1, 3000)

source = randi([0, 1], 1, N);   % 信源
rel_code = zeros(1, N+1) % 相对码初始化
for q = 2:(N+1)
    rel_code(q) = mod((source(q-1)+rel_code(q-1)), 2);
end
source
rel_code
rel_code(rel_code==0) = -1


T = 1/9000;                        % 单个码元持续时间
rate = 2;                          % 载波频率与码元频率之比
fc = rate*(1/T);                   % 载波频率
t = (0:(N+1)*fs*T-1)/fs;           % 时间刻度

t1 = (0:fs*T-1)/fs; % 绝对码时间刻度
t2 = (fs*T:(N+1)*fs*T-1)/fs;


L = fs*T;
t_decode = (0:(N*fs*T+L/2-1))/fs;
decode_2ASK_rel = zeros(1, N+1);
decode_2ASK = zeros(1, N);

carrier = sin(2*pi*fc*t);       % 载波
Ori = reshape((repmat([0, source]', 1, (fs*T)))', 1, (N+1)*fs*T);
% s = reshape((repmat(AMI', 1, (fs*T)))', 1, N*fs*T);
s = reshape((repmat(rel_code', 1, (fs*T)))', 1, (N+1)*fs*T);
ASK_res = carrier.*s;

s_awgn = awgn(ASK_res, 15);
x = abs(s_awgn);    % 加噪，全波整流
Rp=5;               % 通带纹波系数
Rs=6;               % 阻带衰减系数
Wp=14000;           % 通带截止频率
Fs=12000;           % 阻带起始频率

Wp1=[12000 24000];               
Fs1=[16000 20000];
Wp1=Wp1/(fs/2);
Fs1=Fs1/(fs/2);
[n1,Wn1]=buttord(Wp1,Fs1,Rp,Rs);
[b1,a1]=butter(n1,Wn1);
de_car1=filter(b1,a1,s_awgn);

m = de_car1.*carrier

Wp=Wp/(fs/2);    % 通带截止频率
Fs=Fs/(fs/2);    % 阻带起始频率
 
% 计算滤波器的最低阶数n和3dB截止频率Wn
[n,Wn]=buttord(Wp,Fs,Rp,Rs);
% 计算系统函数H(z)的分子、分母多项式系数
[b,a]=butter(n,Wn);
y=filter(b,a,m);      % filter函数实现IIR滤波器的直接形式

figure
%subplot(4,1,1);plot(t, de_car1);axis([0,10*T,-1.5,1.5]);
subplot(3,1,1);
plot(t1, de_car1(1: fs*T), 'b:'); hold on; plot(t2, de_car1(fs*T+1: end)); axis([0,10*T,-1.5,1.5]);
subplot(3,1,2);
plot(t1, m(1: fs*T), 'b:'); hold on; plot(t2, m(fs*T+1: end)); axis([0,10*T,-1.5,1.5]);
subplot(3,1,3);
plot(t1, y(1: fs*T), 'b:'); hold on; plot(t2, y(fs*T+1: end)); axis([0,10*T,-1.5,1.5]);


for i = 2:(N+1)
    if y(L*i-L/2)>=0
        decode_2ASK_rel(i)=1;
    else
        decode_2ASK_rel(i)=0;
    end
end
source;
decode_2ASK_rel;

decode_2ASK = xor(decode_2ASK_rel(1:end-1), decode_2ASK_rel(2:end))

decode_2ASK;

decode_2DPSK_ = [0 decode_2ASK];

de_2ASK = [zeros(1, L/2), reshape(repmat(decode_2ASK_rel', 1, L)', 1, (N+1)*L)];
result = reshape(repmat(decode_2DPSK_', 1, L)', 1, (N+1)*L);
k = xor(source, decode_2ASK);
error = length(find(k(1,:)==1))/N;

disp('误码率：')
disp(error)

figure;
subplot(3, 1, 1); 
plot(t1, Ori(1: fs*T), ':'); hold on; plot(t2, Ori(fs*T+1: end)); axis([0, 10*T, -0.5, 1.5])
subplot(3, 1, 2); 
plot(t1, ASK_res(1: fs*T), ':'); hold on; plot(t2, ASK_res(fs*T+1: end)); axis([0,10*T,-1.5,1.5])

figure;
subplot(4, 1, 1); 
plot(t1, Ori(1: fs*T), '-.'); hold on; plot(t2, Ori(fs*T+1: end), 'b'); axis([0, 10*T, -0.5, 1.5])
subplot(4, 1, 2); 
plot(t1, result(1: fs*T), '-.'); hold on; plot(t2, result(fs*T+1: end), 'r'); axis([0, 10*T, -0.5, 1.5])
subplot(4, 1, 3); 
plot(t1, ASK_res(1: fs*T), '-.'); hold on; plot(t2, ASK_res(fs*T+1: end), 'b'); axis([0,10*T,-1.5,1.5])
subplot(4,1,4);
plot(t1, y(1: fs*T), '-.'); hold on; plot(t2, y(fs*T+1: end), 'r'); axis([0,10*T,-1.5,1.5]);
sgtitle('15dB')

figure;
subplot(3, 1, 1); 
plot(t1, Ori(1: fs*T), '-.'); hold on; plot(t2, Ori(fs*T+1: end), 'b'); axis([0, 10*T, -0.5, 1.5])
subplot(3, 1, 2); 
plot(t1, s(1: fs*T), '-.'); hold on; plot(t2, s(fs*T+1: end), 'b'); axis([0, 10*T, -1.5, 1.5])
subplot(3, 1, 3); 
plot(t1, ASK_res(1: fs*T), '-.'); hold on; plot(t2, ASK_res(fs*T+1: end), 'b'); axis([0,10*T,-1.5,1.5])

误码率曲线

2PSK
clc
clear
N = 500                       % 码元数
fs = 900000                     % 采样频率
error_all_2PSK = zeros(1, 3000)
error_all_2DPSK = zeros(1, 3000)

laji_2PSK = zeros(1, 3000)
laji_2DPSK = zeros(1, 3000)

T = 1/9000;                     % 单个码元持续时间
rate = 2;                       % 载波频率与码元频率之比
fc = rate*(1/T);                % 载波频率
t = (0:N*fs*T-1)/fs;            % 时间刻度
L = fs*T;
t_decode = (0:(N*fs*T+L/2-1))/fs;
decode_2PSK = zeros(1, N);
carrier = sin(2*pi*fc*t);       % 载波

rel_code = zeros(1, N+1); % 相对码初始化

t_2DPSK = (0:(N+1)*fs*T-1)/fs;            % 时间刻度


%***************滤波器参数*****************
Rp=5;               % 通带纹波系数
Rs=6;               % 阻带衰减系数

Wp1=[12000 24000];               
Fs1=[16000 20000];
Wp1=Wp1/(fs/2);                     
Fs1=Fs1/(fs/2);                     
[n1,Wn1]=buttord(Wp1,Fs1,Rp,Rs);    
[b1,a1]=butter(n1,Wn1);
% de_car1=filter(b1,a1,s_awgn);   % 2PSK

% m = de_car1.*carrier;

Wp=14000;           % 通带截止频率
Fs=12000;           % 阻带起始频率     
Wp=Wp/(fs/2);       % 通带截止频率     
Fs=Fs/(fs/2);       % 阻带起始频率     
% 计算滤波器的最低阶数n和3dB截止频率Wn
[n,Wn]=buttord(Wp,Fs,Rp,Rs);
% 计算系统函数H(z)的分子、分母多项式系数
[b,a]=butter(n,Wn);
% filter函数实现IIR滤波器的直接形式
% y=filter(b,a,m);  % 2PSK    



for z = 1:3000%3000  
    z;
    SNR = -15+0.01*z;
 
source = randi([0, 1], 1, N);   % 信源

AMI = source;
AMI(AMI==0) = -1;


2DPSK初始化

% rel_code = zeros(1, N+1); % 相对码初始化
for q = 2:(N+1)
    rel_code(q) = mod((source(q-1)+rel_code(q-1)), 2);
end
rel_code(rel_code==0) = -1;

decode_2DPSK_rel = zeros(1, N+1);
decode_2DPSK = zeros(1, N);


carrier_2DPSK = sin(2*pi*fc*t_2DPSK);       % 载波
Ori_2DPSK = reshape((repmat([0, source]', 1, (fs*T)))', 1, (N+1)*fs*T);
s_2DPSK = reshape((repmat(rel_code', 1, (fs*T)))', 1, (N+1)*fs*T);
DPSK_res = carrier_2DPSK.*s_2DPSK;

s_awgn_2DPSK = awgn(DPSK_res, SNR);


Ori = reshape((repmat(source', 1, (fs*T)))', 1, N*fs*T);
s = reshape((repmat(AMI', 1, (fs*T)))', 1, N*fs*T);
PSK_res = carrier.*s;

s_awgn = awgn(PSK_res, SNR);    % 加高斯白噪声
x = abs(s_awgn);                % 全波整流

de_car1=filter(b1,a1,s_awgn);   % 2PSK
m = de_car1.*carrier;
y=filter(b,a,m);  % 2PSK  

2DPSK

de_car1_2DPSK=filter(b1,a1,s_awgn_2DPSK);
m_2DPSK = de_car1_2DPSK.*carrier_2DPSK;

y_2DPSK=filter(b,a,m_2DPSK);      % filter函数实现IIR滤波器的直接形式


for i_2DPSK = 2:(N+1)
    if y_2DPSK(L*i_2DPSK-L/2)>=0
        decode_2DPSK_rel(i_2DPSK)=1;
    else
        decode_2DPSK_rel(i_2DPSK)=0;
    end
end
% source;
decode_2DPSK_rel;
for p = 2:(N+1)
    decode_2DPSK(p-1) = mod((decode_2DPSK_rel(p)+decode_2DPSK_rel(p-1)), 2);
end
% decode_2DPSK;

de_2DPSK = [zeros(1, L/2), reshape(repmat(decode_2DPSK_rel', 1, L)', 1, (N+1)*L)];
k_2DPSK = xor(source, decode_2DPSK);
error_2DPSK = length(find(k_2DPSK(1,:)==1))/N;
laji_2DPSK(z) = length(find(k_2DPSK(1,:)==1));

for i = 1:N
    if y(L*i-L/2)>=0
        decode_2PSK(i)=1;
    else
        decode_2PSK(i)=0;
    end
end
source;
decode_2PSK;
de_2ASK = [zeros(1, L/2), reshape(repmat(decode_2PSK', 1, L)', 1, N*L)];
k = xor(source, decode_2PSK);
error = length(find(k(1,:)==1))/N;

error_all_2PSK(z) = error;
error_all_2DPSK(z) = error_2DPSK;
laji_2PSK(z) = length(find(k(1,:)==1));
end

figure;
plot((-15+0.01):0.01:15, error_all_2PSK, 'b-'); hold on;
plot((-15+0.01):0.01:15, error_all_2DPSK, 'r-'); hold on;
xlabel('信噪比/dB'); ylabel('误码率');
rate = sum(laji_2PSK)/sum(laji_2DPSK)

figure;
semilogy((-15+0.01):0.01:15, error_all_2PSK); hold on;
semilogy((-15+0.01):0.01:15, error_all_2DPSK); hold on;
legend({'\it2PSK', '\it2DPSK'}, 'Location', 'northeast');
xlabel('信噪比/dB'); ylabel('误码率');