how to plot BER curve for VLC-OFDM

Pessoal.

Estou fazendo o script VLC-OFDM baseado em BPSK Watermark.

Alguém poderia me ajudar a encontrar a curva de BER para o código listado abaixo?

Aceito sugestões.

% Número de bits
N = 1e4;     

% tamanho da constelação
M = 2;

% Número máximo de iterações para um único SNR
max_run = 100;              

% Relação Sinal-Ruído (em dB)
SNRdB = 0:12;                             
SNR = 10. ^ (SNRdB / 10);    

% desvio padrão do sinal (potência)
sigmaX2 = 0,15; 

% Modelo LED
imax = 0,5; 

% Modelo LED
K = 2; 

% desvio padrão do ruído (potência)
% sigmaV2 = 0,1;
limiar = 150;
thresholdv = linspace (50,170,20);
thresholdv = thresholdv (fim: -1: 1);
sigmaV2 = (imax / limite) ^ 2; 

SNRL1 = 10 * log10 (sigmaX2 / sigmaV2);

% Modelo LED
Vt = 0,5;

% Modelo LED
Is = 0,12;

sigmaX2v = logspace (log10 (0,001), log10 (0,5), 20);
indsigmaX2 = 1;
indsigmaV2 = 1;


 para contagem = 1: comprimento (SNRL1)                  
     
% para limiar = limiarv       
%     
% sigmaV2vmedian (indsigmaV2) = sigmaV2;
%     
%     SNRideal(indsigmaV2) = 10*log10(sigmaX2/sigmaV2);
% 
%     indsigmaV2 = indsigmaV2+1;
% 
%      end
     
     
%      
%      for sigmaX2=sigmaX2v
%     
%          sigmaX2vmedian(indsigmaX2) = mean(abs( sqrt(sigmaX2)*randn(N,1)));      
%     
%          indsigmaX2 = indsigmaX2+1;
%             
%      end

     
     
     for run_time=1:max_run          
                
         % erro inicial
         Error=0;
                  
         % Gera bits aleatórios ( matriz de variaveis binarias)
         Bits_ale = randsrc(1,N,[0:1]);
         
         % Modulação BPSK 0 -> -1; 1 -> 1
         Bits_bpsk = 2 * Bits_ale-1;  
         
         % Gera bits de marca d'água aleatórios
         Bit_wat = randsrc(1,N,[0:1]);

         % Ângulo da marca d'água
         Theta = pi/4;


         % Fórmula para obter a equação geral.
         for k = 1:N 
      
             if Bit_wat(k)==1        
        
                 Bit_enviado(k) = Bits_bpsk(k) * exp(1i * Bits_bpsk(k) * Theta);
        
             else
                                  
                 Bit_enviado(k) = Bits_bpsk(k) * exp(-1i * Bits_bpsk(k) * Theta);    
    
             end             
             
         end
         
                       
            
         Y = Bit_enviado; 
         
         stdS = std(Y);
         
         X = zeros(1,4*N);
         
         X(2:2:end) = [Y; conj(Y(end:-1:1))];
         
         xTD = ifft(X);
         
         xTD = (sqrt(sigmaX2))*xTD/std(xTD);
         
         x= xTD; 
         
         ind = find(xTD<0);
         
         x(ind) = zeros(size(ind)); 
         
         % Laser-related equations
         s = x;
         
         iLED = Is*(exp(s/Vt)-1);
         xLED = iLED./(  (1 + (iLED/imax).^(2*K)).^(1/(2*K)));
         
         nErrors = 0;
         cont = 0;
         
         while nErrors < 200    
     
             cont=cont+1;
         
             v = sqrt(sigmaV2)*randn(size(x));

    
              y = xLED + v;

              Nframe = length(y);
                
              Xhat = (fft(y));
              
              Shat = Xhat(2:2:Nframe/2);
                
              Shat = stdS*Shat/std(Shat);
                
              RxData1 = pskdemod (Shat, M);
              RxData = Shat;
     
     
             nErrorsFrame = soma (sinal (abs (RxData1 - Bits_bpsk)));     
     
             nErrors = nErrors + nErrorsFrame;   
        

          fim
         
         
         BER = nErros / (N * cont);        
                  
    
     fim                                  
     

% Fim do loop para SNR diferente
 
 fim                                      

 % Calcular BER analítico (para função normal
BER_th = (1/2) * erfc (sqrt (SNR));