【通信系统仿真系列】基于延迟相乘鉴相的4Q-DPSK差分通信系统仿真

sys_rst_n

已于 2022-11-04 12:37:38 修改

阅读量1.6k

点赞数 13

分类专栏：仿真文章标签： matlab

于 2020-05-06 23:20:22 首次发布

本文链接：https://blog.csdn.net/weixin_41476562/article/details/105948026

版权

仿真专栏收录该内容

10 篇文章 23 订阅

订阅专栏

【通信系统仿真系列】基于延迟相乘鉴相的4Q-DPSK差分通信系统仿真

前言
原理
实验结果
仿真代码
后语

前言

之前做出了2Q-DPSK系统，感兴趣4Q系统的解调判决如何进行。发现网上关于这方面的原理仿真资料太少了，《数字通信原理》课本也只是简单地介绍了可以进行延迟相乘的非相干解调，或者提取同步信号做相干解调，具体怎么做并没有说法。经过五一黄金周这几天的研究，加上上一篇文章《2Q-DPSK通信系统》
的经验，用了三天就仿真出来了，放上来给大家学习参考。
仿真代码基于Matlab编写，笔者仿真使用的是R2016a版。

原理

码元

码元是四进制的，用0、1、2、3表示，基于随机数分区法生成。

载波

载波采用1kHz正弦波，通过sin函数生成，针对0、1、2、3码元，采用不同的初始相位，载波波形图如下所示。
四相位载波
每个载波取正弦波的一个完整周期，初始相位分别为0，π/2，π以及3*π/2，分别对应0、1、2、3四种码元。

差分处理

假如不使用差分传输，在接收端进行的将是差分判决，即判决输出的结果跟上一次的输出结果有关，此时假如某一位判决出错，就会导致以后的判决都出错，所以要避免使用差分判决。
因此，发送前需要进行差分处理，如此，判决输出结果则只会与当前时刻的解调结果有关。
差分前后的对比如下图所示。
差分效果图
上图中差分输入：[1 2 3 0 3 1 0 2]，差分输出：[0 1 2 0 3 3 0 1 1 3]。差分输出比输入多了两个码元，且必须以0和1开始。
差分分为两部分进行，以下分别说明。

一次差分

一次差分针对奇数码元1和3，首先给输入码元开始添加码元1，然后开始往后遍历码元数组，在遇到奇数码元时，如果该码元是1，那么令该码元与上一次出现的奇数码元相同；如果该码元是3，那么令该码元等于上一次出现的奇数码元+2，然后除以4取余（即如果该码元是3且上一次出现的奇数码元是1，则令该码元为3，否则令该码元为1）。对于偶数码元0和2，一次差分不做处理。迭代公式如下所示：
$\begin{cases} 1 &\text{if }n = 1 \\ (diff1(k)+2)\%4 &\text{if }input(n-1) = 3 \\ diff1(k) &\text{if }input(n-1) = 1 \\ input(n-1) &\text{else } \\ \end{cases} \\ 其中,n \subset Z^+,input表示用户四进制码元,diff1表示一次差分输出。\\ k表示diff1不大于n且最接近n的奇数位的索引，即：\\ ∀x＜n且diff1(x)\%2=1,k=max(x)。$
示例：一次差分前：[3 2 0 1 1 2 3]，一次差分后：[1 3 2 0 3 3 2 1]

二次差分

二次差分针对一次差分的结果进行，并规定第一位为0，从第二位开始，每一位是其前一位加上一次差分的该位的结果除以4取余，迭代公式如下：
$\begin{cases} 0 &\text{if }n = 1 \\ (diff1(n-1)+res(n-1))\%4 &\text{if }n > 1 \\ \end{cases} \\ 其中,n \subset Z^+,diff1表示一次差分结果$
示例：
二次差分前：[0 1 2 3 2 1 0]
二次差分后：[0 0 1 3 2 0 1 1]

调制

调制根据差分处理结果与载波进行，差分结果与初相位的对应关系如下表所示：

差分后码元	初相位（弧度制）
0	0
1	pi/2
2	pi
3	3*pi/2
下图是针对差分后码元：[0 1 2 0 3 3 0 1 1 3]的调制结果

