【通信系统仿真系列】基于延迟相乘鉴相的2QPSK调相通信系统仿真

1 前言

祝大家五一劳动节快乐！因为新冠肺炎疫情，我过了一个比暑假还长的寒假。是时候补一下当年数字通信原理的坑了。一直弄不清楚为什么延迟相乘可以鉴相，在实验的过程中出现了一些问题，前前后后想了一个星期才解决了。然后我整理了一下仿真原理和实现代码，发布在此。
仿真基于Matlab，我使用的版本是R2016a。

2 原理

2.1 码元生成

用户码元采用双极性码，以-1和1表示，通过随机数方法产生。实验时使用了一万个码元。

2.2 差分处理

2.2.1 不使用差分码

在第一版程序的时候，因为缺少了这个步骤，查找原因花费了我两个星期，给个图大家体会一下。以下实验结果是错误的，正确的请到第三节：实验结果观看。

这个问题困扰了我两个星期：为什么-2dB到2dB之间那一段会出现这样的震荡？一开始以为是实验用的码元数量不够，增加码元数量，问题依旧；我接着以为是实验次数不够，但多次实验后，发现-2dB到2dB这段的震荡还是随机的！
昨天打开代码回顾，忽然顿悟，会出现以上震荡现象是因为鉴相判决时使用了差分判决，每一个新的码元均由前一个已经判决出的码元和相乘结果共同决定，当其中某一位出现判决错误，那么后面的判决结果就会全部错误，正所谓一步错步步错！

2.2.2 使用差分码

为了避免以上差分判决的弊端，在传输前首先对用户码元做差分处理，差分对比图如下图所示。图中为了方便区分，绘制时对差分后的波形幅度乘了2，实际仿真时，差分处理不会改变幅度。

差分码生成规则：

1. 以-1开始；
2. 如果用户码元为1，那么差分输出为上一位输出跳变；
3. 如果用户码元为-1，差分输出与上一位一致。

差分处理示例：输入[-1 -1 1 -1 -1 1 1]，输出[-1 -1 1 1 1 -1 1]。
进过差分处理后，鉴相判决时可以不使用差分判决，判决器的输出仅与当前时刻有关，不再与过去的上一位码元有关，从而有效避免“一步错步步错”这种现象。

2.3 载波

仿真使用1kHz的正弦波作为载波，采样率为25kHz，载波采样点数为25个，采一个周期。载波采样结果如下图所示。

2.4 调制

调制的被调信号是上述差分码，经过相位调制的结果如下图所示。正弦波的一个周期代表一个码元，当差分码为-1时，调制结果与载波反相；为1时则与载波同相。

2.5 经过高斯信道

为了仿真出误码率与信噪比的关系曲线，模拟信号通过高斯信道，在调制结束后，需要给信号加入高斯白噪声，上述2.4调制结果加入噪声后结果如下图所示，图中信号的信噪比为2dB。

在仿真中，加噪声后的信号作为接收端接收到的信号。

2.6 接收滤波

在鉴相前，需要进行带通滤波，移除带外无用信号。

2.6.1 带通滤波器

仿真使用的带通滤波器是第一类FIR滤波器，采用凯赛尔窗函数法设计生成，它的分析如下图所示。

该滤波器的通带为500 ~ 1500Hz，过渡带为1500Hz ~ 2000Hz与0Hz ~ 500Hz，可见其通带具有线性相位的特性。

2.6.2 滤波结果

滤波前后的频谱对比如下图所示。

可见信号的有效部分在滤波后得以保留，而无用部分被过滤掉了，经过滤波，有效提高信号信噪比。
滤波前后的时域图如下图所示。

2.7 延迟&相乘

对上面滤波后的结果进行延迟一拍，因为25个点对应一个载波，所以实际延迟（整体左移）了25个点，最右侧补25个0。把延迟后的信号与延迟前的型号信号相乘，对比图如下所示。

2.8 判决滤波

因为信号相乘后，主要信息都集中在低频部分，特别是直流部分，此时信号的高频部分对于判决是无益的，需要滤除。

2.8.1 低通滤波器

