Reed Solomon(リード ソロモン符号)を実装します。
今回は実装のみで理論の解説はありません。
また、コードも流用コードを含みます。
使い方はこちら
こちらを参考にしております。理論などの解説も下記書籍からどうぞ。
ガロア体原始多項式実装ユニット
unit iPentecGalois;
interface
type
TiPentecGalois = class
private
PP : Array[0..16] of Byte ; { specify irreducible polynomial coeffts }//ガロア体原始多項式
pPP : Array[2..16] of Pointer ;
FTableLength:integer;
//procedure init_exp_table();
procedure init_exp_table(Dimension:integer);
procedure SetPrimitive(var PP; nIdx : Integer ) ;
public
gexp: array of integer;
glog: array of integer;
//gexp: array [0..512] of integer;
//glog: array [0..512] of integer;
constructor Create();
procedure InitTables(Dimension:integer);
function gmult(a:integer; b:integer):integer;
function ginv (elt:integer):integer;
end;
implementation
const
//ガロア体原始多項式
PP2 : Array[0..2] of Byte = ( 1 , 1 , 1 ) ;
PP3 : Array[0..3] of Byte = ( 1 , 1 , 0 , 1 ); { 1 + x + x^3 }
PP4 : Array[0..4] of Byte = ( 1 , 1 , 0 , 0 , 1 ) ; { 1 + x + x^4 }
PP5 : Array[0..5] of Byte = ( 1 , 0 , 1 , 0 , 0 , 1 ) ; { 1 + x^2 + x^5 }
PP6 : Array[0..6] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 1 ) ; { 1 + x + x^6 }
PP7 : Array[0..7] of Byte = ( 1, 0, 0, 1, 0, 0, 0, 1 ) ; { 1 + x^3 + x^7 }
PP8 : Array[0..8] of Byte = ( 1 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x^2+x^3+x^4+x^8 }
PP9 : Array[0..9] of Byte = ( 1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^4+x^9 }
PP10 : Array[0..10] of Byte = ( 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^3+x^10 }
PP11 : Array[0..11] of Byte = ( 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^2+x^11 }
PP12 : Array[0..12] of Byte = ( 1 , 1 , 0 , 0 , 1 , 0 , 1 , 0 , 0 ,
0 , 0 , 0 , 1 ) ; { 1+x+x^4+x^6+x^12 }
PP13 : Array[0..13] of Byte = ( 1 , 1 , 0 , 1 , 1 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 1 ) ; { 1+x+x^3+x^4+x^13 }
PP14 : Array[0..14] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 1 , 0 , 0 ,
0 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x+x^6+x^10+x^14 }
PP15 : Array[0..15] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x+x^15 }
PP16 : Array[0..16] of Byte = ( 1 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x+x^3+x^12+x^16 }
constructor TiPentecGalois.Create();
begin
inherited;
pPP[2] := @PP2 ;
pPP[3] := @PP3 ;
pPP[4] := @PP4 ;
pPP[5] := @PP5 ;
pPP[6] := @PP6 ;
pPP[7] := @PP7 ;
pPP[8] := @PP8 ;
pPP[9] := @PP9 ;
pPP[10] := @PP10 ;
pPP[11] := @PP11 ;
pPP[12] := @PP12 ;
pPP[13] := @PP13 ;
pPP[14] := @PP14 ;
pPP[15] := @PP15 ;
pPP[16] := @PP16 ;
end;
procedure TiPentecGalois.InitTables(Dimension:integer);
begin
//* initialize the table of powers of alpha */
//init_exp_table();
init_exp_table(Dimension);
end;
procedure TiPentecGalois.init_exp_table(Dimension:integer);
var
i : Integer ;
mask : smallint ;
begin
FTableLength := (1 shl Dimension);
SetLength(gexp,FTableLength);
SetLength(glog,FTableLength);
SetPrimitive(PP , Dimension) ;
mask := 1 ;
gexp[Dimension] := 0 ;
for i := 0 to (Dimension - 1) do begin
gexp[i] := mask ;
glog[gexp[i]] := i ;
If pp[i] <> 0 then gexp[Dimension] := gexp[Dimension] xor mask ;
mask := mask shl 1 ;
end ;
glog[gexp[Dimension]] := Dimension ;
mask := mask shr 1 ;
for i := (Dimension + 1) to (FTableLength - 1) do Begin
if gexp[i-1] >= mask
then gexp[i] := gexp[Dimension] xor ((gexp[i-1] xor mask) shl 1)
else gexp[i] := gexp[i-1] shl 1 ;
glog[gexp[i]] := i ;
end ;
end;
procedure TiPentecGalois.SetPrimitive(var PP; nIdx : Integer ) ;
begin
Move(pPP[nIdx]^ , PP , (nIdx + 1)) ;
end ;
(*
// 底が8のみのバージョン //かつ gexpテーブル長が2倍になる
procedure TiPentecGalois.init_exp_table();
var
i,z:integer;
pinit,p1,p2,p3,p4,p5,p6,p7,p8:integer;
begin
pinit:=0; p2:=0; p3:=0; p4:=0; p5:=0; p6:=0; p7:=0; p8:=0;
p1:=1;
gexp[0] := 1;
gexp[255] := gexp[0];
glog[0] := 0; //* shouldn't log[0] be an error? */
for i := 1 to 256-1 do begin
pinit := p8;
p8 := p7;
p7 := p6;
p6 := p5;
p5 := p4 xor pinit;
p4 := p3 xor pinit;
p3 := p2 xor pinit;
p2 := p1;
p1 := pinit;
gexp[i] := p1 + p2*2 + p3*4 + p4*8 + p5*16 + p6*32 + p7*64 + p8*128;
gexp[i+255] := gexp[i];
end;
for i := 1 to 256-1 do begin
for z := 0 to 256-1 do begin
if (gexp[z] = i) then begin
glog[i] := z;
break;
end;
end;
end;
end;
*)
// /* multiplication using logarithms */
function TiPentecGalois.gmult(a:integer; b:integer):integer;
var
i,j:integer;
begin
if (a=0) or (b=0) then begin
result:=0;
end else begin
i := glog[a];
j := glog[b];
//result:=gexp[i+j]; //テーブルの長さが２倍のとき
result:=gexp[(i+j) mod (FTableLength-1)]; //テーブルの長さがぴったりのとき
end;
end;
function TiPentecGalois.ginv (elt:integer):integer;
begin
result:= gexp[255-glog[elt]];
end;
end.
Reed Solomon符号化実装ユニット
unit iPentecReedSolomon;
interface
uses
Windows,SysUtils,Classes;
Const
DEFAULT_MM = 8;
DEFAULT_NP = 32;
type
TiPentecReedSolomon = class(TComponent)
private
FGFDimension:integer;
FRSLength:integer;
FDataLength:integer;
FParityLength:integer;
FMaxRSLength:integer;
PP : Array[0..16] of Byte ; { specify irreducible polynomial coeffts }//ガロア体原始多項式
ExpToVector : Array of smallint; //alpha to
VectorToExp : Array of smallint; //index of
gg : Array of smallint ; // 誤り訂正生成多項式
bb : Array of smallint ;
data : Array of smallint ;
recd : Array of smallint ;
CodeWord : Array of Byte ;
procedure InitParameters(GFDimension:integer;ParityLength:integer);
procedure ReInitParameters(GFDimension:integer;ParityLength:integer;RSLength:integer);
procedure SetGFDimension(value:integer);
procedure SetParityLength(value:integer);
procedure SetRSLength(value:integer);
{ Private 宣言 }
protected
procedure generate_gf(mm:integer) ;
procedure gen_poly ;
procedure SetPrimitive(Var PP; nIdx : Integer ) ;
{ Protected 宣言 }
public
constructor Create(AOnwer:TCOmponent);override;
destructor Destroy;override;
procedure InitBuffers ;
procedure EncodeRS(var xData;var xEncoded);
function DecodeRS(var xData;var xDecoded):Integer ;
procedure EncodeRSEx(var xData:array of byte;var xEncoded:array of byte;var xParity:array of byte);
function DecodeRSEx(var xData:array of byte; var xParity:array of byte; var xDecoded:array of byte):Integer;
procedure SetRSMaxLength();
property DataLength:integer read FDataLength;
{ Public 宣言 }
published
property GFDimension:integer read FGFDimension Write SetGFDimension;
property ParityLength:integer read FParityLength write SetParityLength;
property RSLength:integer read FRSLength write SetRSLength;
{ Published 宣言 }
end;
procedure Register;
implementation
var
pPP : Array[2..16] of Pointer ;
const
//ガロア体原始多項式
PP2 : Array[0..2] of Byte = ( 1 , 1 , 1 ) ;
PP3 : Array[0..3] of Byte = ( 1 , 1 , 0 , 1 ); { 1 + x + x^3 }
PP4 : Array[0..4] of Byte = ( 1 , 1 , 0 , 0 , 1 ) ; { 1 + x + x^4 }
PP5 : Array[0..5] of Byte = ( 1 , 0 , 1 , 0 , 0 , 1 ) ; { 1 + x^2 + x^5 }
PP6 : Array[0..6] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 1 ) ; { 1 + x + x^6 }
PP7 : Array[0..7] of Byte = ( 1, 0, 0, 1, 0, 0, 0, 1 ) ; { 1 + x^3 + x^7 }
PP8 : Array[0..8] of Byte = ( 1 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x^2+x^3+x^4+x^8 }
PP9 : Array[0..9] of Byte = ( 1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^4+x^9 }
PP10 : Array[0..10] of Byte = ( 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^3+x^10 }
PP11 : Array[0..11] of Byte = ( 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x^2+x^11 }
PP12 : Array[0..12] of Byte = ( 1 , 1 , 0 , 0 , 1 , 0 , 1 , 0 , 0 ,
0 , 0 , 0 , 1 ) ; { 1+x+x^4+x^6+x^12 }
PP13 : Array[0..13] of Byte = ( 1 , 1 , 0 , 1 , 1 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 1 ) ; { 1+x+x^3+x^4+x^13 }
PP14 : Array[0..14] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 1 , 0 , 0 ,
0 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x+x^6+x^10+x^14 }
PP15 : Array[0..15] of Byte = ( 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 1 ) ; { 1+x+x^15 }
PP16 : Array[0..16] of Byte = ( 1 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 ) ; { 1+x+x^3+x^12+x^16 }
procedure Register;
begin
RegisterComponents('ipentec Component', [TiPentecReedSolomon]);
end;
constructor TiPentecReedSolomon.Create(AOnwer:TCOmponent);
begin
inherited;
InitBuffers ; { Clear the internal buffers }
{ generate the Galois Field GF(2**mm) }
//Generate_gf(DEFAULT_MM);
//FGFDimension:=DEFAULT_MM;
InitParameters(DEFAULT_MM,DEFAULT_NP);
end;
destructor TiPentecReedSolomon.Destroy ;
begin
Finalize(gg);
Finalize(bb);
Finalize(data);
Finalize(ExpToVector);
Finalize(VectorToExp);
Finalize(CodeWord);
Finalize(recd);
inherited;
end;
procedure TiPentecReedSolomon.InitParameters(GFDimension:integer;ParityLength:integer);
begin
FGFDimension:=GFDimension;
FParityLength:=ParityLength;
FMaxRSLength:=(1 shl FGFDimension) - 1; {nn=2**mm -1 Max length of codeword}
FRSLength:=(1 shl FGFDimension) - 1; {nn=2**mm -1 Max length of codeword}
FDataLength := FRSLength - FParityLength ; { data symbols, kk = nn-2*MaxErrors }
SetLength(gg,FParityLength+1);
SetLength(bb,FParityLength-1+1);
SetLength(data,FDataLength+1);
SetLength(ExpToVector,FRSLength+1);
SetLength(VectorToExp,FRSLength+1);
SetLength(CodeWord,FRSLength-1+1);
SetLength(recd,FRSLength-1+1);
generate_gf(FGFDimension);
gen_poly;
end;
procedure TiPentecReedSolomon.ReInitParameters(GFDimension:integer;ParityLength:integer;RSLength:integer);
begin
FGFDimension:=GFDimension;
FParityLength:=ParityLength;
FMaxRSLength:=(1 shl FGFDimension) - 1; {nn=2**mm -1 Max length of codeword}
FRSLength:=RSLength;
FDataLength := FRSLength - FParityLength ; { data symbols, kk = nn-2*MaxErrors }
SetLength(gg,FParityLength+1);
SetLength(bb,FParityLength-1+1);
SetLength(data,FDataLength+1);
SetLength(ExpToVector,FMaxRSLength+1);
SetLength(VectorToExp,FMaxRSLength+1);
SetLength(CodeWord,FRSLength-1+1);
SetLength(recd,FRSLength-1+1);
generate_gf(FGFDimension);
gen_poly;
end;
{ generate GF(2**mm) from the irreducible polynomial p(X)
in pp[0]..pp[mm]
lookup tables:
index->polynomial form alpha_to[] contains j=alpha**i ;
polynomial form -> index form index_of[j=alpha**i] = i
alpha = 2 is the primitive element of GF(2**mm)
}
{
ベクトル表現・べき表現変換テーブルの作成
}
Procedure TiPentecReedSolomon.Generate_GF(mm:integer) ;
Var
i : Integer ;
mask : smallint ;
Begin
SetPrimitive(PP , mm) ;
mask := 1 ;
ExpToVector[mm] := 0 ;
For i := 0 to (mm - 1) do Begin
ExpToVector[i] := mask ;
VectorToExp[ExpToVector[i]] := i ;
If pp[i] <> 0 then ExpToVector[mm] := ExpToVector[mm] xor mask ;
mask := mask shl 1 ;
End ;
VectorToExp[ExpToVector[mm]] := mm ;
mask := mask shr 1 ;
For i := (mm + 1) to (FRSLength - 1) do Begin
If ExpToVector[i-1] >= mask
then ExpToVector[i] := ExpToVector[mm] xor ((ExpToVector[i-1] xor mask) shl 1)
Else ExpToVector[i] := ExpToVector[i-1] shl 1 ;
VectorToExp[ExpToVector[i]] := i ;
End ;
VectorToExp[0] := -1 ;
End ;
{ Obtain the generator polynomial of the tt-error correcting, length
nn=(2**mm -1) Reed Solomon code from the product of
(X+alpha**i), i=1..2*tt
}
Procedure TiPentecReedSolomon.gen_poly ;
var
i, j : smallint;
begin
gg[0] := 2 ; { primitive element alpha = 2 for GF(2**mm) }
gg[1] := 1 ; { g(x) = (X+alpha) initially }
for i := 2 to (FRSLength - FDataLength) do begin
gg[i] := 1 ;
for j := (i - 1) downto 1 do begin
if (gg[j] <> 0)
then gg[j] := gg[j-1] xor ExpToVector[(VectorToExp[gg[j]]+i) mod FMaxRSLength]
else gg[j] := gg[j-1] ;
end;
{ gg[0] can never be zero }
gg[0] := ExpToVector[(VectorToExp[gg[0]] + i) mod FMaxRSLength] ;
end ;
{ Convert gg[] to index form for quicker encoding. }
for i := 0 to FParityLength do gg[i] := VectorToExp[gg[i]] ;
end ;
procedure TiPentecReedSolomon.SetPrimitive(Var PP; nIdx : Integer ) ;
begin
Move(pPP[nIdx]^ , PP , (nIdx + 1)) ;
end ;
{ TReedSolomon.InitBuffers }
procedure TiPentecReedSolomon.InitBuffers ;
begin
FillChar(data , SizeOf(data) , 0) ;
FillChar(recd , SizeOf(recd) , 0) ;
FillChar(CodeWord , SizeOf(CodeWord) , 0) ;
end ;
{ take the string of symbols in data[i], i=0..(k-1) and encode
systematically to produce 2*tt parity symbols in bb[0]..bb[2*tt-1]
data[] is input and bb[] is output in polynomial form. Encoding is
done by using a feedback shift register with appropriate connections
specified by the elements of gg[], which was generated above.
Codeword is c(X) = data(X)*X**(np)+ b(X) }
procedure TiPentecReedSolomon.EncodeRS(var xData; var xEncoded) ;
type
bArrType = Array[0..16383] of Byte ;
var
e:integer;
rswork: array[0..123 - 1] of Byte; // RS符号計算用作業領域
nI, i, j:Integer;
feedback:smallint;
axData:bArrType absolute xData ;
begin
for i:=0 to Length(bb)-1 do bb[i]:=0;
for nI := 0 to (FRSLength - 1) do Data[nI] := axData[nI] ;
for i := (FDataLength - 1) downto 0 do begin
feedback := VectorToExp[data[i] xor bb[FParityLength - 1]] ;
if (feedback <> -1) then begin
for j := (FParityLength - 1) downto 1 do begin
if (gg[j] <> -1)
then bb[j] := bb[j-1] xor ExpToVector[(gg[j] + feedback) mod FMaxRSLength]
else bb[j] := bb[j-1] ;
end;
bb[0] := ExpToVector[(gg[0] + feedback) mod FMaxRSLength] ;
end else begin
for j := (FParityLength - 1) downto 1 do bb[j] := bb[j-1] ;
bb[0] := 0 ;
end ;
end ;
{ put the transmitted codeword, made up of data }
{ plus parity, in CodeWord }
for nI := 0 to (FParityLength - 1) do Recd[nI] := bb[nI] ;
for nI := 0 to (FDataLength - 1) do Recd[nI + FParityLength] := Data[nI] ;
for nI := 0 to (FRSLength - 1) do CodeWord[nI] := Recd[nI] ;
Move(CodeWord[0] , xEncoded , FRSLength) ;
end ;
procedure TiPentecReedSolomon.EncodeRSEx(var xData:array of byte;var xEncoded:array of byte;var xParity:array of byte);
var
nI, i, j:Integer;
feedback:smallint;
e:integer;
invertgg:array of smallint;
begin
SetLength(invertgg,Length(gg));
for i:=0 to Length(gg)-1 do begin
invertgg[i]:=gg[Length(gg)-1-1-i];
end;
for i:=0 to Length(bb)-1 do bb[i]:=0;
for nI := 0 to (FRSLength - 1) do data[nI] := xData[nI] ;
for i := (FDataLength - 1) downto 0 do begin
feedback := VectorToExp[data[i] xor bb[FParityLength - 1]] ;
if (feedback <> -1) then begin
for j := (FParityLength - 1) downto 1 do begin
if (gg[j] <> -1)
then bb[j] := bb[j-1] xor ExpToVector[(gg[j] + feedback) mod FMaxRSLength]
else bb[j] := bb[j-1] ;
end;
bb[0] := ExpToVector[(gg[0] + feedback) mod FMaxRSLength] ;
end else begin
for j := (FParityLength - 1) downto 1 do bb[j] := bb[j-1] ;
bb[0] := 0 ;
end ;
end ;
{ put the transmitted codeword, made up of data }
{ plus parity, in CodeWord }
for nI := 0 to (FParityLength - 1) do begin
Recd[nI] := bb[nI] ;
xParity[ni] := bb[ni];
end;
for nI := 0 to (FDataLength - 1) do Recd[nI + FParityLength] := Data[nI] ;
for nI := 0 to (FRSLength - 1) do CodeWord[nI] := Recd[nI] ;
Move(CodeWord[0] , xEncoded , FRSLength) ;
end ;
{ assume we have received bits grouped into mm-bit symbols in recd[i],
i=0..(nn-1), and recd[i] is index form (ie as powers of alpha).
We first compute the 2*tt syndromes by substituting alpha**i into
rec(X) and evaluating, storing the syndromes in s[i], i=1..2tt
(leave s[0] zero). Then we use the Berlekamp iteration to find the
error location polynomial elp[i]. If the degree of the elp is >tt,
we cannot correct all the errors and hence just put out the information
symbols uncorrected. If the degree of elp is <=tt, we substitute
alpha**i , i=1..n into the elp to get the roots, hence the inverse
roots, the error location numbers. If the number of errors located
does not equal the degree of the elp, we have more than tt errors and
cannot correct them. Otherwise, we then solve for the error value at
the error location and correct the error. The procedure is that found
in Lin and Costello. For the cases where the number of errors is known
to be too large to correct, the information symbols as received are
output (the advantage of systematic encoding is that hopefully some
of the information symbols will be okay and that if we are in luck, the
errors are in the parity part of the transmitted codeword). Of course,
these insoluble cases can be returned as error flags to the calling
routine if desired. }
function TiPentecReedSolomon.DecodeRS(var xData; var xDecoded):Integer;
type
bArrType = Array[0..16383] of Byte ;
var
axData : bArrType absolute xData ;
axDecoded : bArrType absolute xDecoded ;
cStr : String ;
nI , nJ , nK : Integer ;
i , j : Integer ;
u , q : smallint ;
elp : array of array of smallint; //Array[0..np+1 , 0..np - 1] of smallint ;
d : array of smallint;//Array[0..np+1] of smallint ;
l : array of smallint;//Array[0..np+1] of smallint ;
u_lu : array of smallint;//Array[0..np+1] of smallint ;
s : array of smallint;//Array[0..np] of smallint ;
count : smallint ;
syn_error : smallint ;
lpc:integer;
var
root : array of smallint;//Array[0..MaxErrors - 1] of smallint ;
loc : array of smallint;//Array[0..MaxErrors - 1] of smallint ;
z : array of smallint;//Array[0..MaxErrors] of smallint ;
err : array of smallint;//Array[0..nn-1] of smallint ;
reg : array of smallint;//Array[0..MaxErrors] of smallint ;
begin { DecodeRS }
SetLength(elp,FParityLength+1+1);
for lpc:=0 to Length(elp)-1 do begin
SetLength(elp[lpc],FParityLength-1+1);
end;
SetLength(d,FParityLength+1+1);
SetLength(l,FParityLength+1+1);
SetLength(u_lu,FParityLength+1+1);
SetLength(s,FParityLength+1);
//
SetLength(root,Trunc(FParityLength/2)-1+1);
SetLength(loc,Trunc(FParityLength/2)-1+1);
SetLength(z,Trunc(FParityLength/2)+1);
SetLength(err,FRSLength-1+1);
SetLength(reg,Trunc(FParityLength/2)+1);
for nI := 0 to (FRSLength - 1) do Recd[nI] := axData[nI] ;
for i := 0 to (FRSLength - 1) do recd[i] := VectorToExp[recd[i]] ; { put recd[i] into index form }
count := 0 ;
syn_error := 0 ;
Result := 0 ;
{ first form the syndromes }
for i := 1 to FParityLength do begin
s[i] := 0 ;
for j := 0 to (FRSLength - 1) do begin
if recd[j] <> -1 then begin { recd[j] in index form }
s[i] := s[i] xor ExpToVector[(recd[j] + i * j) mod FMaxRSLength] ;
end ;
end;
{ convert syndrome from polynomial form to index form }
if (s[i] <> 0) then begin
syn_error := 1 ; { set flag if non-zero syndrome => error }
end ;
s[i] := VectorToExp[s[i]] ;
end ;
if syn_error <> 0 then begin { if errors, try and correct }
{ Compute the error location polynomial via the Berlekamp }
{ iterative algorithm, following the terminology of Lin and }
{ Costello : d[u] is the 'mu'th discrepancy, where u='mu'+1 }
{ and 'mu' (the Greek letter!) is the step number ranging from }
{ -1 to 2*tt (see L&C), l[u] is the degree of the elp at that }
{ step, and u_l[u] is the difference between the step number }
{and the degree of the elp. }
{ Initialize table entries }
d[0] := 0 ; { index form }
d[1] := s[1] ; { index form }
elp[0][0] := 0 ; { index form }
elp[1][0] := 1 ; { polynomial form }
for i := 1 to (FParityLength - 1) do begin
elp[0][i] := -1 ; { index form }
elp[1][i] := 0 ; { polynomial form }
end ;
l[0] := 0 ;
l[1] := 0 ;
u_lu[0] := -1 ;
u_lu[1] := 0 ;
u := 0 ;
while ((u < FParityLength) and (l[u+1] <= Trunc(FParityLength/2))) do begin
Inc(u) ;
if (d[u] = -1) then begin
l[u+1] := l[u] ;
for i := 0 to l[u] do begin
elp[u+1][i] := elp[u][i] ;
elp[u][i] := vectorToExp[elp[u][i]] ;
end ;
end else begin
{ search for words with greatest u_lu[q] for which d[q]!=0 }
q := u - 1 ;
while ((d[q] = -1) and (q > 0)) do dec(q) ;
{ have found first non-zero d[q] }
if (q > 0) then begin
j := q ;
while j > 0 do begin
Dec(j) ;
If ((d[j] <> -1) and (u_lu[q] < u_lu[j])) then q := j ;
end ;
end ;
{ have now found q such that d[u]!=0 and u_lu[q] is maximum }
{ store degree of new elp polynomial }
if (l[u] > l[q]+u-q) then l[u+1] := l[u] else l[u+1] := l[q] + u - q ;
{ form new elp(x) }
for i := 0 to (FParityLength - 1) do elp[u+1][i] := 0 ;
for i := 0 to l[q] do begin
//if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
end;
for i := 0 to l[u] do begin
elp[u+1][i] := elp[u+1][i] xor elp[u][i] ;
{ convert old elp value to index }
elp[u][i] := VectorToExp[elp[u][i]] ;
end ;
end ;
u_lu[u+1] := u - l[u+1] ;
{ form (u+1)th discrepancy }
if u < FParityLength then begin { no discrepancy computed on last iteration }
if (s[u+1]<> -1) then d[u+1] := ExpToVector[s[u+1]] else d[u+1] := 0 ;
for i := 1 to l[u+1] do begin
if ((s[u+1-i] <> -1) and (elp[u+1][i] <> 0))
then d[u+1] := d[u+1] xor ExpToVector[(s[u+1-i]+VectorToExp[elp[u+1][i]]) mod FMaxRSLength] ;
end;
{ put d[u+1] into index form }
d[u+1] := VectorToExp[d[u+1]] ;
end ;
end ; { end While }
Inc(u) ;
if l[u] <= Trunc(FParityLength/2) then begin { can correct error }
{ put elp into index form }
for i := 0 to l[u] do elp[u][i] := VectorToExp[elp[u][i]] ;
{ find roots of the error location polynomial }
for i := 1 to l[u] do begin
reg[i] := elp[u][i] ;
end ;
//for i := 1 to FRSLength do begin
for i := 1 to FMaxRSLength do begin
q := 1 ;
for j := 1 to l[u] do begin
if reg[j] <> -1 then begin
reg[j] := (reg[j] + j) mod FMaxRSLength ;
q := q xor ExpToVector[reg[j]] ;
end ;
end;
if q = 0 then begin { store root and error location number indices }
root[count] := i ;
//loc[count] := FRSLength - i ;
loc[count] := FMaxRSLength - i ;
Inc(count) ;
end ;
end ;
if count = l[u] then begin { no. roots = degree of elp hence <= tt errors }
Result := count ;
{ form polynomial z(x) }
for i := 1 to l[u] do begin { Z[0] = 1 always - do not need }
if ((s[i]<> -1) and (elp[u][i]<> -1)) then z[i] := ExpToVector[s[i]] xor ExpToVector[elp[u][i]]
else if ((s[i] <> -1) and (elp[u][i] = -1)) then z[i] := ExpToVector[s[i]]
else if ((s[i] = -1) and (elp[u][i] <> -1)) then z[i] := ExpToVector[elp[u][i]]
else z[i] := 0 ;
for j := 1 to (i - 1) do begin
if ((s[j] <> -1) and (elp[u][i-j] <> -1)) then z[i] := z[i] xor ExpToVector[(elp[u][i-j] + s[j]) mod FMaxRSLength ] ;
end;
{ put into index form }
z[i] := VectorToExp[z[i]] ;
end ;
{ evaluate errors at locations given by }
{ error location numbers loc[i] }
for i := 0 to (FRSLength - 1) do begin
//for i := 0 to (FMaxRSLength - 1) do begin
err[i] := 0 ;
{ convert recd[] to polynomial form }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] else recd[i] := 0 ;
end ;
for i := 0 to (l[u] - 1) do begin { compute numerator of error term first }
err[loc[i]] := 1 ; { accounts for z[0] }
for j := 1 to l[u] do begin
if z[j] <> -1 then err[loc[i]] := err[loc[i]] xor ExpToVector[(z[j]+j * root[i]) mod FMaxRSLength] ;
end;
if err[loc[i]] <> 0 then begin
err[loc[i]] := VectorToExp[err[loc[i]]] ;
q := 0 ; { form denominator of error term }
for j := 0 to (l[u] - 1) do begin
if j <> i then q := q + VectorToExp[1 xor ExpToVector[(loc[j] + root[i]) mod FMaxRSLength]] ;
end;
q := q mod FMaxRSLength ;
//err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] := ExpToVector[(err[loc[i]] - q + FMaxRSLength) mod FMaxRSLength] ;
{ recd[i] must be in polynomial form }
recd[loc[i]] := recd[loc[i]] xor err[loc[i]] ;
end ;
end ;
end else begin
{ no. roots != degree of elp => >tt errors and cannot solve }
Result := -1 ; { Signal an error. }
for i := 0 to (FRSLength - 1) do begin { could return error flag if desired }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ; { just output received codeword as is }
end;
end ;
end else begin { if l[u] <= tt then }
{ elp has degree has degree >tt hence cannot solve }
for i := 0 to (FRSLength - 1) do begin { could return error flag if desired }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ; { just output received codeword as is }
end;
end;
end else begin; { If syn_error <> 0 then }
{ no non-zero syndromes => no errors: output received codeword }
for i := 0 to (FRSLength - 1) do begin
if recd[i] <> -1
then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ;
Result := 0 ; { No errors ocurred. }
end ;
end;
for nI := 0 to (FRSLength - 1) do axDecoded[nI] := Recd[nI] ;
end;
{ TReedSolomon.DecodeRS }
function TiPentecReedSolomon.DecodeRSEx(var xData:array of byte; var xParity:array of byte;
var xDecoded:array of byte):Integer;
var
cStr : String ;
nI , nJ , nK : Integer ;
i , j : Integer ;
u , q : smallint ;
elp : array of array of smallint; //Array[0..np+1 , 0..np - 1] of smallint ;
d : array of smallint;//Array[0..np+1] of smallint ;
l : array of smallint;//Array[0..np+1] of smallint ;
u_lu : array of smallint;//Array[0..np+1] of smallint ;
s : array of smallint;//Array[0..np] of smallint ;
count : smallint ;
syn_error : smallint ;
lpc:integer;
var
root : array of smallint;//Array[0..MaxErrors - 1] of smallint ;
loc : array of smallint;//Array[0..MaxErrors - 1] of smallint ;
z : array of smallint;//Array[0..MaxErrors] of smallint ;
err : array of smallint;//Array[0..nn-1] of smallint ;
reg : array of smallint;//Array[0..MaxErrors] of smallint ;
begin { DecodeRS }
SetLength(elp,FParityLength+1+1);
for lpc:=0 to Length(elp)-1 do begin
SetLength(elp[lpc],FParityLength-1+1);
end;
SetLength(d,FParityLength+1+1);
SetLength(l,FParityLength+1+1);
SetLength(u_lu,FParityLength+1+1);
SetLength(s,FParityLength+1);
//
SetLength(root,Trunc(FParityLength/2)-1+1);
SetLength(loc,Trunc(FParityLength/2)-1+1);
SetLength(z,Trunc(FParityLength/2)+1);
SetLength(err,FRSLength-1+1);
SetLength(reg,Trunc(FParityLength/2)+1);
for nI := 0 to (FParityLength - 1) do Recd[nI] := xParity[nI] ;
for nI := 0 to (FRSLength - 1) do Recd[FParityLength + nI] := xData[nI] ;
for i := 0 to (FRSLength - 1) do Recd[i] := VectorToExp[recd[i]] ; { put recd[i] into index form }
count := 0 ;
syn_error := 0 ;
Result := 0 ;
{ first form the syndromes }
for i := 1 to FParityLength do begin
s[i] := 0 ;
for j := 0 to (FRSLength - 1) do begin
if recd[j] <> -1 then begin { recd[j] in index form }
s[i] := s[i] xor ExpToVector[(recd[j] + i * j) mod FMaxRSLength] ;
end ;
end;
{ convert syndrome from polynomial form to index form }
if (s[i] <> 0) then begin
syn_error := 1 ; { set flag if non-zero syndrome => error }
end ;
s[i] := VectorToExp[s[i]] ;
end ;
if syn_error <> 0 then begin { if errors, try and correct }
{ Compute the error location polynomial via the Berlekamp }
{ iterative algorithm, following the terminology of Lin and }
{ Costello : d[u] is the 'mu'th discrepancy, where u='mu'+1 }
{ and 'mu' (the Greek letter!) is the step number ranging from }
{ -1 to 2*tt (see L&C), l[u] is the degree of the elp at that }
{ step, and u_l[u] is the difference between the step number }
{and the degree of the elp. }
{ Initialize table entries }
d[0] := 0 ; { index form }
d[1] := s[1] ; { index form }
elp[0][0] := 0 ; { index form }
elp[1][0] := 1 ; { polynomial form }
for i := 1 to (FParityLength - 1) do begin
elp[0][i] := -1 ; { index form }
elp[1][i] := 0 ; { polynomial form }
end ;
l[0] := 0 ;
l[1] := 0 ;
u_lu[0] := -1 ;
u_lu[1] := 0 ;
u := 0 ;
while ((u < FParityLength) and (l[u+1] <= Trunc(FParityLength/2))) do begin
Inc(u) ;
if (d[u] = -1) then begin
l[u+1] := l[u] ;
for i := 0 to l[u] do begin
elp[u+1][i] := elp[u][i] ;
elp[u][i] := vectorToExp[elp[u][i]] ;
end ;
end else begin
{ search for words with greatest u_lu[q] for which d[q]!=0 }
q := u - 1 ;
while ((d[q] = -1) and (q > 0)) do dec(q) ;
{ have found first non-zero d[q] }
if (q > 0) then begin
j := q ;
while j > 0 do begin
Dec(j) ;
If ((d[j] <> -1) and (u_lu[q] < u_lu[j])) then q := j ;
end ;
end ;
{ have now found q such that d[u]!=0 and u_lu[q] is maximum }
{ store degree of new elp polynomial }
if (l[u] > l[q]+u-q) then l[u+1] := l[u] else l[u+1] := l[q] + u - q ;
{ form new elp(x) }
for i := 0 to (FParityLength - 1) do elp[u+1][i] := 0 ;
for i := 0 to l[q] do begin
//if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
if (elp[q][i] <> -1) then elp[u+1][i+u-q] := ExpToVector[(d[u] + FMaxRSLength - d[q] + elp[q][i]) mod FMaxRSLength] ;
end;
for i := 0 to l[u] do begin
elp[u+1][i] := elp[u+1][i] xor elp[u][i] ;
{ convert old elp value to index }
elp[u][i] := VectorToExp[elp[u][i]] ;
end ;
end ;
u_lu[u+1] := u - l[u+1] ;
{ form (u+1)th discrepancy }
if u < FParityLength then begin { no discrepancy computed on last iteration }
if (s[u+1]<> -1) then d[u+1] := ExpToVector[s[u+1]] else d[u+1] := 0 ;
for i := 1 to l[u+1] do begin
if ((s[u+1-i] <> -1) and (elp[u+1][i] <> 0))
then d[u+1] := d[u+1] xor ExpToVector[(s[u+1-i]+VectorToExp[elp[u+1][i]]) mod FMaxRSLength] ;
end;
{ put d[u+1] into index form }
d[u+1] := VectorToExp[d[u+1]] ;
end ;
end ; { end While }
Inc(u) ;
if l[u] <= Trunc(FParityLength/2) then begin { can correct error }
{ put elp into index form }
for i := 0 to l[u] do elp[u][i] := VectorToExp[elp[u][i]] ;
{ find roots of the error location polynomial }
for i := 1 to l[u] do begin
reg[i] := elp[u][i] ;
end ;
//for i := 1 to FRSLength do begin
for i := 1 to FMaxRSLength do begin
q := 1 ;
for j := 1 to l[u] do begin
if reg[j] <> -1 then begin
reg[j] := (reg[j] + j) mod FMaxRSLength ;
q := q xor ExpToVector[reg[j]] ;
end ;
end;
if q = 0 then begin { store root and error location number indices }
root[count] := i ;
//loc[count] := FRSLength - i ;
loc[count] := FMaxRSLength - i ;
Inc(count) ;
end ;
end ;
if count = l[u] then begin { no. roots = degree of elp hence <= tt errors }
Result := count ;
{ form polynomial z(x) }
for i := 1 to l[u] do begin { Z[0] = 1 always - do not need }
if ((s[i]<> -1) and (elp[u][i]<> -1)) then z[i] := ExpToVector[s[i]] xor ExpToVector[elp[u][i]]
else if ((s[i] <> -1) and (elp[u][i] = -1)) then z[i] := ExpToVector[s[i]]
else if ((s[i] = -1) and (elp[u][i] <> -1)) then z[i] := ExpToVector[elp[u][i]]
else z[i] := 0 ;
for j := 1 to (i - 1) do begin
if ((s[j] <> -1) and (elp[u][i-j] <> -1)) then z[i] := z[i] xor ExpToVector[(elp[u][i-j] + s[j]) mod FMaxRSLength ] ;
end;
{ put into index form }
z[i] := VectorToExp[z[i]] ;
end ;
{ evaluate errors at locations given by }
{ error location numbers loc[i] }
for i := 0 to (FRSLength - 1) do begin
//for i := 0 to (FMaxRSLength - 1) do begin
err[i] := 0 ;
{ convert recd[] to polynomial form }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] else recd[i] := 0 ;
end ;
for i := 0 to (l[u] - 1) do begin { compute numerator of error term first }
err[loc[i]] := 1 ; { accounts for z[0] }
for j := 1 to l[u] do begin
if z[j] <> -1 then err[loc[i]] := err[loc[i]] xor ExpToVector[(z[j]+j * root[i]) mod FMaxRSLength] ;
end;
if err[loc[i]] <> 0 then begin
err[loc[i]] := VectorToExp[err[loc[i]]] ;
q := 0 ; { form denominator of error term }
for j := 0 to (l[u] - 1) do begin
if j <> i then q := q + VectorToExp[1 xor ExpToVector[(loc[j] + root[i]) mod FMaxRSLength]] ;
end;
q := q mod FMaxRSLength ;
//err[loc[i]] := ExpToVector[(err[loc[i]] - q + FRSLength) mod FMaxRSLength] ;
err[loc[i]] := ExpToVector[(err[loc[i]] - q + FMaxRSLength) mod FMaxRSLength] ;
{ recd[i] must be in polynomial form }
recd[loc[i]] := recd[loc[i]] xor err[loc[i]] ;
end ;
end ;
end else begin
{ no. roots != degree of elp => >tt errors and cannot solve }
Result := -1 ; { Signal an error. }
for i := 0 to (FRSLength - 1) do begin { could return error flag if desired }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ; { just output received codeword as is }
end;
end ;
end else begin { if l[u] <= tt then }
{ elp has degree has degree >tt hence cannot solve }
for i := 0 to (FRSLength - 1) do begin { could return error flag if desired }
if recd[i] <> -1 then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ; { just output received codeword as is }
end;
end;
end else begin; { If syn_error <> 0 then }
{ no non-zero syndromes => no errors: output received codeword }
for i := 0 to (FRSLength - 1) do begin
if recd[i] <> -1
then recd[i] := ExpToVector[recd[i]] { convert recd[] to polynomial form }
else recd[i] := 0 ;
Result := 0 ; { No errors ocurred. }
end ;
end;
//for nI := 0 to (FRSLength - 1) do axDecoded[nI] := Recd[nI] ;
for nI := FParityLength to (FRSLength - 1) do xDecoded[nI-FParityLength] := Recd[nI] ;
end; { TReedSolomon.DecodeRS }
procedure TiPentecReedSolomon.SetGFDimension(value:integer);
begin
FGFDimension:=Value;
//ガロア体の次元が変わったらすべて作成しなおす必要がある
ReInitParameters(FGFDimension,FParityLength,FRSLength);
end;
procedure TiPentecReedSolomon.SetParityLength(value:integer);
begin
FParityLength:=Value;
FDataLength := FRSLength - FParityLength ; { data symbols, kk = nn-2*MaxErrors }
//パリティの長さが変わった場合は エラー訂正多項式を再作成
SetLength(gg,FParityLength+1);
SetLength(bb,FParityLength);
SetLength(data,FDataLength+1);
gen_poly;
end;
procedure TiPentecReedSolomon.SetRSLength(value:integer);
begin
FRSLength:=Value;
FDataLength := FRSLength - FParityLength ; { data symbols, kk = nn-2*MaxErrors }
//データの長さが変わった場合は エラー訂正多項式を再作成
SetLength(CodeWord,FRSLength-1+1);
SetLength(recd,FRSLength-1+1);
gen_poly;
end;
procedure TiPentecReedSolomon.SetRSMaxLength();
var
ml:integer;
begin
ml:=(1 shl FGFDimension) - 1 ; {nn=2**mm -1 length of codeword};
SetRSLength(ml);
end;
initialization
{ Initialize the Primitive Polynomial vector. }
pPP[2] := @PP2 ;
pPP[3] := @PP3 ;
pPP[4] := @PP4 ;
pPP[5] := @PP5 ;
pPP[6] := @PP6 ;
pPP[7] := @PP7 ;
pPP[8] := @PP8 ;
pPP[9] := @PP9 ;
pPP[10] := @PP10 ;
pPP[11] := @PP11 ;
pPP[12] := @PP12 ;
pPP[13] := @PP13 ;
pPP[14] := @PP14 ;
pPP[15] := @PP15 ;
pPP[16] := @PP16 ;
Randomize ;
end.
Reed Solomon符号化実装ユニット(別版)
unit iPentecReedSolomonCodec;
interface
uses
SysUtils, Classes,iPentecGalois;
type
TiPentecReedSolomonCodec = class(TComponent)
private
NPAR:integer;
MAXDEG:integer;
FDimension:integer;
GT:TiPentecGalois;
FNErasures:integer;
ErasureLocs:array [0..256-1] of integer;
ErrorLocs:array [0..256-1] of integer;
NErrors:integer;
Lambda:array of integer;//[MAXDEG];
Omega:array of integer;//[MAXDEG];
procedure SetDimension(value:integer);
{ Private 宣言 }
protected
procedure compute_genpoly (nbytes:integer; out genpoly:array of integer);
procedure build_codeword (msg:array of char; nbytes:integer; out dst:array of char);
procedure Modified_Berlekamp_Massey;
procedure init_gamma(out gamma:array of integer);
procedure mul_z_poly(var src:array of integer);
function compute_discrepancy(lambda:array of integer;S:array of integer;L:integer;n:integer):integer;
procedure compute_modified_omega();
procedure Find_Roots();
{ Protected 宣言 }
public
pBytes:array of integer;
synBytes:array of integer;
genPoly: array of integer; //[MAXDEG*2];
constructor Create(Aowner:TComponent);override;
destructor Destroy;override;
procedure initialize_ecc (ParityLength:integer);
procedure decode_data (var data: array of char; nbytes:integer);
procedure encode_data (msg:array of char; nbytes:integer; out dst:array of char);
procedure encode_data_parity(msg:array of char; nbytes:integer;out Parity: array of char);
function check_syndrome ():integer;
function crc_ccitt(msg:Pchar; leng:integer):smallint;
function correct_errors_erasures (var codeword:array of char; csize:integer; nerasures:integer;
erasures:array of integer):integer;
procedure add_polys(out dst:array of integer; src:array of integer);
procedure scale_poly(k:integer; out poly:array of integer);
procedure mult_polys(out dst:array of integer; p1:array of integer; p2:array of integer);
procedure copy_poly(out dst:array of integer; src:array of integer);
procedure zero_poly(out poly:array of integer);
{ Public 宣言 }
published
property Dimension:integer read FDimension write SetDimension;
{ Published 宣言 }
end;
procedure Register;
implementation
procedure Register;
begin
RegisterComponents('iPentec Component', [TiPentecReedSolomonCodec]);
end;
constructor TiPentecReedSolomonCodec.Create(Aowner:TComponent);
begin
inherited;
FDimension:=8;
GT:=TiPentecGalois.Create;
end;
destructor TiPentecReedSolomonCodec.destroy;
begin
GT.Free;
inherited;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
プロパティセット関数
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
procedure TiPentecReedSolomonCodec.SetDimension(value:integer);
begin
if FDimension <> value then begin
FDimension:=value;
end;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
初期化
}
procedure TiPentecReedSolomonCodec.initialize_ecc (ParityLength:integer);
begin
NPAR:=ParityLength;
MAXDEG := NPAR * 2;
//* Initialize the galois field arithmetic tables */
GT.InitTables(FDimension); //init_galois_tables();
//* Compute the encoder generator polynomial */
SetLength(genPoly,MAXDEG*2);
compute_genpoly(NPAR, genPoly);
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Decode
デコード関数
デコード関数を呼び出したのち check_syndrome でシンドロームを求める
}
procedure TiPentecReedSolomonCodec.decode_data (var data: array of char; nbytes:integer);
var
i,j,sum:integer;
begin
SetLength(synBytes,MAXDEG);
for j:=0 to NPAR-1 do begin
//for j:=0 to 8-1 do begin
sum := 0;
for i := 0 to nbytes-1 do begin
sum := Byte(data[i]) xor GT.gmult(GT.gexp[j+1], sum);
end;
synBytes[j] := sum;
end;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Check Syndrome
シンドローム算出
Decodeを事前に呼び出す必要がある
}
function TiPentecReedSolomonCodec.check_syndrome ():integer;
var
i,nz:integer;
begin
nz:=0;
for i:=0 to NPAR-1 do begin
if synBytes[i] <>0 then nz := 1;
end;
result:=nz;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Encode
エンコード関数 コードワード版
}
procedure TiPentecReedSolomonCodec.encode_data (msg:array of char; nbytes:integer;out dst: array of char);
var
i,dbyte,j:integer;
LFSR:array of integer;
begin
SetLength(LFSR,NPAR+1);
for i:=0 to NPAR do LFSR[i]:=0;
for i := 0 to nbytes-1 do begin
dbyte := Byte(msg[i]) xor LFSR[NPAR-1];
for j:=NPAR-1 downto 1 do begin
LFSR[j] := LFSR[j-1] xor GT.gmult(genPoly[j], dbyte);
end;
LFSR[0] := GT.gmult(genPoly[0], dbyte);
end;
SetLength(pBytes,NPAR);
for i := 0 to NPAR-1 do pBytes[i] := LFSR[i];
build_codeword(msg, nbytes, dst);
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Encode Parity
エンコード関数 パリティ分離版
}
procedure TiPentecReedSolomonCodec.encode_data_parity(msg:array of char; nbytes:integer;out Parity: array of char);
var
i,dbyte,j:integer;
LFSR:array of integer;
begin
SetLength(LFSR,NPAR+1);
for i:=0 to NPAR do LFSR[i]:=0;
for i := 0 to nbytes-1 do begin
dbyte := Byte(msg[i]) xor LFSR[NPAR-1];
for j:=NPAR-1 downto 1 do begin
LFSR[j] := LFSR[j-1] xor GT.gmult(genPoly[j], dbyte);
end;
LFSR[0] := GT.gmult(genPoly[0], dbyte);
end;
SetLength(pBytes,NPAR);
for i := 0 to NPAR-1 do Parity[i] := Char(LFSR[i]);
//for i := 0 to NPAR-1 do Pareity := pBytes[i];
end;
function TiPentecReedSolomonCodec.crc_ccitt(msg:Pchar; leng:integer):smallint;
begin
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
correct errors erasures
誤り訂正処理
}
function TiPentecReedSolomonCodec.correct_errors_erasures (var codeword:array of char;
csize:integer; nerasures:integer; erasures:array of integer):integer;
var
r,i,j,err:integer;
num,denom:integer;
begin
// If you want to take advantage of erasure correction, be sure to
// set NErasures and ErasureLocs[] with the locations of erasures.
FNErasures := nerasures;
for i:=0 to NErasures-1 do ErasureLocs[i] := erasures[i];
Modified_Berlekamp_Massey();
Find_Roots();
if (NErrors <= NPAR) and (NErrors > 0) then begin
// /* first check for illegal error locs */
for r:=0 to NErrors-1 do begin
if ErrorLocs[r] >= csize then begin
result:=0;
exit;
end;
end;
for r:=0 to NErrors-1 do begin
i := ErrorLocs[r];
// /* evaluate Omega at alpha^(-i) */
num := 0;
for j:=0 to MAXDEG-1 do num := num xor GT.gmult(Omega[j], GT.gexp[((255-i)*j) mod 255]);
// /* evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear */
denom := 0;
j:=1;
while (j<MAXDEG) do begin
denom := denom xor GT.gmult(Lambda[j], GT.gexp[((255-i)*(j-1)) mod 255]);
inc(j,2);
end;
err := GT.gmult(num, GT.ginv(denom));
codeword[csize-i-1] := Char( Byte(codeword[csize-i-1]) xor err );
end;
result:=1;
exit;
end else begin
result:=0;
end;
end;
//
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Compute GenPoly
誤り訂正関数作成
誤り訂正関数はエンコード時のみ使用される
}
procedure TiPentecReedSolomonCodec.compute_genpoly (nbytes:integer; out genpoly:array of integer);
var
i:integer;
tp:array [0..255] of integer;
tp1:array [0..255] of integer;
begin
// /* multiply (x + a^n) for n = 1 to nbytes */
zero_poly(tp1);
tp1[0] := 1;
for i := 1 to nbytes do begin
zero_poly(tp);
tp[0] := GT.gexp[i]; //* set up x+a^n */
tp[1] := 1;
mult_polys(genpoly, tp, tp1);
copy_poly(tp1, genpoly);
end;
//QRコード用 誤り訂正多項式
//for i := 0 to nbytes do begin
// genpoly[i]:=GT.gexp[GT.glog[genpoly[i]] - (nbytes-i)];
//end;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Build Codeword
データの末尾にパリティを加える処理関数
}
procedure TiPentecReedSolomonCodec.build_codeword (msg:array of char; nbytes:integer;
out dst:array of char);
var
i:integer;
begin
for i := 0 to nbytes-1 do dst[i] := msg[i];
for i := 0 to NPAR-1 do dst[i+nbytes] := Char(pBytes[NPAR-1-i]);
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
エラー訂正関数
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Compute OMEGA
}
procedure TiPentecReedSolomonCodec.Modified_Berlekamp_Massey;
var
n,l,l2,k,d,i:integer;
psi:array of integer;
psi2:array of integer;
DA:array of integer;
gamma:array of integer;
begin
SetLength(psi,MAXDEG);
SetLength(psi2,MAXDEG);
SetLength(DA,MAXDEG);
SetLength(gamma,MAXDEG);
// /* initialize Gamma, the erasure locator polynomial */
init_gamma(gamma);
// /* initialize to z */
copy_poly(DA, gamma);
mul_z_poly(DA);
copy_poly(psi, gamma);
k := -1;
L := FNErasures;
for n:=FNErasures to NPAR-1 do begin
//for n:=FNErasures to 8-1 do begin
d := compute_discrepancy(psi, synBytes, L, n);
if d <> 0 then begin
// /* psi2 = psi - d*D */
for i:= 0 to MAXDEG-1 do psi2[i] := psi[i] xor GT.gmult(d, DA[i]);
if L < (n-k) then begin
L2 := n-k;
k := n-L;
// /* D = scale_poly(ginv(d), psi); */
for i:=0 to MAXDEG-1 do DA[i] := GT.gmult(psi[i], GT.ginv(d));
L := L2;
end;
// /* psi = psi2 */
for i:=0 to MAXDEG-1 do psi[i] := psi2[i];
end;
mul_z_poly(DA);
end;
//λの初期化
SetLength(Lambda,MAXDEG);
for i:=0 to MAXDEG-1 do Lambda[i] := psi[i];
compute_modified_omega();
end;
//* gamma = product (1-z*a^Ij) for erasure locs Ij */
procedure TiPentecReedSolomonCodec.init_gamma (out gamma:array of integer);
var
e:integer;
tmp:array of integer;
begin
SetLength(tmp,MAXDEG);
zero_poly(gamma);
zero_poly(tmp);
gamma[0] := 1;
for e:=0 to FNErasures-1 do begin
copy_poly(tmp, gamma);
scale_poly(GT.gexp[ErasureLocs[e]], tmp);
mul_z_poly(tmp);
add_polys(gamma, tmp);
end;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Compute Discrepancy
}
function TiPentecReedSolomonCodec.compute_discrepancy(lambda:array of integer;S:array of integer;L:integer;n:integer):integer;
var
i,sum:integer;
begin
sum:=0;
for i:=0 to L do sum := sum xor GT.gmult(lambda[i], S[n-i]);
result:=sum;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Compute OMEGA
}
// given Psi (called Lambda in Modified_Berlekamp_Massey) and synBytes,
// compute the combined erasure/error evaluator polynomial as
// Psi*S mod z^4
procedure TiPentecReedSolomonCodec.compute_modified_omega();
var
i:integer;
product:array of integer;
begin
SetLength(product, MAXDEG*2);
mult_polys(product, Lambda, synBytes);
//ωの初期化 (0でクリア後 productで初期値を代入)
SetLength(Omega,MAXDEG);
zero_poly(Omega);
for i:=0 to NPAR-1 do Omega[i] := product[i];
Finalize(product);
end;
//_/
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
Find Root
}
procedure TiPentecReedSolomonCodec.Find_Roots();
var
sum,r,k:integer;
begin
NErrors := 0;
for r:=1 to 256-1 do begin
sum := 0;
// /* evaluate lambda at r */
for k:=0 to NPAR+1-1 do begin
sum := sum xor GT.gmult(GT.gexp[(k*r) mod 255], Lambda[k]);
end;
if (sum = 0) then begin
ErrorLocs[NErrors] := 255-r;
inc(NErrors);
end;
end;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
配列処理関数
}
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
0による初期化
}
procedure TiPentecReedSolomonCodec.zero_poly(out poly:array of integer);
var
i:integer;
begin
for i := 0 to MAXDEG-1 do poly[i] := 0;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
コピー
}
procedure TiPentecReedSolomonCodec.copy_poly(out dst:array of integer;src:array of integer);
var
i:integer;
begin
for i := 0 to MAXDEG-1 do dst[i] := src[i];
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
加算
}
procedure TiPentecReedSolomonCodec.add_polys(out dst:array of integer; src:array of integer);
var
i:integer;
begin
for i := 0 to MAXDEG-1 do dst[i] := dst[i] xor src[i];
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
シフト積
}
// /* multiply by z, i.e., shift right by 1 */
procedure TiPentecReedSolomonCodec.mul_z_poly (var src:array of integer);
var
i:integer;
begin
for i:=MAXDEG-1 downto 1 do src[i] := src[i-1];
src[0] := 0;
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
定数との積
}
procedure TiPentecReedSolomonCodec.scale_poly(k:integer; out poly:array of integer);
var
i:integer;
begin
for i := 0 to MAXDEG-1 do poly[i] := GT.gmult(k, poly[i]);
end;
//_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
{
多項式との積
}
procedure TiPentecReedSolomonCodec.mult_polys(out dst:array of integer; p1:array of integer; p2:array of integer);
var
i,j:integer;
tmp1:array of integer;
begin
SetLength(tmp1,MAXDEG*2);
for i:=0 to (MAXDEG*2)-1 do dst[i] := 0;
for i:=0 to MAXDEG-1 do begin
for j:=MAXDEG to (MAXDEG*2)-1 do tmp1[j]:=0;
// /* scale tmp1 by p1[i] */
for j:=0 to MAXDEG-1 do tmp1[j] := GT.gmult(p2[j], p1[i]);
// /* and mult (shift) tmp1 right by i */
j:=(MAXDEG*2)-1;
while (j >=i) do begin
tmp1[j] := tmp1[j-i];
dec(j);
end;
for j:=0 to i-1 do tmp1[j] := 0;
///* add into partial product */
for j:=0 to (MAXDEG*2)-1 do dst[j] := dst[j] xor tmp1[j];
end;
end;
end.