加入噪声

为了更好地演示原理，在原理展示代码（theoryDisp.m）中没有这一步，但为了仿真出信噪比与误码率的关系曲线，此步存在于实际的仿真代码（main.m）中。
加入高斯白噪声使用的是由Matlab提供的awgn函数，调用时传入三个参数：

需要加噪的信号
信噪比（dB）
原信号功率（dBW）

延迟相乘解调

接收端接收到的是含有噪声的信号，对这信号延迟一拍，然后与延迟前的进行相乘，结果如下图所示：
延迟相乘结果
解调前与解调后的频谱对比图如下图所示：
解调前后频谱对比图

判决

把解调结果以每采样点数个点（默认25个）进行分组。判决时，对于偶数码元0和2，只需要看正弦信号的偏置，往上偏置（整个波形基本在y=0以上）的为码元0，往下偏置（整体基本在y=0以下）的为码元2；对于奇数码元1和3，它们没有偏置，都出现在y=0附近，此时就要与上一次出现奇数码元做对比（这也是一次差分时需要以1开始的原因），如果本码元周期的信号与上一次同相，那么本码元为3，如果与上一次反相，那么本码元为1。
如下图所示的解调结果：

已知1 ~ 25的部分表示码元1，因为26 ~ 50部分与1 ~ 25部分反相，所以它是1；51 ~ 75部分偏置到了0以下，所以为2；76 ~ 100部分因为与上一次0偏置部分（26 ~ 50）同相，所以是3。以此类推，最终判决出来为：[1 2 3 0 3 1 2 0]。
判决结果与输入码元对比如下：
判决对比
可见，判决结果完全正确。

实验结果

仿真结果如下图所示。
信噪比&误码率关系曲线
可见在误码率小于某个值时趋于平缓，渐进于0.75，大于某个值时，趋于0；在趋于平缓前，呈现指数型下降，这与理论是符合的。

仿真代码

仿真代码有两条主线：

仿真代码(main.m)
仿真原理展示代码(theoryDisplay.m)

可修改的参数

codeSize：仿真码元数量，默认1万个
carrier_freq：载波频率，默认1kHz，根据奈奎斯特定理，为了不出现频谱混叠，必须小于采样率的一半
carrier_freq：采样率，默认25kHz，同样的原因，必须大于两倍的载波频率
SamplePoint：载波采样点数，最好能采到一个完整的正弦波周期
偶数码元判决阈值，请到judge.m文件中修改

下载链接

Github：https://github.com/highskyno1/4Q-DPSK-simulation
度盘：链接：https://pan.baidu.com/s/1-NJrjLofM33YzgQ5A-if_Q?pwd=bnhl
提取码：bnhl

主函数

仿真代码

main.m

%{
    本代码用于对四进制差分移相键控系统（4Q-DPSK）进行仿真，
    绘制其信噪比与误码率的关系曲线
%}
close all;
%仿真码元数量
codeSize = 1e5;
%载波频率
carrier_freq = 1e3;
%载波采样率
SampleRate = 25*1e3;
%载波采样点数
SamplePoint = 25;

%生成四种相位的载波
carrier = [carrierGen(carrier_freq,SampleRate,SamplePoint,0);
    carrierGen(carrier_freq,SampleRate,SamplePoint,pi/2);
    carrierGen(carrier_freq,SampleRate,SamplePoint,pi);
    carrierGen(carrier_freq,SampleRate,SamplePoint,3*pi/2)];

%生成用户码元
%sourceCode = [2 3 3 2 2];
sourceCode = getSourceCode(codeSize);
%生成差分码
diffCode = myDiff(sourceCode);
%调制
modu = modulate(diffCode,carrier,SamplePoint);
send = modu;

snr_start = -30;
snr_end = 10;
snr_div = 0.1;
snr_size = floor((snr_end - snr_start)/snr_div);
errorRate = zeros(1,snr_size);
parfor snr_index = 1:snr_size
    snr = snr_start + snr_div * snr_index;
    receive = awgn(send,snr,powerCnt(send));
    %解调
    delay = receive(SamplePoint+1:end);
    receive2 = receive(1:(length(diffCode)-1)*SamplePoint);
    phase = delay.*receive2;
    %判决
    res = judge(phase,SamplePoint);
    errorRate(snr_index) = 1 - rightRateCnt(res,sourceCode);
