2DPSK调制和差分解调的Matlab实现

------------12.25更新--------------
我果然是半桶水，题目稍微改一改就做不出来了。实验考试考了2ASK调制解调系统的，当信噪比很大的时候，误码率死活都降不下来，我还一直以为是调制解调哪里出了问题，最后没想到是抽烟判决点

remod_bits(val_decision > 0.5) = 1;     

2ASK的判决门限是0.5啊～～～
---------------------分割线----------------
一早上的时间都用来写2DPSK调制和差分解调，ctrl一时爽，debug一早上。

我为啥会调一早上呢，主要是我以为2DPSK的调制解调应当和2PSK的是一样的，于是就把之前的2PSK的代码直接ctrl过来，只是中间增加了一个取差分码的过程，改掉一些错误后发现最后解调出来的误码率高达0.5，然后我就走上了调代码的不归路，调完以后觉得写代码还是要一步一步写，当代码的结果运行出想要的结果是再往下写，效率才会高；写完一堆代码，发现结果不对，再回去调会不可控因素就变多了，效率反而更低。

记录一下踩过的坑，希望以后不要再踩了

具体问题：

绘制2DPSK调制与差分解调的流程及对应的波形，仿真绘制误码率曲线，并与理论误码曲线进行比较(SNR=-3dB…8dB)。

实现的流程及其对应的波形
根据流程图去看代码一目了然，但有些细节还是需要再解释下。

一、调制过程

1、码型变换：由绝对码变成相对码
2DPSK的编码规则是当码元为1时，其波形的相位就要发生翻转，对应的相对码就要发生翻转，所以就是当绝对码当前码元为1时，相对码下一时刻要发生翻转

代码实现：相对码的下一码元=绝对码的当前码元+相对码的当前码元 %2.
DiffBits(i) =mod( DiffBits(i - 1) + bits(i-1) , 2 ) ;
注意：此时码元长度增加了1

2、时间序列的产生
t = (0:sample_bit(nsym+1) - 1)’/fs; %时间长度*
注意：因为码元数增加了1，则对应的时间序列也要加长

二、差分解调
1、产生延迟波形
delay_t 的第二个码元往后与r_t的第一个码元开始一一对应
delay_t = zeros(length(r_t),1);
delay_t(sample_bit+1:length(r_t)) = r_t(1:length(r_t)-sample_bit);

2、滤波
使用butterworth filter：
截止频率点：码元传输速率Rb
[b, a] = butter(4,Rb/(fs/2));

3、 抽样
1）系统延迟：从下图可以观察出，包络检波后系统是发生了延时的，通过对比两张图的对应点，我取了延迟 group_delay = 30;

2）最佳抽样点：中间时刻+系统延时,round表示取最接近的整数
round(1+sample_bit+group_delay+sample_bit/2):sample_bit:length(r_t);

3）且因为第一个码元没有意义，所以第一个采样点在第二个采样点处
round(1+sample_bit+group_delay+sample_bit/2)

4、判决
remod_bits = ones(length(val_decision)+1,1);
remod_bits(val_decision > 0) = 0;

1）因为采样时加了大于sample_bit/2的系统延时，所以会导致码元数小于原码元数，所以长度要加1:length(val_decision)+1

2）因为后面判决的是对应的0码元，一开始产生的则就应该为1的向量；如果后面判决为对应的1码元，一开始产生的则就应该为0的向量；
代码的另一种形式：
remod_bits = ze ros(length(val_decision)+1,1);
remod_bits(val_decision < 0) = 1;

完整的代码：

clear all
clc

%%
%第二题
%2PSK 调制与解调模型
Rb = 1e3;
fc = 5e3;
fs = 50e3;
nsym = 8000;      %信息长度
Ac = 1;                %载波幅度

%随机产生绝对码
bits = randi(2,nsym,1);
bits = bits -1;

%码型变换-相对码
DiffBits = zeros(nsym+1 ,1);        %码元长度增加了
for i = 2:nsym+1
    DiffBits(i) =mod( DiffBits(i - 1) + bits(i-1) , 2 ) ;
end

%码型变换至双极性不归零码
 Diffbb = 2*DiffBits- 1;                  


