OTFS学习记录3——MRC接收机

我将研究最大比率组合（MRC）检测，它提供了与MP检测相似的性能，但复杂度要低得多。所提出的MRC检测方法是由码分多址（CDMA）系统中传统的耙接收机驱动的，因为OTFS可以解释为信息符号在时间和频率上扩展的二维CDMA，因此可以将MRC与OTFS相结合，同时基于该MRC接收机我也会给出相应算法对应的matlab仿真代码。

Background

直接序列CDMA在多路径衰落信道中运行时，rake接收机组合传输符号的延迟分量（或回波），可使用调谐到各自延迟移位的匹配滤波器提取。同样，在OTFS中，可以提取在接收端接收到的不同延迟分支的信道受损信号分量采用MRC等分集组合技术进行相干组合，以提高累积信号的信噪比（SNR）。MRC接收机图如下：

众所周知，MRC是独立AWGN信道的最优组合器。除此之外，如果根据的信号振幅-调谐功率比来选择分支权值，则对于相关通道也是最优的。对于OTFS，每个分支中的噪声加干扰（NPI）是相关的，因为它依赖于信道，这引入了跨分支的干扰。在MRC检测中，我们可以抵消组合选择的分支的估计干扰，迭代提高MRC后信干扰加噪比（SINR）。

下面，我们提出了ZP-OTFS的延迟-多普勒MRC检测和降低复杂度的延迟-时间MRC检测。这些MRC检测方法可以很容易地扩展到其他的OTFS变体。为了便于推导，我们考虑N个ZPs，每个长度为L_ZP = lmax，是一个大小为M’N的扩展OTFS框架的一部分，其中M'= M + lmax。在接收机上，ZP不被丢弃，而是用于MRC检测。

Delay-Doppler domain MRC detection

在每个块中都有一个ZP，我们将传输和接收的符号向量设为

x,y矩阵大小为N*（M'-1）且矩阵中的元素Xm=0当m=M，....M-1时（因为他们是ZP序列的元素，不是信息）

请注意，（6.32）中的符号向量包括ZP。根据ZP的影响，（6.3）中的多普勒输入-输出关系可以修改为

接下来，我们进行了迭代的检测过程，以进行估计Xm（仅需要估计m=0~M-1）

为了便于表示法，我们在下面的步骤中省略了迭代索引

在该方案中，我们不是从（6.35）中的每个方程中单独估计x^m，而是对通道受损分量bl m进行MRC计算，运算如下：

然后是逐符号最大似然检测（MLD），得到由硬估计结果：

其中aj是来自大小为Q的QAM字母表A，其中a=1~Q，n=0~N-1。硬决策函数D（·）由（6.39）中的MLD准则给出。

一旦我们更新了估计x^m[n]，我们增加m并重复相同的过程来估计所有M符号向量xm，使用之前解码的估计量xm（见图6.2）。请注意，DFE操作会导致顺序更新，而只使用之前的迭代估计值会导致并行更新。我们后来通过实验验证了并行更新会导致较慢的收敛速度。算法2提供了详细的延迟-多普勒MRC操作，如下图所示：

在算法2中，初始计算包括生成H的所有条目包括多普勒频移vml，以及i = 1,2，...，P的整数延迟多普勒信道参数（hi，li，ki），这些都可以利用公式生成。具体来说，对于算法2中的延迟-多普勒域MRC操作，设Ki为每个向量νml（或同一延迟抽头具有不同多普勒频移的路径）中的非零信道系数的个数，使OTFS接收机看到的信道系数或传播路径的总数等于所有Ki的和。

Reduced complexity delay-time domain implementation

在上一节中对于每个符号向量xm，我们需要计算l∈L的所有向量bl m。现在提出了一种基于MRC的简化复杂度算法，其算法基本流程如下：

同样的，在算法3中，初始计算包括生成H的所有条目，这需要计算ν对于m=0，...，M‘−1和l∈L，以及不同传播路径下的信道参数 (hi,li,ki)。另一方面，对于分数多普勒情况下，初始计算的复杂性不受影响的延迟时域探测器˜νm，l可以直接从通道增益，延迟，和不同路径下的多普勒转移（hi，li，ki）得出

算法优化

在上面两个算法中，我们均在初始条件下对于所有m使用接下来，我们考虑使用一个更好的OTFS符号向量的初始估计，而不是xm = 0N，来减少MRC迭代的次数并提高了检测和解码的收敛性。假设采用理想的脉冲整形波形，我们在时频域采用单抽头均衡器来提供一个改进的低复杂度的初始估计。