仿真的低通滤波器使用了第一类FIR滤波器，该滤波器采用凯赛尔窗函数法设计，它的分析图如下图所示。

本低通滤波器通带为0Hz ~ 1500Hz，过渡带为1500Hz ~ 2000Hz，其余是阻带，由图中可以看出，在通带和过渡带部分具有线性相位特性。

2.8.2 滤波结果

滤波前后的时域图如下图所示。

为了更好地理解，给出滤波前后频谱对比图。

通过滤波，可以提高码元判决的正确率。

2.9 码元判决

下图为判决结果与用户码元的对比图，可见两个图一模一样。判决正常，仿真模型可以正常工作。

3 实验结果

实验得到的误码率曲线如下图所示。这才是正常的曲线，随着信噪比的提高，误码率呈现指数型下降，并最终达到0误码率。

4 仿真代码

4.1 说明

代码一共有3条主线：

1. 使用差分码的仿真代码(main.m)，结果与理论一致
2. 不使用差分码的仿真代码(undiff.m)，结果与理论不一致
3. 使用了差分码的仿真原理展示代码(theoryDisplay.m)

这些代码共用了多个模块，请注意区分。

4.2 完整代码下载

Github：https://github.com/highskyno1/2QPSK-simulation
链接：https://www.highskyno1.cn/phaseDetector.zip
度盘资源老被吞，放自己服务器上了

4.3 主函数

4.3.1 使用差分码的仿真代码

main.m

%{
    本代码用于延迟差分鉴相通信系统的仿真
%}
%仿真码元数量
codeSize = 1e5;
%载波频率
carrier_freq = 1e3;
%载波采样率
SampleRate = 25*1e3;
%载波采样点数
SamplePoint = 25;
%用户码元
source = bipolarGen(codeSize);
%码元差分处理
source_diff = diffSource(source); 
%载波
carrier = carrierGen(carrier_freq,SampleRate,SamplePoint);
%调制
modulate = modulation(source_diff,carrier);

%接收滤波器
fcuts = [0 500 1500 2000]; %定义通带和阻带的频率
mags = [0 1 0];             %定义通带和阻带
devs = [0.1 0.1 0.1];      %定义通带或阻带的纹波系数
send_filter = getFilter(fcuts,mags,devs,SampleRate);

%鉴相滤波器
fcuts = [1500 2000]; %定义通带和阻带的频率
mags = [1 0];             %定义通带和阻带
devs = [0.1 0.1];      %定义通带或阻带的纹波系数
phase_filter = getFilter(fcuts,mags,devs,SampleRate);

snr_start = -20;
snr_end = 20;
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_index-1) * snr_div;
    %加噪
    send = awgn(modulate,snr,powerCnt(modulate)/length(modulate)*SamplePoint);
    %接收端
    %接收滤波
    len = length(send);
    send = conv(send_filter,send);
    send = arrayCut(send,len);
    %延迟一个码元周期
    delay = [send(SamplePoint+1:end),zeros(1,SamplePoint)];
    %相乘
    phase = delay.*send;
    %鉴相滤波
    len = length(phase);
    phase = conv(phase,phase_filter);
    phase = arrayCut(phase,len);
    %鉴相
    res = detector(phase,SamplePoint);
    rightR = rightRateCnt(res,source);
    errorRate(snr_index) = 1-rightR;
    disp([num2str(snr),'dB:',num2str(rightR)]);
end
foo_x = linspace(snr_start,snr_end,snr_size);
figure
semilogy(foo_x,errorRate);
xlabel('信噪比dB')
ylabel('误码率%')
title('误码率与信噪比的关系')

4.3.2 不使用差分码的仿真代码

以下代码的仿真结果跟理论不一致，仅用于错误演示。
undiff.m

%{
    以下是第一版代码，是错误的
因为没有使用差分鉴相。
错误！
错误！
错误！
%}
codeSize = 1e5;
carrier_freq = 1e3;
SampleRate = 4e4;
SamplePoint = 80;

