ber_ofdm_64s_1e4_withclipping_snr20_differcr_1_2_3_nofilter.m

限幅方法降低OFDM系统误码率

echo  off;clear all;close all;clc;
fprintf( 'OFDM仿真\n') ;
tic
%无限幅
% --------------------------------------------- %
%                   参数定义                     %
% --------------------------------------------- %
% Initialize the parameters
NumLoop = 10000;
NumSubc = 64;
NumCP = 8;
SyncDelay = 0;
%-----------------------------------------------

%-----------------------------------------------
% Generate the random binary stream for transmit test
%-----------------------------------------------
BitsTx = floor(rand(1,NumLoop*NumSubc)*2);%用rand产生0：1的随机数，不会产生负值，如果用randn会产生负值；乘以2是为了使用floor函数

% --------------------------------------------- %
%                   QAM MODULATION              %
% --------------------------------------------- %
% Modulate (Generates QAM symbols)
% input: BitsTx(1,NumLoop*NumSubc); output:SymQAM(NumLoop,NumSubc/2)
SymQAMtmp = reshape(BitsTx,2,NumLoop*NumSubc/2).';
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');%/?????(bi2de)
QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1);

% --------------------------------------------- %
%                   IFFT   input: SymQAM(NumLoop,NumSubc/2); output: SymIFFT(NumSubc,NumLoop)
% --------------------------------------------- %
SymIFFT = zeros(NumSubc,NumLoop);
SymIFFTtmp = reshape(SymQAM,NumSubc/2,NumLoop)
SymIFFTtmptmp = zeros(NumSubc,NumLoop);
SymIFFTtmptmp(1,:) = real(SymIFFTtmp(1,:)); % 实数
SymIFFTtmptmp(NumSubc/2+1,:) = imag(SymIFFTtmp(1,:)); % 实数
SymIFFTtmptmp(2:NumSubc/2,:) = SymIFFTtmp(2:NumSubc/2,:);
SymIFFTtmptmp((NumSubc/2+2):NumSubc,:) = flipdim(conj(SymIFFTtmp(2:NumSubc/2,:)),1);

SymIFFT = ifft(SymIFFTtmp,NumSubc,1);%(后边这个1是干什么的？)

% --------------------------------------------- %
%             Add cyclic prefix                 %
% --------------------------------------------- %
% input: SymIFFT(NumSubc,NumLoop); output: SymCP(NumSubc + NumCP,NumLoop)
NumAddPrefix = NumSubc + NumCP;
SymCP = zeros(NumAddPrefix,NumLoop);
RowPrefix = (NumSubc - NumCP + 1):NumSubc;
SymCP = [SymIFFT(RowPrefix,:);SymIFFT];%(把保护间隔放在符号的前面)

% --------------------------------------------- %
%             Go through the channel            %
% --------------------------------------------- %
% input: SymCP(NumSubc + NumCP,NumLoop); output: SymCh(1,(NumSubc + NumCP)*NumLoop)
SymChtmp = zeros(1,(NumSubc + NumCP)*NumLoop);
SymChtmp = SymCP(:).';   % 进行这个转置操作之后就成了一个矢量;       并串变换
% 相当于把矩阵的列向量依次排列 改变为一个行向量

%--------------------------------------------------
% Add the AWGN
%--------------------------------------------------
snrIndex = 0;
BerSnrTable_NOCP = zeros(7,2);
for snr=0:2:13;  % = SNR + 10*log10(log2(2));
    snrIndex = snrIndex+1;
    BerSnrTable_NOCP(snrIndex,1) = snr;
    SymCh = awgn(SymChtmp,snr,'measured');

    % --------------------------------------------- %
    %            Remove Guard Intervals             %
    % --------------------------------------------- %
    % input: SymCh(1,(NumSubc + NumCP)*NumLoop); output: SymDeCP(NumSubc,NumLoop)
    SymDeCP = zeros(NumSubc,NumLoop);
    SymDeCPtmp = reshape(SymCh,NumSubc + NumCP,NumLoop);
    SymDeCP = SymDeCPtmp((NumCP+1+SyncDelay):NumAddPrefix+SyncDelay,:);

    % --------------------------------------------- %
    %                     FFT                       %
    % --------------------------------------------- %
    % input: SymDeCP(NumSubc,NumLoop); output: SymFFT(NumSubc,NumLoop)
    SymFFT = fft(SymDeCP,NumSubc,1);

    % --------------------------------------------- %
    %        Make Decision(Include DeQAM)           %
    % --------------------------------------------- %
    % SymFFT(NumSubc,NumLoop); output: SymDec(NumSubc,NumLoop)
    SymDec = zeros(NumSubc,NumLoop);
    SymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:);
    SymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:);%???????????(SymEqtmp是哪里来的？)
    for m = 1:NumLoop
        for n = 1:NumSubc/2
            Real = real(SymEqtmp(n,m));
            Imag = imag(SymEqtmp(n,m));
            if( abs((Real -1)) < abs((Real +1)))%????????每个点到坐标轴上1和-1的距离，哪个离1近就把1付给它，同理得-1
                SymDec(2*n-1,m) = 1;
            else
                SymDec(2*n-1,m) = 0;%把实部放在奇数行
            end

            if( abs((Imag -1)) < abs((Imag +1  )) )
                SymDec(2*n,m) = 1;

            else
                SymDec(2*n,m) = 0;%把虚部放在偶数行
            end
        end
    end
    
    %----------------------------------------------------------------------
    % Bit Error
    %----------------------------------------------------------------------
    BitsRx = zeros(1,NumSubc*NumLoop);
    BitsRx = SymDec(:).';
    [Num,Ber] = symerr(BitsTx,BitsRx)%(BitsTx是原始数据符号)
    %BerSnrTable(snr+1,2) = Num ;
    BerSnrTable_NOCP(snrIndex,2) = Ber ;%BerSnrTable(snr+1,3) = Ber ;