我们设H_dd∈M×N为带有元件的理想脉冲成形波形的德拉多普勒信道矩阵

在这里我们使用的是满足3GPP标准的DD域信道

function [chan_coef,delay_taps,Doppler_taps,taps]=Generate_delay_Doppler_channel_parameters(N,M,car_fre,delta_f,T,max_speed)
one_delay_tap = 1/(M*delta_f);
one_doppler_tap = 1/(N*T);%delays = [0 30 70 90 110 190 410]*10^(-9);%EPA model
delays = [0 30 150 310 370 710 1090 1730 2510]*10^(-9);%EVA model
% delays = [0 50 120 200 230 500 1600 2300 5000]*10^(-9);%ETU modeltaps = length(delays);% number of delay taps
delay_taps = round(delays/one_delay_tap);%assuming no fraction for the delay%pdp = [0 -1.0 -2.0 -3.0 -8.0 -17.2 -20.8];% EPA power delay profile
pdp = [0 -1.5 -1.4 -3.6 -0.6 -9.1 -7.0 -12.0 -16.9];% EVA power delay profile
%pdp = [-1 -1 -1 0 0 0  -3 -5 -7];% ETU power delay profilepow_prof = 10.^(pdp/10);
pow_prof = pow_prof/sum(pow_prof);%normalization of power delay profile
chan_coef = sqrt(pow_prof).*(sqrt(1/2) * (randn(1,taps)+1i*randn(1,taps)));%channel coef. for each path
max_UE_speed = max_speed*(1000/3600);
Doppler_vel = (max_UE_speed*car_fre)/(299792458);
max_Doppler_tap = Doppler_vel/one_doppler_tap;
Doppler_taps = (max_Doppler_tap*cos(2*pi*rand(1,taps)));%Doppler taps using Jake's spectrum
end

通过对分层多普勒信道进行ISFFT运算，得到了理想脉冲整形波形对应的时频信道响应