%用户码元
%规定发送码元必须以-1开始
source = [-1,bipolarGen(codeSize)];
%载波
carrier = carrierGen(carrier_freq,SampleRate,SamplePoint);
%调制
modulate = modulation(source,carrier);

snr_start = -10;
snr_end = 10;
snr_div = 0.1;
snr_size = floor((snr_end-snr_start)/snr_div);
rightRate = zeros(1,snr_size);
parfor snr_index = 1:snr_size
    snr = snr_start + (snr_index-1) * snr_div;
    %加噪
    send = awgn(modulate,snr,powerCnt(modulate)/length(modulate)*SamplePoint);
    %接收端
    %延迟一个码元周期
    delay = [send(SamplePoint+1:end),zeros(1,SamplePoint)];
    %相乘
    phase = delay.*send;
    %鉴相
    len = length(phase)/SamplePoint-1;
    res = [-1,zeros(1,len)];
    for i = 1:len
        foo = sum(phase((i-1)*SamplePoint+1:i*SamplePoint));
        if foo > 0
            res(i+1) = res(i);
        else
            res(i+1) = -res(i);
        end
    end
    rightR = rightRateCnt(res,source);
    rightRate(snr_index) = 1-rightR;
    disp([num2str(snr),'dB:',num2str(rightR)]);
end
plot(linspace(snr_start,snr_end,snr_size),rightRate);

4.3.3 原理展示代码

theoryDisplay.m

%{
    本代码用于仿真原理展示，
    仿真代码放到了main.m
%}
carrier_freq = 1e3;
SampleRate = 25*1e3;
SamplePoint = 25;

close all;
%用户码元
source = [-1 1 1 -1 -1];
figure;
stairs(source);
hold on;
%码元差分处理
source_diff = diffSource(source);
stairs(source_diff.*2);
title('用户码元')
legend('差分前','差分后');
set(gca,'YLim',[-3 3]);%Y轴的数据显示范围
%载波
carrier = carrierGen(carrier_freq,SampleRate,SamplePoint);
figure;
plot(carrier);
title('载波')
%调制
modulate = modulation(source_diff,carrier);
figure;
plot(modulate);
title('调制结果')

%接收滤波器
fcuts = [0 500 1500 2000]; %定义通带和阻带的频率
mags = [0 1 0];             %定义通带和阻带
devs = [0.1 0.1 0.1];      %定义通带或阻带的纹波系数
receive_filter = getFilter(fcuts,mags,devs,SampleRate);
figure;
freqz(modulate);
title('接收滤波器分析')

%鉴相滤波器
fcuts = [1500 2000]; %定义通带和阻带的频率
mags = [1 0];             %定义通带和阻带
devs = [0.1 0.1];      %定义通带或阻带的纹波系数
phase_filter = getFilter(fcuts,mags,devs,SampleRate);
figure;
freqz(phase_filter);
title('鉴相滤波器分析')

snr = 2;
%加噪
receive = awgn(modulate,snr,powerCnt(modulate)/length(modulate)*SamplePoint);
figure;
plot(receive);
title('加噪结果')

%接收端

%接收滤波
figure;
len = length(receive);
plot(abs(fft(receive)))
title('滤波频谱图');
hold on;
receive2 = conv(receive_filter,receive);
receive2 = arrayCut(receive2,len);
plot(abs(fft(receive2)));
legend('滤波前','滤波后');
figure;
plot(receive);
hold on;
plot(receive2);
title('接收滤波结果')
legend('滤波前','滤波后');

%延迟一个码元周期
delay = [receive2(SamplePoint+1:end),zeros(1,SamplePoint)];
%相乘
phase = delay.*receive2;
figure;
subplot(3,1,1);
plot(receive2);
title('相乘前')
subplot(3,1,2);
plot(delay);
title('相乘前延迟一拍')
subplot(3,1,3);
plot(phase);
title('相乘后')

%鉴相滤波
len = length(phase);
phase2 = conv(phase,phase_filter);
phase2 = arrayCut(phase2,len);
figure;
plot(phase);
hold on;
plot(phase2);
title('鉴相滤波');
legend('滤波前','滤波后');
figure;
plot(abs(fft(phase)));
hold on;
plot(abs(fft(phase2)));
title('鉴相滤波频谱图');
legend('滤波前','滤波后');

