物理层网络编码的相位同步matlab仿真一、理论基础二、核心程序三、测试结果

122 阅读 0 评论 81 点赞

我是靠谱客的博主和谐高山，这篇文章主要介绍物理层网络编码的相位同步matlab仿真一、理论基础二、核心程序三、测试结果，现在分享给大家，希望可以做个参考。

欢迎订阅《FPGA学习入门100例教程》、《MATLAB学习入门100例教程》

一、理论基础

二、核心程序

三、测试结果

一、理论基础

网络编码技术通过在网络的中间节点处对数据流进行编码来提高无线网络的吞吐量.基于网络编码的理论,物理层网络编码(Physical-Layer Network Coding,PNC)技术利用了电磁波的叠加性,允许来自源/目的节点的无线信号在中继节点处叠加,进而"自然地"完成网络编码;PNC进一步提升了网络吞吐量. 适用于PNC的相位同步方案设计,同步误差影响下M进制正交幅度调制(M-ary Quadrature Amplitude Modulation,M-QAM)PNC的平均误符号率分析以及在二者基础上的同步方案与数据传输的相互作用研究.

目前PNC的实现依然处于研究阶段,相位同步问题是实现PNC的关键因素及主要障碍.本文首先提出了适用于PNC的相位同步方案,随后分别研究了同步方案采用基于最大似然估计(Maximum Likelihood Estimation,MLE)的频率相位估计算法和基于锁相环(Phase-Locked Loop,PLL)的频率相位跟踪的同步精度.

最大似然估计(Maximum Likelihood Estimation)是一种可以生成拟合数据的任何分布的参数的最可能估计的技术。它是一种解决建模和统计中常见问题的方法——将概率分布拟合到数据集。

例如，假设数据来自泊松(λ)分布，在数据分析时需要知道λ参数来理解数据。这时就可以通过计算MLE找到给定数据的最有可能的λ，并将其用作对参数的良好估计。

MLE是用于拟合或估计数据集概率分布的频率法。这是因为MLE从不计算假设的概率，而贝叶斯解会同时使用数据和假设的概率。MLE假设在计算方法之前，所有的解决方案(分布的参数)都是等可能的，而贝叶斯方法(MAP)不是这样，它使用了关于分布参数的先验信息。

MLE之所以有效，是因为它将寻找数据分布的参数视为一个优化问题。通过最大化似然函数，找到了最可能的解。

二、核心程序

复制代码

clc;
clear;
close all;
warning off;
addpath 'func'

%数据长度
Len    = 504;
%网络数据包长度
Npkt   = 1000;  
%网络译码迭代次数
Niter  = 5;    
per    = randperm(2*Len);
SNRs   = [0:2:12];
ij     = 0;
    
for i = SNRs
    i
    sigma = 1/sqrt(0.5*10^(i/10));
    Error = 0;
    ij    = ij+1;
    for k = 1:Npkt
        %QPSK
        Id     = randint(1,Len);%I路
        Qd     = randint(1,Len);%Q路
        %网络编码
        Id_enc = func_enc(Id,per);
        Qd_enc = func_enc(Qd,per);
        Id_enc2= 1-2*Id_enc;
        Qd_enc2= 1-2*Qd_enc;
        Rec_ref= Id_enc2+Qd_enc2;
        
        theta1 = pi/6;
        Id_enc2= Id_enc2.*exp(-sqrt(-1)*theta1)+sigma*randn(1,2*Len);
        Qd_enc2= Qd_enc2.*exp(-sqrt(-1)*theta1)+sigma*randn(1,2*Len);
        %通过干扰
        Rec0   = Id_enc2  + Qd_enc2; 
        %接收
        %MLE进行相位估计
        if k <= 10
           thest0  = atan(sum(imag(conj(Rec_ref).*Rec0))/sum(real(conj(Rec_ref).*Rec0))); 
           thest   = thest0;
        else
           %相位跟踪
           Rec   = Rec0.*exp(sqrt(-1)*thest);
           err   = real(mean(Rec0-Rec));
           C1    = 0.00015;
           thest = thest+C1*err;
        end
 
        Rec    = Rec0.*exp(-sqrt(-1)*thest);
        h      = func_H(Id,per);
        Y      = func_dec(Rec,h,1,1,sigma,Niter);
        Error  = Error + sum(abs(mod(Id+Qd,2) - Y));
    end
    ber(ij) = Error/Npkt/Len;
end

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
clc;
clear;
close all;
warning off;
addpath 'func'

%数据长度
Len    = 504;
%网络数据包长度
Npkt   = 1000;  
%网络译码迭代次数
Niter  = 5;    
per    = randperm(2*Len);
SNRs   = [0:2:12];
ij     = 0;
    