function [H_t_f]=Generate_time_frequency_channel_OTFSvariants(N,M,gs,L_set,L_cp,variant)
H_t_f=zeros(N,M); % Time-frequency single tap channel matrix
Fm=dftmtx(M);
Fm=Fm./norm(Fm);
Gn=zeros(M,M);
for n=1:Nfor m=1:Mfor l=L_set+1if(strcmp(variant,'CP'))Gn(m,mod(m-l,M)+1)=gs(l,m+(n-1)*(M+L_cp));elseif(m>=l)Gn(m,m-l+1)=gs(l,m+(n-1)*M);  %equation(42) in [R1]endendendendH_t_f(n,1:M)=diag(Fm*Gn*Fm').';
end
end

%References
%  [R1]. Y. Hong, T. Thaj, E. Viterbo, ``Delay-Doppler Communications: Principles and Applications'', Academic Press, 2022, ISBN:9780323850285%----------------------------------------------------------------------------------------
%                                MRC function
%----------------------------------------------------------------------------------------function [est_bits,ite,x_data] = MRC_delay_time_detector_CP(N,M,M_mod,no,data_grid,Y_tilda,H_tf,n_ite_MRC,omega,r,Fn,decision,L_set,nu_ml_tilda,init_estimate)
%% Initial assignments
%Number of symbols per frame
N_syms_perfram=sum(sum((data_grid>0)));
%Arranging the delay-Doppler grid symbols into an array
data_array=reshape(data_grid,1,N*M);
%finding position of data symbols in the array
[~,data_index]=find(data_array>0);
M_bits=log2(M_mod);
N_bits_perfram = N_syms_perfram*M_bits;
L_set=L_set+1;
error=zeros(n_ite_MRC);%% initial estimate using single tap TF equalizer
if(init_estimate==1)Y_tf=fft(Y_tilda).'; % delay-time to frequency-time domain                                                                     % Section 6.5.5 in Chapter 6, [R1]                                  X_tf=conj(H_tf).*Y_tf./(H_tf.*conj(H_tf)+no); % single tap equalizer                                                            X_est = ifft(X_tf.')*Fn; % SFFT                                                                                                 X_est=qammod(qamdemod(X_est,M_mod,'gray'),M_mod,'gray');X_est=X_est.*data_grid;X_tilda_est=X_est*Fn';
elseX_est=zeros(M,N);X_tilda_est=X_est*Fn';
end
x_m=X_est.';
x_m_tilda=X_tilda_est.';%% MRC detector    %% Algorithm 2 in [R1] (or Algotithm 3 in Chapter 6, [R2])
%% initial computation
d_m_tilda=zeros(N,M);
y_m_tilda=reshape(r,M,N).';
delta_y_m_tilda=y_m_tilda;
for m=1:Mfor l=L_setd_m_tilda(:,m)=d_m_tilda(:,m)+abs(nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l).^2);                              % equation (6.49) in Chapter 6, [R1]end
end
for m=1:Mfor l=L_set       delta_y_m_tilda(:,m)=delta_y_m_tilda(:,m)-nu_ml_tilda(:,m,l).*x_m_tilda(:,mod(m-(l-1)-1,M)+1);         % Line 6 of Algorithm 3 in Chapter 6, [R1]end
end
x_m_tilda_old=x_m_tilda;
c_m_tilda=x_m_tilda;%% iterative computation
for ite=1:n_ite_MRCdelta_g_m_tilda=zeros(N,M);for m=1:Mfor l=L_setdelta_g_m_tilda(:,m)=delta_g_m_tilda(:,m)+conj(nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l)).*delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1);    % Line 9 of Algorithm 3 in Chapter 6, in [R1]endc_m_tilda(:,m)=x_m_tilda_old(:,m)+delta_g_m_tilda(:,m)./d_m_tilda(:,m);                                                             % Line 10 of Algorithm 3 in Chapter 6, in [R1]if(decision==1)x_m(:,m)=qammod(qamdemod(Fn*(c_m_tilda(:,m)),M_mod,'gray'),M_mod,'gray');                                                       % Line 11 of Algorithm 3 in Chapter 6, in [R1]x_m_tilda(:,m)=(1-omega)*c_m_tilda(:,m)+omega*Fn'*x_m(:,m);elsex_m_tilda(:,m)=c_m_tilda(:,m);endfor l=L_set                                                                                                                         % Line 12 of Algorithm 3 in Chapter 12, in [R1]delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1)=delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1)-nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l).*(x_m_tilda(:,m)-x_m_tilda_old(:,m));  % Line 13 of Algorithm 3 in Chapter 6, in [R1]end                                                                                                                                 % Line 14 of Algorithm 3 in Chapter 6, in [R1]x_m_tilda_old(:,m)=x_m_tilda(:,m);end%% convergence criteriaerror(ite)=norm(delta_y_m_tilda);if(ite>1)if(error(ite)>=error(ite-1))break;  % Line 16 of Algorithm 3 in Chapter 6, in [R1]endend
end
if(n_ite_MRC==0)ite=0;
end
%% detector output bits
X_est=(Fn*x_m_tilda).';
x_est=reshape(X_est,1,N*M);
x_data=x_est(data_index);
est_bits=reshape(qamdemod(x_data,M_mod,'gray','OutputType','bit'),N_bits_perfram,1);  % Line 18 of Algorithm 3 in Chapter 6, in [R1]end

接下来我们逐段分析代码

第一段

%% Initial assignments
%Number of symbols per frame
N_syms_perfram=sum(sum((data_grid>0)));%N_syms_perfram=64*64=4096 symbols
%Arranging the delay-Doppler grid symbols into an array
data_array=reshape(data_grid,1,N*M);% generate 4096*1 array
%finding position of data symbols >0 in the array
[~,data_index]=find(data_array>0);
M_bits=log2(M_mod);
N_bits_perfram = N_syms_perfram*M_bits;
L_set=L_set+1;%L_set is the unique delay taps =0,1,2
error=zeros(n_ite_MRC);% initialize error; n_ite_MRC：maximum number of MRC detector iterations

第二段

%% initial estimate using single tap TF equalizer
% Alternatively, the initial x estimate can be computed as described in Section 6.5.5.
if(init_estimate==1) % 1-use the TF single tap estimate as the initial estimate for MRC detection, 0-initialize the symbol estimates to 0 at the start of MRC iterationY_tf=fft(Y_tilda).'; % delay-time to frequency-time domain see Section 6.5.5 in Chapter 6 (6.52)                                 X_tf=conj(H_tf).*Y_tf./(H_tf.*conj(H_tf)+no); % single tap equalizer see Section 6.5.5 in Chapter 6 (6.53)                                                             X_est = ifft(X_tf.')*Fn; % SFFT                                                                                                 X_est=qammod(qamdemod(X_est,M_mod,'gray'),M_mod,'gray');X_est=X_est.*data_grid;X_tilda_est=X_est*Fn';  %see Section 6.5.5 in Chapter 6 (6.54)
elseX_est=zeros(M,N);   %Line 1 of Algorithm 2 in Chapter 6, [R1], all the initial x equals to 0X_tilda_est=X_est*Fn';  %all equals to 0
end
x_m=X_est.';
x_m_tilda=X_tilda_est.';

