## Copyright (C) 2017 Stephen Bates
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## Author: Stephen Bates <batesste@ennis.local>
## Created: 2017-08-20

##
## sdc_ecc.m
## ---------
##
## This simple script generates the "Bates Conjecture" curve that shows what
## the media Raw Bit Error Rate (RBER) needs to be to reach a target 
## Uncorrectable Bit Error Rate (UBER) as specified by the user.
##
## Basically the code works by assuming a constant code rate (eCodeRate) across
## a range of ECC codeword sizes and determines (via a simple binary search) the
## media RBER needed to target the provided UBER. We assume ECC is done via 
## a Hamming code type structure and the symbol size and correction capabilities
## are computed accordingly.

close all; clear all; clc

eCodeRate = 0.9;   # Code rate, should be between 0 and 1.
eStartBer = 1e-3;  # The start RBER for the binary search
eUber     = 1e-18; # The target UBER.
pnN       = 64:32:8*4096;
                   # The ECC codeword sizes to compute over


  # Calculate the K, M and T for the ECC codewords based on Hamming code 
  # type assumptions.
                   
pnK = floor(eCodeRate*pnN);
pnM = floor(log2(pnN));
pnT = floor((pnN.-pnK)./pnM);

  # Now we enter a search loop for each ECC codeword size using the 
  # incremental Beta function to determine what RBER results in the target 
  # UBER for that size. Stop when we get within 1% of eUber.

peProbFecError = [];
peInputBer     = [];

for i=1:length(pnT)
  
  eOutputBer = 0.5;
  eInputBer  = 24/1080/8;
  bGoingDown = true;
  eScale     = 0.5;
  eTol       = 0.01;
  bUseFer    = true;

  while abs((eOutputBer-eUber)/eUber)>eTol

    eOutputFer = betainc(eInputBer, pnT(i)+1, pnN(i)-pnT(i)-1);
    
    if bUseFer
        eOutputBer = eOutputFer;
    end

    if eOutputBer>eUber
        if ~bGoingDown
            eScale     = (1+eScale)/2;
            eScale     = 1/eScale;
            bGoingDown = true;
        end
        eInputBer = eInputBer*eScale;
    else
        if bGoingDown
            eScale     = (1+eScale)/2;
            eScale     = 1/eScale;
            bGoingDown = false;
        end
        eInputBer = eInputBer*eScale;
    end
  end

  peInputBer = [ peInputBer eInputBer ];
  
  
end

  # Plot the results with Bytes as the X-Axis. Also save to a PNG (this line 
  # works in Octave by might break Matlab).
  
figErr = figure();
hold on ; grid on ; zoom on
loglog(pnN./8, peInputBer);
ylabel('Required Media RBER')
xlabel('Media Codeword Size (Bytes)')
title('Media RBER vs Codeword Size for 1e-18 UBER')
print(figErr,"sdc_ecc.png","-dpng");