end;
BerSnrTable_NOCP_plot=zeros(7,1);
BerSnrTable_NOCP_plot=BerSnrTable_NOCP(:,2);
save('NO_CR','BerSnrTable_NOCP_plot');
    semilogy(BerSnrTable_NOCP(:,1),BerSnrTable_NOCP_plot,'o-');
    legend('orignal');
    xlabel('SNR');
    ylabel('BER');
    grid on;
    hold on;
    % -------------------------------------------------------CR-----------------
    time_of_sim = toc
% end;




% %CR=1.2
% % --------------------------------------------- %
% %                   参数定义                     %
% % --------------------------------------------- %
% % Initialize the parameters
% % NumLoop = 10000;
% % NumSubc = 64;
% % NumCP = 8;
% % SyncDelay = 0;
% %-----------------------------------------------

%-----------------------------------------------
% Generate the random binary stream for transmit test
%-----------------------------------------------
BitsTx = floor(rand(1,NumLoop*NumSubc)*2);%用rand产生0：1的随机数，不会产生负值，如果用randn会产生负值；乘以2是为了使用floor函数

% --------------------------------------------- %
%                   QAM MODULATION              %
% --------------------------------------------- %
% Modulate (Generates QAM symbols)
% input: BitsTx(1,NumLoop*NumSubc); output:SymQAM(NumLoop,NumSubc/2)
SymQAMtmp = reshape(BitsTx,2,NumLoop*NumSubc/2).';
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');%/?????(bi2de)
%--------------------------------------------------------------------------

QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1);

% --------------------------------------------- %
%                 IFFT   input: SymQAM(NumLoop,NumSubc/2); output: SymIFFT(NumSubc,NumLoop)
% --------------------------------------------- %

SymIFFT = zeros(NumSubc,NumLoop);
SymIFFTtmp = reshape(SymQAM,NumSubc/2,NumLoop);%（为什么载波数除以2）
SymIFFTtmptmp = zeros(NumSubc,NumLoop);
SymIFFTtmptmp(1,:) = real(SymIFFTtmp(1,:)); % 实数
SymIFFTtmptmp(NumSubc/2+1,:) = imag(SymIFFTtmp(1,:)); % 实数

SymIFFTtmptmp(2:NumSubc/2,:) = SymIFFTtmp(2:NumSubc/2,:);
SymIFFTtmptmp((NumSubc/2+2):NumSubc,:) = flipdim(conj(SymIFFTtmp(2:NumSubc/2,:)),1);

SymIFFT = ifft(SymIFFTtmp,NumSubc,1);%(后边这个1是干什么的？)

% --------------------------------------------- %
%             Add cyclic prefix                 %
% --------------------------------------------- %
% input: SymIFFT(NumSubc,NumLoop); output: SymCP(NumSubc + NumCP,NumLoop)
NumAddPrefix = NumSubc + NumCP;
SymCP = zeros(NumAddPrefix,NumLoop);
RowPrefix = (NumSubc - NumCP + 1):NumSubc;
SymCP = [SymIFFT(RowPrefix,:);SymIFFT];%(把保护间隔放在符号的前面)

% --------------------------------------------- %
%                  限幅                         %
% --------------------------------------------- %
CR = 1.2;
SymCP_Clip_last =  zeros(NumAddPrefix,NumLoop);
SymCP_Clip_temp3 = zeros(NumAddPrefix,1);
%SymCP_Clip_temp2 = zeros()