%鉴相
res = detector(phase2,SamplePoint);
figure;
subplot(2,1,1);
stairs(res);
title('鉴相结果')
set(gca,'YLim',[-2 2]);%Y轴的数据显示范围
subplot(2,1,2);
stairs(source);
title('用户码元')
set(gca,'YLim',[-2 2]);%Y轴的数据显示范围

4.4 模块代码

4.4.1 双极性随机码产生模块

bipolarGen.m

function [code] = bipolarGen(size)
%bipolarGen 产生指定长度的二进制双极性码元数组
%   码元的产生采用随机函数
    foo = rand(1,size);
    for i = 1:length(foo)
        if foo(i) > 0.5
            foo(i) = 1;
        else
            foo(i) = -1;
        end
    end
    code = foo;
end

4.4.2 载波生成模块

carrierGen.m

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

4.4.3 差分处理模块

diffSource.m

function res = diffSource(source)
%{
本函数用于把绝对码转换成差分码,
差分方法：输出结果第一个必须为-1，如果码元是-1，那么输出结果与上一个相同；
    如果码元是1，那么输出结果为对上一个输出取反。
示例：输入[-1 -1 1 -1 -1 1 1]
示例：输出[-1 -1 1 1 1 -1 1]
%}
    res = zeros(1,length(source));
    res(1) = -1;
    for i = 2:length(source);
        if source(i) < 0
            res(i) = res(i-1);
        else
            res(i) = -res(i-1);
        end
    end
end

4.4.4 载波生成模块

modulation.m

function [res] = modulation(code,carrier)
    %调相调制器
    %code:需要调相的码元
    %carrier:载波
    res = zeros(1,length(code) * length(carrier));
    for i = 1:length(code)
        if code(i) > 0
            res((i-1)*length(carrier)+1:i*length(carrier)) = carrier;
        else
            res((i-1)*length(carrier)+1:i*length(carrier)) = -carrier;
        end
    end
end

4.4.5 FIR滤波器生成模块

getFilter.m

function result = getFilter(fcuts,mags,devs,freq)
%getFilter 获取滤波器，滤波器由凯赛尔窗函数法生成
%   fcuts:定义通带和阻带的频率
%   mags:定义通带和阻带
%   devs:定义通带或阻带的纹波系数
%   freq:采样率
%{
    demo
    fcuts = [200 900 1100 2000]; %定义通带和阻带的频率
    mags = [1 0 1];             %定义通带和阻带
    devs = [0.1 0.01 0.1];      %定义通带或阻带的纹波系数
%}
    [n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,freq);  %计算出凯塞窗N，beta的值
    result = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');   %生成滤波器
end

4.4.6 居中截取模块

arrayCut.m

function res = arrayCut(array,size)
%arrayCut 根据需要的长度，居中截取数组
%   因为在FIR滤波之后，结果会增长，对于解调来说，
%   滤波卷积结果两端的部分是多余的，应该去除。
    start = floor((length(array)-size)/2+1);
    res = array(start:start+size-1);
end

4.4.7 能量计算模块

powerCnt.m

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

4.4.8 正确率统计模块

rightRateCnt.m

function res = rightRateCnt(arg0,arg1)
%rightRateCnt 统计两个数组的相识程度
    rightCnt = 0;
    for i = 1:length(arg0)
        if(arg0(i) == arg1(i))
            rightCnt = rightCnt + 1;
        end
    end
    res = rightCnt / length(arg0);
end

5 后语

前几个星期出了考研国家线，忽然发现自己专业课刚好压线，获得了一张复试的门票！当然以我这样的成绩去竞争，能考上的概率依然很低，但有了一次复试的机会，可以给以后重新考研增加经验。我觉得我们这一年的学生总能遇上波折，高考改用全国卷的我们是第一届；到了大学，旧的教学任务大纲我们是最后一届；新冠肺炎来了，我们刚好是大四毕业！
看到这了，不点个赞再走吗？