for i = SNRs
    i
    sigma = 1/sqrt(0.5*10^(i/10));
    Error = 0;
    ij    = ij+1;
    for k = 1:Npkt
        %QPSK
        Id     = randint(1,Len);%I路
        Qd     = randint(1,Len);%Q路
        %网络编码
        Id_enc = func_enc(Id,per);
        Qd_enc = func_enc(Qd,per);
        Id_enc2= 1-2*Id_enc;
        Qd_enc2= 1-2*Qd_enc;
        Rec_ref= Id_enc2+Qd_enc2;
        
        theta1 = pi/6;
        Id_enc2= Id_enc2.*exp(-sqrt(-1)*theta1)+sigma*randn(1,2*Len);
        Qd_enc2= Qd_enc2.*exp(-sqrt(-1)*theta1)+sigma*randn(1,2*Len);
        %通过干扰
        Rec0   = Id_enc2  + Qd_enc2; 
        %接收
        %MLE进行相位估计
        if k <= 10
           thest0  = atan(sum(imag(conj(Rec_ref).*Rec0))/sum(real(conj(Rec_ref).*Rec0))); 
           thest   = thest0;
        else
           %相位跟踪
           Rec   = Rec0.*exp(sqrt(-1)*thest);
           err   = real(mean(Rec0-Rec));
           C1    = 0.00015;
           thest = thest+C1*err;
        end
 
        Rec    = Rec0.*exp(-sqrt(-1)*thest);
        h      = func_H(Id,per);
        Y      = func_dec(Rec,h,1,1,sigma,Niter);
        Error  = Error + sum(abs(mod(Id+Qd,2) - Y));
    end
    ber(ij) = Error/Npkt/Len;
end

复制代码

function Dout = func_dec(Din,h0,h1,h2,sigma,Iter)

[M,N] = size(h0);
P0    = exp(-abs(Din+h1+h2).^2/2/sigma^2) + exp(-abs(Din-h1-h2).^2/2/sigma^2); 
P1    = exp(-abs(Din+h1-h2).^2/2/sigma^2) + exp(-abs(Din-h1+h2).^2/2/sigma^2); 
LP    = log(P0./P1);
Lrji  = zeros(M, N);
Lqij  = repmat(LP, M, 1);
%译码迭代
for n = 1:Iter
   for i = 1:M
      col = find(h0(i, :));
      for k = 1:length(col)
         prodOfdrji = 1;
         for l = 1:length(col)
            if l~= k
               prodOfdrji = prodOfdrji*tanh(0.5*Lqij(i,col(l)));
            end
         end % for l
         Lrji(i, col(k)) = 2*atanh(prodOfdrji);
      end % for k
   end % for i
   for j = 1:N
      row = find(h0(:, j));
      for k = 1:length(row)
         dLrij = 0;
         for l = 1:length(row)
            if l~= k
               dLrij = dLrij + Lrji(row(l), j);
            end
         end % for l
         Lqij(row(k), j) = LP(j) + dLrij; 
      end % for k
      LQi = LP(j) + sum(Lrji(row, j));      
      if LQi > 0
         vHat(j) = 0;
      else
         vHat(j) = 1;
      end
   end % for j
   if mod(h0*vHat',2)==0
      break;
   end
end % for n
for r = 1:0.5*N
    Dout(r) = vHat(r);
end

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
function Dout = func_dec(Din,h0,h1,h2,sigma,Iter)

[M,N] = size(h0);
P0    = exp(-abs(Din+h1+h2).^2/2/sigma^2) + exp(-abs(Din-h1-h2).^2/2/sigma^2); 
P1    = exp(-abs(Din+h1-h2).^2/2/sigma^2) + exp(-abs(Din-h1+h2).^2/2/sigma^2); 
LP    = log(P0./P1);
Lrji  = zeros(M, N);
Lqij  = repmat(LP, M, 1);
%译码迭代
for n = 1:Iter
   for i = 1:M
      col = find(h0(i, :));
      for k = 1:length(col)
         prodOfdrji = 1;
         for l = 1:length(col)
            if l~= k
               prodOfdrji = prodOfdrji*tanh(0.5*Lqij(i,col(l)));
            end
         end % for l
         Lrji(i, col(k)) = 2*atanh(prodOfdrji);
      end % for k
   end % for i
   for j = 1:N
      row = find(h0(:, j));
      for k = 1:length(row)
         dLrij = 0;
         for l = 1:length(row)
            if l~= k
               dLrij = dLrij + Lrji(row(l), j);
            end
         end % for l
         Lqij(row(k), j) = LP(j) + dLrij; 
      end % for k
      LQi = LP(j) + sum(Lrji(row, j));      
      if LQi > 0
         vHat(j) = 0;
      else
         vHat(j) = 1;
      end
   end % for j
   if mod(h0*vHat',2)==0
      break;
   end
end % for n
for r = 1:0.5*N
    Dout(r) = vHat(r);
end