sample_bit = fs/Rb;                     %每个符号的采样点数
t = (0:sample_bit*(nsym+1) - 1)'/fs;  %时间长度


 %产生不归零码的时域波形
 %bb是消息的长度，要处理的应该是经过采样率采样后的序列
 %repmat矩阵的复制；reshape:改变矩阵排布，读取方式是一列一列往后读取
Diffbb_t = reshape(repmat(Diffbb',sample_bit,1) ,sample_bit*(nsym+1),1);                      

%调制
c_t = Ac * cos(2* pi * fc * t);         
s_t = Diffbb_t .* c_t;                            

figure (1), plot(s_t),title('s_t');axis([0,350,-1,1]);

%%
%绘制发送端差分编码前后的波形
bb = 2*(bits ) - 1;
bb_t = reshape(repmat(bb',sample_bit,1) ,sample_bit*nsym,1);  
figure,subplot(211),plot(t(1:length(t)-sample_bit),bb_t);title('差分编码前的波形');
subplot(212),plot(t,Diffbb_t);title('差分编码后的波形');


%%
SNR=(-3:1:8)';
B = 2*Rb;           %信道带宽
N = length(s_t);
SimuRat = zeros(length(SNR),1);
TheorRat = zeros(length(SNR),1);

for i = 1:length(SNR)
    snr = SNR(i);
    n_pow = Ac*Ac/2/10^(0.1*snr);       %求解噪声功率
    n_t = bandlimit_noise(N,fs,fc,B,n_pow);
    r_t = n_t + s_t;                                    %信号+信道的噪声
    [remod_t,remod_bits] = BDPSK_remodulate(r_t,Rb,fs,nsym);

   [number,ratio] = symerr(bits,remod_bits);
   SimuRat(i) = ratio;                                      %仿真误码率
   TheorRat(i) = 0.5*exp(-1*10^(0.1*snr));%理论误码率


end
figure,subplot(121),plot( SNR,SimuRat);title('2DPSK仿真误码率曲线');
subplot(122),plot(SNR,TheorRat);title('2DPSK理论误码曲线');

差分解调的函数

function [remod_t,remod_bits] = BDPSK_remodulate(r_t,Rb,fs,nsym)
%差分解调
%nsym 信息长度

    sample_bit = fs/Rb;                     %每个符号的采样点数

    t = (0:sample_bit*(nsym+1) - 1)'/fs;  %时间长度

    %产生延迟波形
    delay_t = zeros(length(r_t),1);
    delay_t(sample_bit+1:length(r_t)) = r_t(1:length(r_t)-sample_bit); 

    %相乘器
    multi_t = delay_t.*r_t;

    %滤波-包络检波
     [b, a] = butter(4,Rb/(fs/2));
    remod_t = filter(b, a, multi_t);   
 
     
    %抽样
     group_delay = 30;%系统延迟
    idx_decision = round(1+sample_bit+group_delay+sample_bit/2):sample_bit:length(r_t); %抽样点下标
                         %最佳抽样点：中间时刻+系统延时,round表示取最接近的整数
                         %第一个抽样点+sample_bit，因为差分码会有一个周期的延迟
    val_decision = remod_t(idx_decision);%取得抽样点

    %判决
    remod_bits = ones(length(val_decision)+1,1);
    remod_bits(val_decision > 0) = 0;                  

end

产生带限噪声的函数

function n_t = bandlimit_noise(N,fs,fc,B,n_pow)
% 产生特定功率的带限高斯白噪声
%  输入： N  - 产生的噪声的长度
%        fs - 仿真采样率
%        fc - 带限噪声中心频率
%        B  - 带限噪声的带宽
%        n_pow - 带限噪声的功率
%  输出： n_t - 带限噪声

    %带通滤波器
    filter_b_L = 4000;
    filter_b = fir1(filter_b_L, [(fc-B/2)/(fs/2), (fc+B/2)/(fs/2)]);

    n_fullband_pow = n_pow/(B/(fs/2));
    noise_t = sqrt(n_fullband_pow)*randn(N+filter_b_L,1);
    n_t = filter(filter_b,1,noise_t);
    n_t = n_t(end-N+1:end);

end