end
plot(linspace(snr_start,snr_end,snr_size),errorRate);
title('信噪比&误码率关系曲线');
xlabel('信噪比(dB)');
ylabel('误码率');

原理展示代码

theoryDisp.m

%{
    本代码用户原理展示，不含加入白噪声的代码
    仿真代码放到了main.m
%}
close all;
%仿真码元数量
codeSize = 5;
%载波频率
carrier_freq = 1e3;
%载波采样率
SampleRate = 25*1e3;
%载波采样点数
SamplePoint = 25;

%生成四种相位的载波
carrier = [carrierGen(carrier_freq,SampleRate,SamplePoint,0);
    carrierGen(carrier_freq,SampleRate,SamplePoint,pi/2);
    carrierGen(carrier_freq,SampleRate,SamplePoint,pi);
    carrierGen(carrier_freq,SampleRate,SamplePoint,3*pi/2)];

figure;
plot(carrier(1,:));
hold on;
plot(carrier(2,:));
hold on;
plot(carrier(3,:));
hold on;
plot(carrier(4,:));
hold on;
title('载波');
legend('初相位0','初相位pi/2','初相位pi','初相位3*pi/2');

%生成用户码元
sourceCode = [1 2 3 0 3 1 0 2];

%生成差分码
diffCode = myDiff(sourceCode);
figure;
subplot(2,1,1);
stairs(sourceCode);
title('差分前');
set(gca,'XLim',[0 11]); %X轴的数据显示范围
set(gca,'YLim',[-1 4]); %Y轴的数据显示范围
subplot(2,1,2);
stairs(diffCode);
title('差分后');
set(gca,'XLim',[0 11]); %X轴的数据显示范围
set(gca,'YLim',[-1 4]); %Y轴的数据显示范围

%调制
modu = modulate(diffCode,carrier,SamplePoint);
figure;
plot(modu);
title('调制结果');

%演示不加噪声
receive = modu;

%解调
%延迟一个码元周期
delay = receive(SamplePoint+1:end);
%移除最后一个码元周期
modu2 = receive(1:(length(diffCode)-1)*SamplePoint);
%相乘
phase = delay.*modu2;
figure;
subplot(3,1,1);
plot(modu2);
title('接收信号')
subplot(3,1,2);
plot(delay);
title('延迟结果');
subplot(3,1,3);
plot(phase);
title('相乘结果');
figure;
fft_point = 512;
plot(linspace(0,1,fft_point),abs(fft(modu,fft_point)));
hold on;
plot(linspace(0,1,fft_point),abs(fft(phase,fft_point)));
title('频谱图');
xlabel('数字频率(π/rad)')
legend('解调前','解调后');

%判决
res = judge(phase,SamplePoint);
figure
plot(phase);
figure;
subplot(2,1,1);
stairs(sourceCode);
title('用户码元');
set(gca,'XLim',[0 9]); %X轴的数据显示范围
set(gca,'YLim',[-1 4]); %Y轴的数据显示范围
subplot(2,1,2);
stairs(res);
title('判决结果');
set(gca,'XLim',[0 9]); %X轴的数据显示范围
set(gca,'YLim',[-1 4]); %Y轴的数据显示范围

调用函数

码元生成模块

getSourceCode.m

function res = getSourceCode(size)
%getSourceCode 获取四极性码元
%size:码元数量
    res = rand(1,size);
    for i = 1:length(res)
        foo = res(i);
        if foo > 0.75
            res(i) = 3;
        elseif foo > 0.5
            res(i) = 2;
        elseif foo > 0.25
            res(i) = 1;
        else
            res(i) = 0;
        end
    end
end

载波生成模块

carrierGen.m

function [res] = carrierGen(freq,sampleRate,size,phase)
%carrierGen 产出载波
%   参数：freq:频率 sampleRate:采样率 size:载波长度 phase:初相位
    n=0:size-1;
    t=n/sampleRate;%时间序列
    res=sin(2*pi*freq*t+phase);