for nSymbol = 1:NumLoop
    Signal_Power = abs(SymCP(:,nSymbol).^2);
    Peak_Power = max(Signal_Power); %Peak_Power = max(Signal_Power,[],2);
    Mean_Power = mean(Signal_Power);%Mean_Power = mean(Signal_Power,2);
    
    A = CR*sqrt(Mean_Power);
    SymCP_Clip_temp1 = SymCP(:,nSymbol)';%把第nSymbol列的值转换成行向量并付给一个向量
    SymCP_Clip_temp2 = SymCP_Clip_temp1(Signal_Power>A);
    SymCP_Clip = A*SymCP_Clip_temp2./abs(SymCP_Clip_temp2);
    SymCP_Clip_temp1(Signal_Power>A) = SymCP_Clip;
    
    SymCP_Clip_temp3 = SymCP_Clip_temp1'
    SymCP_Clip_last(:,nSymbol) = SymCP_Clip_temp3;
end
    
    
% --------------------------------------------- %
%             Go through the channel            %
% --------------------------------------------- %
% input: SymCP(NumSubc + NumCP,NumLoop); output: SymCh(1,(NumSubc + NumCP)*NumLoop)
SymChtmp = zeros(1,(NumSubc + NumCP)*NumLoop);
SymChtmp = SymCP_Clip_last(:).';   % 进行这个转置操作之后就成了一个矢量;       并串变换
% 相当于把矩阵的列向量依次排列 改变为一个行向量

% % --------------------------------------------- %
% %滤波
% % --------------------------------------------- %
% Ch = [1 1/2 1/4];
% SymChtmptmp = filter(Ch,1,SymChtmp);%滤波
snrIndex4 = 0;
BerSnrTable_CR1P2 = zeros(7,2);%?????
for snr=0:2:13;  % = SNR + 10*log10(log2(2));
    snrIndex4=snrIndex4+1;
    BerSnrTable_CR1P2(snrIndex4,1) = snr;
    SymCh = awgn(SymChtmp,snr,'measured');

    % --------------------------------------------- %
    %            Remove Guard Intervals             %
    % --------------------------------------------- %
    % input: SymCh(1,(NumSubc + NumCP)*NumLoop); output: SymDeCP(NumSubc,NumLoop)
    SymDeCP = zeros(NumSubc,NumLoop);
    SymDeCPtmp = reshape(SymCh,NumSubc + NumCP,NumLoop);
    SymDeCP = SymDeCPtmp((NumCP+1+SyncDelay):NumAddPrefix+SyncDelay,:);

    % --------------------------------------------- %
    %                     FFT                       %
    % --------------------------------------------- %
    % input: SymDeCP(NumSubc,NumLoop); output: SymFFT(NumSubc,NumLoop)
    SymFFT = fft(SymDeCP,NumSubc,1);

    % --------------------------------------------- %
    %        Make Decision(Include DeQAM)           %
    % --------------------------------------------- %
    % SymFFT(NumSubc,NumLoop); output: SymDec(NumSubc,NumLoop)
    SymDec = zeros(NumSubc,NumLoop);
    SymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:);
    SymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:);%???????????(SymEqtmp是哪里来的？)
    for m = 1:NumLoop
        for n = 1:NumSubc/2
            Real = real(SymEqtmp(n,m));
            Imag = imag(SymEqtmp(n,m));
            if( abs((Real -1)) < abs((Real +1)))%????????每个点到坐标轴上1和-1的距离，哪个离1近就把1付给它，同理得-1
                SymDec(2*n-1,m) = 1;
            else
                SymDec(2*n-1,m) = 0;%把实部放在奇数行
            end

            if( abs((Imag -1)) < abs((Imag +1  )) )
                SymDec(2*n,m) = 1;

            else
                SymDec(2*n,m) = 0;%把虚部放在偶数行
            end
        end
    end
    %----------------------------------------------------------------------
    % Bit Error
    %----------------------------------------------------------------------
    BitsRx = zeros(1,NumSubc*NumLoop);
    BitsRx = SymDec(:).';
    [Num,Ber] = symerr(BitsTx,BitsRx)%(BitsTx是原始数据符号)
    %BerSnrTable(snr+1,2) = Num ;
    BerSnrTable_CR1P2 (snrIndex4,2) = Ber ;%BerSnrTable(snr+1,3) = Ber ;
end;
BerSnrTable_CR1P2_plot = zeros(7,1);
BerSnrTable_CR1P2_plot=BerSnrTable_CR1P2(:,2);
save('CR_1P2','BerSnrTable_CR1P2_plot');
    semilogy(BerSnrTable_CR1P2(:,1),BerSnrTable_CR1P2_plot,'+-r');
    legend('CR=3');
%     xlabel('SNR');
%     ylabel('BER');
    grid on;
    hold on;
    % ------------------------------------------------------------------------
    time_of_sim = toc



%CR=0.8
% --------------------------------------------- %
%                   参数定义                     %
% --------------------------------------------- %
% Initialize the parameters
% NumLoop = 10000;
% NumSubc = 64;
% NumCP = 8;
% SyncDelay = 0;
%-----------------------------------------------

%-----------------------------------------------
% Generate the random binary stream for transmit test
%-----------------------------------------------