测试文献Iterative MAP equalization %表达式(18)-(19)

%6. 测试文献Iterative MAP equalization and decoding in wireless mobile coded ofdm %表达式(18)-(19) clc;clear;clc;M = 8;   %在载波数K = MF = 1/sqrt(M) * fft(eye(M));   % unitary FFT matrix of size ML =   2;       % CP length (same as channel order)Q2 =   1;       %对角化用q  =   Q2;lambda = 0.6; pdp = exp( -lambda*(0:L).' ) / sum( exp( -lambda*(0:L).') ); % for exponential power-delay profilenormalized_Doppler_frequency = 0.05;     % with respect to the subcarrier spacingf_d = normalized_Doppler_frequency / M; bit1=[0 1 1 1 1 1 0 1]';bit2=[0 0 0 1 1 1 1 1]';   %数据源d = 1/sqrt(2) * ((2*bit1-1)+sqrt(-1)*(2*bit2-1));  %星座点n = 1/sqrt(2) * (randn(M,1)+sqrt(-1)*randn(M,1)); % generation of the current channel realization ht = ltv_channel_gen_test(f_d,M,32,L+1,pdp); % time-varying CIR; Ht = zeros(M,M); % time-domain channel matrix (pre-allocation) for index_subcarrier = 1:M      index_col = mod(((index_subcarrier-L):index_subcarrier)-1,M)+1;      Ht(index_subcarrier,index_col) = fliplr(ht(1:(L+1),index_subcarrier).'); % time-domain channel matrix end % frequency-domain channel matrices (with mask to enforce banded matrices) Hf       = F * Ht * (F'); %对角化 mask2 = toeplitz([ones(Q2+1,1);zeros(M-(2*Q2+1),1);ones(Q2,1)]); mask2(M-1:M,1) = 0; mask2(M,2)     = 0; mask2(1,M-1:M) = 0; mask2(2,M)     = 0;  Hf_banded     = Hf.* mask2;        %对角后的矩阵  d_hat         = Hf*d + n;          %信道传输  Hf_banded_hat = Hf_banded';  temp1 = d_hat'* Hf_banded *d + d'*  Hf_banded_hat*d_hat;    %虚部和虚部约掉了，只剩下实部  temp2 = 2*real( d'*  Hf_banded_hat*d_hat);                  %式(18) 第一项  temp3 = d'*Hf_banded'*Hf_banded*d;                          %式(18) 第二项       %第一项  temp_sum = 0;  for m=1:K      %mod(m-q:m+q,K)      temp = 0;      for k=(mod(m-q:m+q,K))                            %matlab下标从1开始，模掉K，再把0置为K          if (k == 0)              k = K;          end          temp = temp +  Hf_banded_hat(m,k)*d_hat(k);      end      temp_d = (d(m))'*2*temp;      temp_sum = temp_sum + real(temp_d);  end  %第二项%   G = Hf_banded_hat*Hf_banded;                     %G%   tempg_sum = 0;%   for m=1:K%       %mod(m-q:m+q,K)%       tempg = 0;%       for k=mod(1:2*q,K)                            %matlab下标从1开始，模掉K，再把0置为K%           if (k == 0)%               k = K;%           end%           %           if (m-k>0&&m-k<K)%           tempg = temp +  G(m,m-k)*d(m-k);%           end%       end%       temp_d = G(m,m)*d(m) + 2*tempg;%       temp_gg = (d(m))'*  temp_d%       tempg_sum = tempg_sum + real(temp_gg);%       %   end  G = Hf_banded_hat*Hf_banded;                           %G  tempg_sum = 0;  for m=1:K      %mod(m-q:m+q,K)      tempg = 0;      for k=mod(m-2*q:m+2*q,K)  %matlab下标从1开始，模掉K，再把0置为K          if (k == 0)              k = K;            %用自己推的和原来的计算结果一致，                                %但是用作者的计算与原表达式不一样          end          tempg = tempg +  G(m,k)*d(k);      end          temp_gg = (d(m))'*  tempg;          tempg_sum = tempg_sum + real(temp_gg);  end