end

差分模块

myDiff.m

function res = myDiff(input)
%myDiff 我的差分函数
%input：需要差分处理的四进制信号
    input = [1 input];
    %先对1和3做差分
    %记录上一次的奇数状态
    last_odd = 1;
    for i = 2:length(input)
        if mod(input(i),2) == 1
            if input(i) == 1
                input(i) = last_odd;
            else
                input(i) = mod(last_odd + 2,4);
            end
            last_odd = input(i);
        end
    end
    %对全部码元做2次差分
    res = zeros(1,length(input)+1);
    for i = 1:length(input)
        res(i+1) = mod(res(i) + input(i),4);
    end
end

调制模块

modulate.m

function res = modulate(code,carrier,sampleSize)
%modulate 调制
%code：用户四进制码元
%carrier：四相位载波
%sampleSize：每一码元对应采样点数
    res = zeros(1,length(code)*sampleSize);
    for i = 1:length(code)
        res((i-1)*sampleSize+1:i*sampleSize) = carrier(code(i)+1,:);
    end
end

解调模块

deModulate.m

function res = deModulate(input,SamplePoint)
%deModulate 解调模块
%input：接收到的信号
%解调基于延迟相乘的非相干解调
    %延迟（左移）一拍
    delay = input(SamplePoint+1:end);
    %移除最后一个码元周期
    input = input(1:length(input)-SamplePoint);
    %点对点相乘
    res = delay.*input;
end

判决模块

judge.m

function res = judge(input,SamplePoint)
%judge 码元判决
%input：解调后信号
%SamplePoint：每一个四进制码元对应的采样点数

%记录上一次出现偏置为0的码元周期波形
start_odd = input(1:SamplePoint);
%移除已知结果为1的第一个码元周期
input = input(SamplePoint+1:end);
%初始化res，提高运行速度
res = zeros(1,length(input)/SamplePoint);
for i = 1:length(input)/SamplePoint
    %截取对应的码元周期
    foo = input((i-1)*SamplePoint+1:i*SamplePoint);
    %计算该组码元周期的平均值
    bar = mean(foo);
    if bar > 0.25
        %平均值大于阈值，说明偏置到了0上方
        res(i) = 0;
    elseif bar < -0.25
        %平均值小于阈值，说明偏置到了0下方
        res(i) = 2;
    else
        %没有偏置，说明是奇数码元1或3
        %与上次的偏置为0的码元周期波形做比较
        bar = foo.*start_odd;
        if sum(bar) > 0
            %同相为3
            res(i) = 3;
        else
            %反相为1
            res(i) = 1;
        end
        %记录这一次偏置为0的码元周期波形
        start_odd = foo;
    end
end
end

信号功率计算模块

powerCnt.m

function res = powerCnt(input)
%powerCnt 功率计算函数
%参数 input:需要计算功率的数组
    res = 10*log10(sum(input.*input)/length(input));
end

后语

目前（2020/5/6）已经收到学校的考研复试通知，因为疫情原因要采用远程复试，还要设置双摄像机，一个在面前，另一个在后背45°的位置，另一方面也收到了毕业设计的答辩通知，事情变得有点麻烦。不知道复试结果会怎样，如果有下一篇文章的话，我应该会更新在那篇文章的后语。
看到这了，不点个赞再走吗？

sys_rst_n

关注

13
点赞
踩
15

收藏

觉得还不错? 一键收藏
6
评论
【通信系统仿真系列】基于延迟相乘鉴相的4Q-DPSK差分通信系统仿真

【通信系统仿真系列】基于延迟相乘鉴相的4Q-DPSK差分通信系统仿真前言原理码元载波差分处理一次差分前言之前做出了2Q-DPSK系统，自己感兴趣4Q系统的解调判决如何进行。查资料发现网上关于这方面的原理仿真资料太少了，《数字通信原理》课本也只是简单地介绍了可以进行延迟相乘的非相干解调，或者提取同步信号做相干解调，具体怎么做并没有说法。经过五一黄金周这几天的研究，加上上一篇文章《2Q-DPSK通...
复制链接

扫一扫