第二段中代码对应公式如下

在理想的脉冲整形波形的情况下，由于信道和延迟-多普勒域中的传输符号的圆卷积在时频域中转换为元素级乘积，我们可以每一个抽头下的MMSE均衡器来估计时频域内的传输样本，如下图所示：

其中 m = 0,...,M − 1 and n = 0,...,N − 1，上图中其他信道函数表示如下：

通过对分层多普勒信道进行ISFFT运算，得到了理想脉冲整形波形对应的相应时频信道响应

最后

第三段对应代码如下

%% MRC detector    %% Algorithm 2 in [R1] (or Algotithm 3 in Chapter 6, [R2])
%% initial computationd_m_tilda=zeros(N,M);   % Algorithm 3 line 1
y_m_tilda=reshape(r,M,N).';     % Algorithm 3 line 2
delta_y_m_tilda=y_m_tilda;
for m=1:Mfor l=L_setd_m_tilda(:,m)=d_m_tilda(:,m)+abs(nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l).^2);    % equation (6.49) in Chapter 6, [R1]end
end
for m=1:Mfor l=L_set  %1: max iterations     delta_y_m_tilda(:,m)=delta_y_m_tilda(:,m)-nu_ml_tilda(:,m,l).*x_m_tilda(:,mod(m-(l-1)-1,M)+1);         % Chapter 6 6.5.3 (6.40)end
end
x_m_tilda_old=x_m_tilda;    % Line 10 of Algorithm 3 in Chapter 6, [R1]
c_m_tilda=x_m_tilda;

接下来就是最激动人心的MRC核心代码：

%% iterative computation
for ite=1:n_ite_MRC delta_g_m_tilda=zeros(N,M); % Algorithm 3 line 9for m=1:Mfor l=L_setdelta_g_m_tilda(:,m)=delta_g_m_tilda(:,m)+conj(nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l)).*delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1);    % Line 9 of Algorithm 3 in Chapter 6, in [R1]endc_m_tilda(:,m)=x_m_tilda_old(:,m)+delta_g_m_tilda(:,m)./d_m_tilda(:,m);                                                             % Line 10 of Algorithm 3 in Chapter 6, in [R1]if(decision==1) %1-hard decision, 0-soft decisionx_m(:,m)=qammod(qamdemod(Fn*(c_m_tilda(:,m)),M_mod,'gray'),M_mod,'gray');                                                       % Line 11 of Algorithm 3 in Chapter 6, in [R1]x_m_tilda(:,m)=(1-omega)*c_m_tilda(:,m)+omega*Fn'*x_m(:,m);% set omega to a smaller value (for example: 0.05) for modulation orders greater than 64-QAMelsex_m_tilda(:,m)=c_m_tilda(:,m);%soft decisionendfor l=L_set                                                                                                                         % Line 12 of Algorithm 3 in Chapter 12, in [R1]delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1)=delta_y_m_tilda(:,mod(m+(l-1)-1,M)+1)-nu_ml_tilda(:,mod(m+(l-1)-1,M)+1,l).*(x_m_tilda(:,m)-x_m_tilda_old(:,m));  % Line 13 of Algorithm 3 in Chapter 6, in [R1]end                                                                                                                                 % Line 14 of Algorithm 3 in Chapter 6, in [R1]x_m_tilda_old(:,m)=x_m_tilda(:,m);end%% convergence criteriaerror(ite)=norm(delta_y_m_tilda);if(ite>1)if(error(ite)>=error(ite-1))break;  % Line 16 of Algorithm 3 in Chapter 6, in [R1] EXITendend
end

最后，输出估计符号x

%% detector output bits
X_est=(Fn*x_m_tilda).';
x_est=reshape(X_est,1,N*M);
x_data=x_est(data_index);
est_bits=reshape(qamdemod(x_data,M_mod,'gray','OutputType','bit'),N_bits_perfram,1);
% Line 18 of Algorithm 3 in Chapter 6, in [R1]

%References
% [R1]. Y. Hong, T. Thaj, E. Viterbo, ``Delay-Doppler Communications: Principles and Applications'', Academic Press, 2022, ISBN:9780323850285