OFDM波形探讨
对于通信系统来说,最主要的是解决两个问题:第一是性能的评估标准,也就是信噪比。第二是实现的具体细节,包括编码、调制、均衡等等。

所以首先对信噪比加噪的方法进行探讨,因为加噪的过程直接影响到最后的结果对不对。并且需要明确采样、噪声功率等概念。matlab有两个加噪函数:wgn和awgn。当发射信号的能量为1时,两种加噪都是一致的。

N = 1024;
Tx = ones(1,N);
snr = 10;
Rx = Tx + wgn(1,length(Tx),- snr ,'complex');
Rx = awgn(Tx,snr,'measured');
cau_P(Tx)
cau_P([ Tx , zeros(1,N)])
function P = cau_P(sig)
P = sum(abs(sig).^2)/length(sig);
end
function SNR = cau_SNR(sig,noise)
P1 = sum(abs(sig).^2)/length(sig);
P2 = sum(abs(noise).^2)/length(noise);
SNR = 10*log10(P1 ./ P2);
end

信噪比的理解如下:

参考博客里面有对信噪比比较明确的解释。

在通信系统的接收端,噪声会随着信号一起进入接收机,这时就会判断在信噪比为多少的情况下误码率是多少,这时SNR、Eb/N0、Es/N0都可能用到。SNR 也即信噪比,是接收端模拟信号的重要测量指标,可以通过频谱仪等仪器实际测量接收端的模拟信号得到。而Eb/N0 是指通信系统传输一比特信息所需要的能量和噪声功率谱密度的比值,是衡量整个通信系统性能归一化的一个系统指标。一般情况下,模拟通信系统常采用SNR-BER来衡量通信系统的性能,数字通信系统常采用Eb/N0-BER 来衡量通信系统的性能。由于数字通信系统常采用Eb/N0 作为衡量系统的性能指标,但实际测量Eb/N0 比较困难,故而研究Eb/N0 与SNR 的转化关系变的十分重要。

  • 通过过采样可以提升信噪比SNR,从而系统得到更好的性能;
  • EbNo与EsNo的关系是固定的,只与编码和调制方式有关。当R为4时,说明一个符号表示4个bit,因此EbNo也就越低了。(相当于4bit平分一个符号的能量)
  • 一般对比性能都是用EbN0
close all;clear all;
% rng default

%**********************************
% 仿真参数
SNR_dBlist  = (0:0.25:1.0);
% SNR_dBlist  = (-5.0:0.25:-3.5);
% SNR_dBlist  = (1:0.25:2);
% SNR_dBlist  = 7;
% SlotBitLen  = 1008;
SlotBitLen  = 1344;

%------------------------------------
%2/3编码 
Length = ceil(SlotBitLen/(2/3)/4)*4;
Length1 = 1024;
%------------------------------------
Num1 = 1024;
Num2 = 2048;
fs   = 2e5;          %   码片速率,单位kHz
delta_f = 0;         %   频偏残差,单位kHz
% delta_f = 0;       %   频偏残差,单位kHz
L_delay = 16;        %   
QM = 2;
% --new-add------------------------------------
% ************************************
% SlotBitLen = 5760;
% Length = ceil(3/2*SlotBitLen/2)*2;  %  是2的倍数
%------------------------------------
SubInterBitLen = SlotBitLen + 4;  % 子交织器输入的个数
ColNum      = 32;                 % 行列交织器的列数
RowNum      = ceil(SubInterBitLen/ColNum); % 行列交织器的行数
SubMatEleNum      =  ColNum*RowNum;
PadONum           =  SubMatEleNum - SubInterBitLen;%补充0
ColIndex            = [0  16  8  24  4  20  12  28  2  18  10  26  6  22  14  30  1  17  9  25  5  21  13  29  3  19  11  27  7  23  15  31 ];  % 列交织器
TCRvIdx         = 0;   
RateMatchLenList    = Length;
QPSK_func = 1/sqrt(2)*[1+1i,1-1i,-1+1i,-1-1i];
% -------打印------------------------------
PS = 'delay';
PS = [PS,'\nQPSK 2/3  '];
%***************************************************
[f1, f2]     = getf1f2(SlotBitLen);
InterTbl    = zeros(1,SlotBitLen);
for ii = 0:SlotBitLen-1
   InterTbl(ii+1) = mod(f1*ii+f2*ii^2,SlotBitLen);
end
hTEnc           = comm.TurboEncoder(poly2trellis(4, [13 15], 13), InterTbl + 1);                          % 简化方式
hTDec           = comm.TurboDecoder(poly2trellis(4, [13 15], 13), InterTbl + 1, 10,'Algorithm','Max');    % 简化方式
%**************************************************
tic
Pe_Block = zeros(1,length(SNR_dBlist));err_rate = zeros(1,length(SNR_dBlist));

RcosCoef        = rcosine(1,8,'fir',0.35,6);%97阶
SmplBias        = 6;            % 决定了码片的偏差为SmplBias/8
if (SmplBias==0)
   h_list_1(1:13)    = RcosCoef(SmplBias+1:8:end); 
else
   h_list_1(1:12)    = RcosCoef(SmplBias+1:8:end); 
end
% fid = PRINT_comp(0,PS,0,delta_f,L_delay, SmplBias); %打印
for SnrNum = 1:length(SNR_dBlist)               
   error = 0;Numbit = 0;BlockNum = 0; BlockErr = 0; SNR = SNR_dBlist(SnrNum);
   while( ( (BlockErr < 100)||(BlockNum < 500) ) && (BlockNum < 5000) )
%       for zzzz = 1:100
       %------------------------------------------------------------------
       % 生成基本数据源
       BlockNum   = BlockNum + 1;
       Numbit     = Numbit +  SlotBitLen ;                 % 每次的仿真点数增加SrcLen
       %------------------------------------
       TCSrc           = randi([0 1],1,SlotBitLen);
       % Turbo编码
       TurboCodeOut    = step(hTEnc, TCSrc.').';   % 1/3码率 X0A0B0X1A1B1
% in 4044
       %% -----交织----------------------------------------
       TCBuffMat    = zeros(1,3*SubMatEleNum);
       SubInterMat     = zeros(3,SubMatEleNum);            % 子交织矩阵
       % 交织
       for jj = 1:3
           % 提取第jj路数据
           SubTCSrc        = TurboCodeOut(jj:3:end);
           SigPad          = [NaN*ones(1,PadONum),SubTCSrc];       % LTE标准中以Null进行比特填充, 仿真中以NaN代替
           if (jj==3)
               SigPad      = [SigPad(2:end),SigPad(1)];    % 第二路校验比特要错开1比特
           end
           % 进行子块交织, 行进列出的过程
           SigPadMat       = reshape(SigPad, ColNum, RowNum).';    % 行入
           SigMatColInter  = SigPadMat(:,ColIndex + 1);            % 交织
           SubInterMat(jj,:)  = reshape(SigMatColInter,1,[]);      % 列出
       end
       % 输出得到 3*SubMatEleNum 长度的序列(或者取其他的长度)
       TCBuffMat(1:SubMatEleNum)                     = SubInterMat(1,:);
       TCBuffMat(SubMatEleNum + 1:2:3*SubMatEleNum)  = SubInterMat(2,:);
       TCBuffMat(SubMatEleNum + 2:2:3*SubMatEleNum)  = SubInterMat(3,:);
% 4128 out

       % ------速率匹配----------------------------------------

       TCKw           = 3*SubMatEleNum;             % 理想循环缓存器总长度
       TCNcb          = TCKw;
       % 根据冗余版本号 TCRvIdx 和  TCNcb 确定开始位置
       RMBeginPos     = RowNum*(2*ceil(TCNcb/8/RowNum)*TCRvIdx + 2);
       % 循环读取(mod TCNcb)直到匹配的长度
       kk          = 0;
       jj          = 0;
       MatchOut    = zeros(1,RateMatchLenList);
       while (kk<RateMatchLenList)
           % 候选的比特
           CandiBit       = TCBuffMat( mod(RMBeginPos + jj, TCNcb) + 1);
           if (~isnan(CandiBit))
               MatchOut(kk+1) = CandiBit;
               kk   = kk + 1;
           end
           % 缓存区位置加1
           jj  = jj + 1;
       end
% 2028 out

       TurboCodeOut1d2 = MatchOut;
       %% ---------------------------------------
       % in 2048 TurboCodeOut1d2
       % SC_FDMA
       Tx = zeros(1,(Num2+L_delay)*Length1/Num1);
%         tmp_index = randperm(Num2,Num1);
       tmp_index = [(1:Num1/2),(Num2-Num1/2+1:Num2)];
       for kk = 1:Length1/Num1              %2048/1024
%             tmp = 1-2*TurboCodeOut1d2(Num1*(kk-1)+1:Num1*kk);   % 映射过程  1024-->1024
           if length(TurboCodeOut1d2)/QM<=Num1
               tmp = QPSK_func(2*TurboCodeOut1d2(1:2:end)+TurboCodeOut1d2(2:2:end)+1);
           else
               error('长度超了');
           end
           tmp_fft = 1/sqrt(Num1)*fft(tmp,Num1);
           tmp_fft1 = zeros(1,Num2);
           tmp_fft1(tmp_index) = tmp_fft;
           tmp1 = sqrt(Num2^2/Num1)*ifft(tmp_fft1,Num2);        %%能量归1化
           Tx((Num2+L_delay)*(kk-1)+1:(Num2+L_delay)*kk) = [tmp1(end-L_delay+1:end),tmp1];
       end
       % out 4128  Tx
       %---------------------------------------
       % 发送信号
       h_list = zeros(1,L_delay);
%         h_list(1:5) = [0.29,0.5,0.58,0.5,0.29];
%         h_list(1:5) = [1,0,0,0,0];
       % 这里加入了码片的偏差
       h_list(1:length(h_list_1)) = h_list_1;%这个应该对应的就是信道的冲击响应
       % ------------调制和加入噪声----------------------------------------------------
       SigChann = zeros(1,length(Tx));
       for kk = 1:L_delay
           SigChann = SigChann + h_list(kk)*[zeros(1,kk-1),Tx(1:end+1-kk)];%%其实相当于是卷积?
       end
       SigChann  = exp(1j*2*pi*delta_f/fs*[(1:length(SigChann)/2),(1:length(SigChann)/2)]).*SigChann + wgn(1,length(SigChann),-SNR_dBlist(SnrNum),'complex');
       %------------------------------------      
       %  收端信号
       % 这里用来修改信道冲激响应
       % 是不是这里模拟的是导频能否进行信道的估计呢?
       % 有没有导频进行信道估计会差
%         h_est = [zeros(1,6),1];
       h_est = h_list_1;
       H_est = fft(h_est,Num2);
%         figure;plot(abs(H_est),'r');grid on;
%         figure;plot(20*log10(abs(fftshift(H_est))),'r');grid on;
       % MMSE均衡
       % 说白了就是用MMSE方法得到2048个点的频率响应以便下面的均衡;
       % 但是,如果没有导频的话,这里的信道估计是不够准确的,就模拟出了导频不进行信道估计的情况;
       H_est_inv = conj(H_est)./(H_est.*conj(H_est)+10^(-SNR_dBlist(SnrNum)/10));
       Rx = zeros(1,Length1);
       for kk = 1:Length1/Num1
           Rx_tmp = SigChann((Num2+L_delay)*(kk-1)+1+L_delay:(Num2+L_delay)*kk);%2048
           Rx_tmp_fft = fft(Rx_tmp,Num2);%2048
           Rx_tmp_fft1 = Rx_tmp_fft.*H_est_inv;      % 均衡        
           Rx_tmp_fft2 = Rx_tmp_fft1(tmp_index);
           Rx_tmp1 = ifft(Rx_tmp_fft2,Num1);
           Rx(Num1*(kk-1)+1:Num1*kk) = Rx_tmp1;
       end
%         SigPunc = Rx(1:Length);
       % Rx  2148
       % out 2048 SigPunc
       SigPunc = Rx(1:Length/2);

       y_QPSK_soft = zeros(1,Length);
       y_QPSK_soft(1:2:end) = real(SigPunc);
       y_QPSK_soft(2:2:end) = imag(SigPunc);


%         biterr(y_QPSK_soft<0,MatchOut)/2016
       %% --------解速率匹配-------------------------------------------------
%         y_QPSK_soft = 
       RcvRMMat        = zeros(1,3*SubMatEleNum);
       % 初始化接收缓存
       RcvCirBufTmp    = zeros(1,3*SubMatEleNum);

       % 确定系统比特和校验比特序列中Null的位置(这里应该跟交织有关系)
       % 这里为什么选取三组28bit的数列赋NaN值呢?

       RcvCirBufTmp(ColIndex(1:PadONum)*RowNum + 1)  = NaN;
       RcvCirBufTmp(SubMatEleNum + 2*ColIndex(1:PadONum)*RowNum + 1)  = NaN;
       RcvCirBufTmp(SubMatEleNum + 2*ColIndex(1:PadONum - 1)*RowNum + 2)  = NaN;
       if (PadONum >0 )
           RcvCirBufTmp(end) = NaN;
       end
       % 解速率匹配过程
       kk          = 0;
       jj          = 0;
       % 根据冗余版本号 TCRvIdx 和  TCNcb 确定开始位置
       RMBeginPos   = RowNum*(2*ceil(TCNcb/8/RowNum)*TCRvIdx + 2);
       AddTimeList  = zeros(1, 3*SubMatEleNum);              % 累加次数统计
       while (kk<RateMatchLenList)
           % 候选的比特
           CandiBit       = RcvCirBufTmp(mod(RMBeginPos + jj, TCNcb) + 1);
           %              if CandiBit~=0    %可能是28
           %                  fprintf('youbuwei0\n')
           %              end
           if (~isnan(CandiBit))
               RcvCirBufTmp(mod(RMBeginPos + jj, TCNcb) + 1) = y_QPSK_soft(kk+1) + CandiBit;
               %                  AddTimeList(mod(RMBeginPos + jj, TCNcb) + 1)  = AddTimeList(mod(RMBeginPos + jj, TCNcb) + 1) + 1;   % 累积次数加1
               kk   = kk + 1;
           end
           % 缓存区位置加1
           jj  = jj + 1;
       end
       % --------解交织-------------------------------------------------

       RcvCirBuf = RcvCirBufTmp;
       % ******************************************************
       RcvSubBuffMat       = zeros(3,SubMatEleNum);
       RcvSubBuffMat(1,:)  = RcvCirBuf(1:SubMatEleNum);
       RcvSubBuffMat(2,:)  = RcvCirBuf(SubMatEleNum+1:2:3*SubMatEleNum);
       RcvSubBuffMat(3,:)  = RcvCirBuf(SubMatEleNum+2:2:3*SubMatEleNum);
       %  解交织
       TCDecSoftInfo       = zeros(1,3*SubInterBitLen);
       for jj = 1:3
           RcvTCSoftMatTmp   = reshape(RcvSubBuffMat(jj,:),RowNum,ColNum);          % 列入
           RcvTCSoftMat(1:RowNum,ColIndex + 1) = RcvTCSoftMatTmp;                   % 解交织
           RcvSoftDeInter = reshape(RcvTCSoftMat.',1,[]);
           if (jj==3)
               RcvSoftDeInter  = [RcvSoftDeInter(end),RcvSoftDeInter(1:end-1)];   % 第二路校验比特要错开1比特
           end
           % Turbo译码输入软信息
           TCDecSoftInfo(jj:3:end) = RcvSoftDeInter(PadONum+1:end);
       end
       % ------------------------------------------------------
       %% 译码
       TCDDecOut      = step(hTDec, -TCDecSoftInfo.').';        % APP算法是按照0映射成-1设计的, 与实际的映射0--1 相反,所以要取负数
       TCDecErr       = sum(abs(TCSrc - TCDDecOut));

      % 进行当次译码结果的错误分析统计
       error = error + TCDecErr ;                         % 错误数目累加
       BlockErr = BlockErr + (TCDecErr>0);                % 错误分组统计

    end
   %---------------------------------------------------------    
   % 统计所有的错误概率
   err_rate( SnrNum) =  error / Numbit;                   % 错误概率输出
   BlockNum = Numbit/SlotBitLen;                           % 分组总数
   Pe_Block(SnrNum) = BlockErr/BlockNum;                   % 误分组概率输出

  fprintf('EsN0 = %5.2f,error = %d, NumBit = %d , pe = %7.6f,SlotNum = %d,SlotErr = %d,Pe_Slot = %5.4f ,\n',...
      SNR_dBlist(SnrNum), error, Numbit, err_rate(SnrNum) ,BlockNum ,BlockErr,Pe_Block(SnrNum));
%    fprintf(' %5.2f %d %d %7.6f %d %d %5.4f \n',...
%        SNR_dBlist(SnrNum), error, Numbit, err_rate(SnrNum) ,BlockNum ,BlockErr,Pe_Block(SnrNum));

end
toc