rscode-objc project in GoogleCode
03/26/2010, 16:51 - Open Source
因為想要做點有趣的實驗,會用到 ReedSolomon FEC 演算法,
網路上有 Java 或 C 版本,但沒找到 Objective-C 可直接用的程式碼.
所以就花了一點時間,參考了 Java 和 C 的版本,做了一個 Objective-C 的版本
給 XCode 在 OS X 下使用,並放到 GoogleCode 上作管理.

原來是想直接用 C 的版本,但在改寫的過程中, 發現原來的程式用了很多的 for loop,
及使用了 array 的資料形態.
後來修改為使用 pointer 和 malloc 動態使用記憶體,並用 memset 和 memcpy 減少了 for loop 的使用.

這個 class 應該可以不用修改就可以經由 GNUStep 搬到各個平台下使用.

想要帶回家玩的可以用 svn checkout http://rscode-objc.googlecode.com/svn/trunk/ rscode-objc-read-only 到 GoogleCode 的 svn 下載

原本 C 語言版本的 rscode 是版權是 GPL, 而商業版要收費的, 而我只是修改, 所以也算是 GPL.
不過我想不收錢, 應該也收不到錢, 因為只是修改,所以也不能收錢.
有注意看的話,在程式開頭權利宣告我是用 Copyleft 而不是 Copyright. XD.


參考資料
http://reedsolomon.sourceforge.jp Java 版 ReedSolomon
http://rscode.sourceforge.net C 語言版本.
http://rscode-objc.googlecode.com GoogleCode Project

程式
RSEncodeDecode.h

//
// RSEncodeDecode.h
// RSCode
//
// Created by Tasuka Hsu on 03/25/2010.
// Copyleft 2010 http://Tasuka.IDV.TW.
// Original C Source code is from http://rscode.sourceforge.net
// and reference to Java Source code from http://reedsolomon.sourceforge.jp

#import <Foundation/Foundation.h>
#import <stdio.h>
#import <stdlib.h>

#define PPOLY 0x1d

/*
Below is NPAR, the only compile-time parameter you should have to
modify.
It is the number of parity bytes which will be appended to
your data to create a codeword.
Note that the maximum codeword size is 255, so the
sum of your message length plus parity should be less than
or equal to this maximum limit.
In practice, you will get slooow error correction and decoding
if you use more than a reasonably small number of parity bytes.
(say, 10 or 20)
*/
#define NPAR 17

/* Maximum degree of various polynomials. */
#define MAXDEG (NPAR*2)

#define TABLESIZE 256

@interface RSEncodeDecode : NSObject {
@public
int npar;
int maxdeg;
int length;

int *gExp;
int *gLog;

// Encoder parity bytes
int *pBytes;
// Decoder syndrome bytes
int *synBytes;
// generator polynomial
int *genPoly;
// The Error Locator Polynomial, also known as Lambda or Sigma. Lambda[0]==1
int *Lambda;
// The Error Evaluator Polynomial
int *Omega;

@private
// error locations found using Chien's search
int *ErrorLocs;
int NErrors;

// erasure flags
int *ErasureLocs;
int NErasures;
}

-(void)setNPAR:(int)n;
-(void)setLength:(int)n;
-(BOOL)initTables;
-(void)freeTables;

// Reed Solomon encode/decode routines
-(void)initializeECC;
-(int)checkSyndrome;
-(void)deCode:(unsigned char *)data:(int)l;
-(void)deCode:(unsigned char *)data;
-(void)deCode:(unsigned char *)data:(int)l:(int)n;

-(void)enCode:(unsigned char *)msg:(int)l:(unsigned char *)dst;
-(void)enCode:(unsigned char *)msg:(unsigned char *)dst;
-(void)enCode:(unsigned char *)msg:(int)l:(unsigned char *)dst:(int)nl;

// CRC-CCITT checksum generator
-(short int)crcCCITT:(unsigned char *)msg:(int)len;
-(short int)crcHWare:(short int)data:(short int)genpoly:(short int)accum;

-(void)initGaloisTables;
-(int)gInv:(int)a;
-(int)gMul:(int)a:(int)b;

/* Error location routines */
-(int)correctErrorsErasures:(unsigned char *)codeword:(int)csize:(int)nerasures:(int *)erasures;

-(void)modifiedBerlekampMassey;
/* polynomial arithmetic */
-(void)addPolys:(int *)dst:(int *)src;
-(void)scalePoly:(int)k:(int *)poly;
-(void)mulPolys:(int *)dst:(int *)p1:(int *)p2;
-(void)copyPoly:(int *)dst:(int *)src;
-(void)zeroPoly:(int *)poly;
-(void)findRoots;

/* local ANSI declarations */
-(int)computeDiscrepancy:(int *)lambda:(int *)S:(int)L:(int)n;
-(void)computeGenPoly:(int)l:(int *)genpoly;
-(void)initGamma:(int *)gamma;
-(void)computeModifiedOmega;
-(void)computeNextOmega:(int)d:(int *)A:(int *)dst:(int *)src;
-(void)mulzPoly:(int *)src;
-(void)zeroFillFrom:(unsigned char *)buf:(int)from:(int)to;

// debugging routines
-(void)printParity;
-(void)printSyndrome;
@end


RSEncodeDecode.m

//
// RSEncodeDecode.m
// RSCode
//
// Created by Tasuka Hsu on 03/25/2010.
// Copyrleft 2010 http://Tasuka.IDV.TW.
// Original C Source code from http://rscode.sourceforge.net
// and reference to Java Source code from http://reedsolomon.sourceforge.jp

#import "RSEncodeDecode.h"

@implementation RSEncodeDecode
/* This is one of 14 irreducible polynomials
* of degree 8 and cycle length 255. (Ch 5, pp. 275, Magnetic Recording)
* The high order 1 bit is implicit */
/* x^8 + x^4 + x^3 + x^2 + 1 */

-(id)init
{
if(![super init]){
return nil;
}
npar=NPAR;
maxdeg=npar*2;
length=0;

if([self initTables]!=TRUE){
NSLog(@"Memory allocated error\n");
return nil;
}

[self initializeECC];

return self;
}

-(void)dealloc
{
[self freeTables];
[super dealloc];
}

-(void)setNPAR:(int)n
{
npar=n;
maxdeg=npar*2;

[self freeTables];

if([self initTables]!=TRUE){
NSLog(@"Memory allocated error\n");
return;
}

[self initializeECC];
}

-(void)setLength:(int)n
{
length=n;
}

-(BOOL)initTables
{
if((gExp=(int *)malloc(sizeof(int)*TABLESIZE*2))==NULL){
return FALSE;
}
if((gLog=(int *)malloc(sizeof(int)*TABLESIZE))==NULL){
free(gExp);
return FALSE;
}
if((pBytes=(int *)malloc(sizeof(int)*maxdeg))==NULL){
free(gExp);
free(gLog);
return FALSE;
}
if((synBytes=(int *)malloc(sizeof(int)*maxdeg))==NULL){
free(gExp);
free(gLog);
free(pBytes);
return FALSE;
}
if((Lambda=(int *)malloc(sizeof(int)*maxdeg))==NULL){
free(gExp);
free(gLog);
free(pBytes);
free(synBytes);
return FALSE;
}
if((Omega=(int *)malloc(sizeof(int)*maxdeg))==NULL){
free(gExp);
free(gLog);
free(pBytes);
free(synBytes);
free(Lambda);
return FALSE;
}
if((genPoly=(int *)malloc(sizeof(int)*maxdeg*2))==NULL){
free(gExp);
free(gLog);
free(pBytes);
free(synBytes);
free(Lambda);
free(Omega);
return FALSE;
}
if((ErrorLocs=(int *)malloc(sizeof(int)*maxdeg*2))==NULL){
free(gExp);
free(gLog);
free(pBytes);
free(synBytes);
free(Lambda);
free(Omega);
free(genPoly);
return FALSE;
}
if((ErasureLocs=(int *)malloc(sizeof(int)*maxdeg*2))==NULL){
free(gExp);
free(gLog);
free(pBytes);
free(synBytes);
free(Lambda);
free(Omega);
free(genPoly);
free(ErrorLocs);
return FALSE;
}
return TRUE;
}

-(void)freeTables
{
if(gExp!=NULL){
free(gExp);
}
if(gLog!=NULL){
free(gLog);
}
if(pBytes!=NULL){
free(pBytes);
}
if(synBytes!=NULL){
free(synBytes);
}
if(Lambda!=NULL){
free(Lambda);
}
if(Omega!=NULL){
free(Omega);
}
if(genPoly!=NULL){
free(genPoly);
}
if(ErrorLocs!=NULL){
free(ErrorLocs);
}
if(ErasureLocs!=NULL){
free(ErasureLocs);
}
}

// Initialize lookup tables, polynomials, etc.
-(void)initializeECC
{
// Initialize the galois field arithmetic tables
[self initGaloisTables];
// Compute the encoder generator polynomial
[self computeGenPoly:npar:genPoly];
}

-(void)initGaloisTables
{
// initialize the table of powers of alpha
register int i,d=1;

for(i=0;i<TABLESIZE-1;i++){
gExp=gExp[TABLESIZE-1+i]=d;
gLog[d]=i;
d<<=1;
if((d&0x100)!=0){
d=(d^PPOLY)&0xff;
}
}
}

// multiplication using logarithms
-(int)gMul:(int)a:(int)b
{
return(((a==0)||(b==0))?0:(gExp[gLog[a]+gLog]));
}

-(int)gInv:(int)a
{
return(gExp[TABLESIZE-1-gLog[a]]);
}

// From Cain, Clark, "Error-Correction Coding For Digital Communications", pp. 216.
-(void)modifiedBerlekampMassey
{
register int i,n;
int L,L2,k,d;
int *psi,*psi2,*D,*gamma;

if((psi=(int *)malloc(sizeof(int)*maxdeg))==NULL){
return;
}
if((psi2=(int *)malloc(sizeof(int)*maxdeg))==NULL){
return;
}
if((D=(int *)malloc(sizeof(int)*maxdeg))==NULL){
return;
}
if((gamma=(int *)malloc(sizeof(int)*maxdeg))==NULL){
return;
}

// initialize Gamma, the erasure locator polynomial
[self initGamma:gamma];

// initialize to z
[self copyPoly:D:gamma];
[self mulzPoly:D];

[self copyPoly:psi:gamma];

k=-1;
L=NErasures;

for(n=NErasures;n<npar;n++){
d=[self computeDiscrepancy:psi:synBytes:L:n];
if(d!=0){
// psi2 = psi - d*D
for(i=0;i<maxdeg;i++){
psi2=psi^[self gMul:d:D];
}

if(L<(n-k)){
L2=n-k;
k=n-L;
// D = scalePoly(gInv(d), psi);
for(i=0;i<maxdeg;i++){
D=[self gMul:psi:[self gInv:d]];
}
L=L2;
}

memcpy(psi,psi2,maxdeg*sizeof(int));
}
[self mulzPoly:D];
}

memcpy(Lambda,psi,maxdeg*sizeof(int));

[self computeModifiedOmega];

free(psi);
free(psi2);
free(D);
free(gamma);
}

// given Psi (called Lambda in modifiedBerlekampMassey) and synBytes,
// compute the combined erasure/error evaluator polynomial as Psi*S mod z^4
-(void)computeModifiedOmega
{
int *product;

if((product=(int *)malloc(sizeof(int)*maxdeg*2))==NULL){
return;
}

[self mulPolys:product:Lambda:synBytes];
[self zeroPoly:Omega];

memcpy(Omega,product,npar*sizeof(int));
free(product);
}

/* polynomial multiplication */
-(void)mulPolys:(int *)dst:(int *)p1:(int *)p2
{
register int i,j;
int *tmp;

if((tmp=(int *)malloc(sizeof(int)*maxdeg*2))==NULL){
return;
}

memset(dst,0,maxdeg*2*sizeof(int));

for(i=0;i<maxdeg;i++){
memset((tmp+maxdeg),0,maxdeg*sizeof(int));
// scale tmp1 by p1
for(j=0;j<maxdeg;j++){
tmp[j]=[self gMul:p2[j]:p1];
}
// and mult (shift) tmp1 right by i
for(j=(maxdeg*2)-1;j>=i;j--){
tmp[j]=tmp[j-i];
}
memset(tmp,0,i*sizeof(int));
// add into partial product
for(j=0;j<(maxdeg*2);j++){
dst[j]^=tmp[j];
}
}
free(tmp);
}

// gamma = product (1-z*a^Ij) for erasure locs Ij
-(void)initGamma:(int *)gamma
{
register int i;
int *tmp;

if((tmp=(int *)malloc(sizeof(int)*maxdeg))==NULL){
return;
}

[self zeroPoly:gamma];
[self zeroPoly:tmp];

gamma[0]=1;

for(i=0;i<NErasures;i++){
[self copyPoly:tmp:gamma];
[self scalePoly:gExp[ErasureLocs]:tmp];
[self mulzPoly:tmp];
[self addPolys:gamma:tmp];
}
free(tmp);
}

-(void)computeNextOmega:(int)d:(int *)A:(int *)dst:(int *)src
{
register int i;

for(i=0;i<maxdeg;i++){
dst=src^[self gMul:d:A];
}
}

-(int)computeDiscrepancy:(int *)lambda:(int *)S:(int)L:(int)n
{
register int i,sum=0;

for(i=0;i<=L;i++){
sum^=[self gMul:lambda:S[n-i]];
}

return(sum);
}

// polynomial arithmetic
-(void)addPolys:(int *)dst:(int *)src
{
register int i;

for(i=0;i<maxdeg;i++){
dst^=src;
}
}

-(void)copyPoly:(int *)dst:(int *)src
{
memcpy(dst,src,maxdeg*sizeof(int));
}

-(void)scalePoly:(int)k:(int *)poly
{
register int i;

for(i=0;i<maxdeg;i++){
poly=[self gMul:k:poly];
}
}

-(void)zeroPoly:(int *)poly
{
memset(poly,0,maxdeg*sizeof(int));
}

// multiply by z, i.e., shift right by 1
-(void)mulzPoly:(int *)src
{
register int i;

for(i=maxdeg-1;i>0;i--){
src=src[i-1];
}
src[0]=0;
}

/* Finds all the roots of an error-locator polynomial with coefficients
* Lambda[j] by evaluating Lambda at successive values of alpha.
*
* This can be tested with the decoder's equations case.
*/

-(void)findRoots
{
int sum;
register int i,j;

NErrors=0;

for(i=1;i<TABLESIZE;i++){
sum=0;
// evaluate lambda at i
for(j=0;j<npar+1;j++){
sum^=[self gMul:gExp[(i*j)%(TABLESIZE-1)]:Lambda[j]];
}
if(sum==0){
ErrorLocs[NErrors]=(TABLESIZE-1-i);
NErrors++;
#ifdef DEBUG
fprintf(stderr,"Root found at i = %d, (255-i) = %d\n",i,(TABLESIZE-1-i));
#endif
}
}
}

/* Combined Erasure And Error Magnitude Computation
*
* Pass in the codeword, its size in bytes, as well as
* an array of any known erasure locations, along the number
* of these erasures.
*
* Evaluate Omega(actually Psi)/Lambda' at the roots
* alpha^(-i) for error locs i.
*
* Returns 1 if everything ok, or 0 if an out-of-bounds error is found
*
*/

-(int)correctErrorsErasures:(unsigned char *)codeword:(int)csize:(int)nerasures:(int *)erasures
{
int r,i,j,err;
int num,denom;

// If you want to take advantage of erasure correction, be sure to
// set NErasures and ErasureLocs[] with the locations of erasures.

NErasures=nerasures;

memcpy(ErasureLocs,erasures,NErasures);

[self modifiedBerlekampMassey];
[self findRoots];

if((NErrors<=npar)&&(NErrors>0)){
// first check for illegal error locs
for(r=0;r<NErrors;r++) {
if(ErrorLocs[r]>=csize){
#ifdef DEBUG
fprintf(stderr,"Error loc i=%d outside of codeword length %d\n",i,csize);
#endif
return(0);
}
}

for(r=0;r<NErrors;r++){
i=ErrorLocs[r];
// evaluate Omega at alpha^(-i)
num=0;
for(j=0;j<maxdeg;j++){
num^=[self gMul:Omega[j]:gExp[((TABLESIZE-1-i)*j)%(TABLESIZE-1)]];
}
// evaluate Lambda' (derivative) at alpha^(-i) ; all odd powers disappear
denom=0;
for(j=1;j<maxdeg;j+=2){
denom^=[self gMul:Lambda[j]:gExp[((TABLESIZE-1-i)*(j-1))%(TABLESIZE-1)]];
}

err=[self gMul:num:[self gInv:denom]];
#ifdef DEBUG
fprintf(stderr,"Error magnitude %#x at loc %d\n",err,csize-i);
#endif
codeword[csize-i-1]^=err;
}
return(1);
}else{
#ifdef DEBUG
if(NErrors){
fprintf(stderr,"Uncorrectable codeword\n");
}
#endif
return(0);
}
}

// models crc hardware (minor variation on polynomial division algorithm)
-(short int)crcHWare:(short int)data:(short int)genpoly:(short int)accum
{
register short int i;

data<<=8;

for(i=8;i>0;i--){
if((data^accum)&0x8000){
accum=((accum<<1)^genpoly)&0xFFFF;
}else{
accum=(accum<<1)&0xFFFF;
}
data=(data<<1)&0xFFFF;
}
return(accum);
}

// Computes the CRC-CCITT checksum on array of byte data, length len
-(short int)crcCCITT:(unsigned char *)msg:(int)len
{
register int i;
register short int acc=0;

for(i=0;i<len;i++){
acc=[self crcHWare:(short int)msg:(short int)0x1021:acc];
}
return(acc);
}

-(void)zeroFillFrom:(unsigned char *)buf:(int)from:(int)to
{
register int i;

for(i=from;i<to;i++){
buf=0;
}
}

// debugging routines
-(void)printParity
{
register int i;

printf("Parity Bytes: ");
for(i=0;i<npar;i++){
printf("[%d]:%x, ",i,pBytes);
}
printf("\n");
}

-(void)printSyndrome
{
register int i;

printf("Syndrome Bytes: ");
for(i=0;i<npar;i++){
printf("[%d]:%x, ",i,synBytes);
}
printf("\n");
}

// Append the parity bytes onto the end of the message
-(void)build_codeword:(unsigned char *)msg:(int)l:(unsigned char *)dst
{
register int i;

memcpy(dst,msg,l*sizeof(int));
for(i=0;i<npar;i++){
dst[i+l]=pBytes[npar-1-i];
}
}

/*
* Reed Solomon Decoder
*
* Computes the syndrome of a codeword. Puts the results
* into the synBytes[] array.
*/

-(void)deCode:(unsigned char *)data:(int)l
{
register int i,j;
int sum;

[self setLength:l];

for(j=0;j<npar;j++){
sum=0;
for(i=0;i<l;i++){
sum=data^[self gMul:gExp[j+1]:sum];
}
synBytes[j]=sum;
}
}

-(void)deCode:(unsigned char *)data
{
[self deCode:data:length];
}

-(void)deCode:(unsigned char *)data:(int)l:(int)n
{
[self setNPAR:n];
[self setLength:l];
[self deCode:data:length];
}

// Check if the syndrome is zero
-(int)checkSyndrome
{
register int i;
int nz=0;

for(i=0;i<npar;i++){
if(synBytes!=0){
nz=1;
break;
}
}
return nz;
}

-(void)debugCheckSyndrome
{
register int i;

for(i=0;i<3;i++){
printf(" inv log S[%d]/S[%d] = %d\n",i,i+1, \
gLog[[self gMul:synBytes:[self gInv:synBytes[i+1]]]]);
}
}

/* Create a generator polynomial for an n byte RS code.
* The coefficients are returned in the genPoly arg.
* Make sure that the genPoly array which is passed in is
* at least n+1 bytes long.
*/

-(void)computeGenPoly:(int)l:(int *)genpoly
{
register int i;
int *tp,*tp1;

if((tp=(int *)malloc(sizeof(int)*TABLESIZE))==NULL){
return;
}
if((tp1=(int *)malloc(sizeof(int)*TABLESIZE))==NULL){
return;
}

// multiply (x + a^n) for n = 1 to l
[self zeroPoly:tp1];
tp1[0]=1;

for(i=1;i<=l;i++){
[self zeroPoly:tp];
tp[0]=gExp; /* set up x+a^n */
tp[1]=1;

[self mulPolys:genpoly:tp:tp1];
[self copyPoly:tp1:genpoly];
}
free(tp);
free(tp1);
}

/* Simulate a LFSR with generator polynomial for n byte RS code.
* Pass in a pointer to the data array, and amount of data.
*
* The parity bytes are deposited into pBytes[], and the whole message
* and parity are copied to dest to make a codeword.
*
*/

-(void)enCode:(unsigned char *)msg:(int)l:(unsigned char *)dst
{
register int i,j;
int *LFSR,dbyte;

if((LFSR=(int *)malloc(sizeof(int)*(npar+1)))==NULL){
return;
}

memset(LFSR,0,npar*sizeof(int));

[self setLength:l];

for(i=0;i<l;i++){
dbyte=msg^LFSR[npar-1];

for(j=npar-1;j>0;j--){
LFSR[j]=LFSR[j-1]^[self gMul:genPoly[j]:dbyte];
}
LFSR[0]=[self gMul:genPoly[0]:dbyte];
}
memcpy(pBytes,LFSR,npar*sizeof(int));

[self build_codeword:msg:l:dst];

free(LFSR);
}

-(void)enCode:(unsigned char *)msg:(unsigned char *)dst
{
[self enCode:msg:length:dst];
}

-(void)enCode:(unsigned char *)msg:(int)l:(unsigned char *)dst:(int)n
{
[self setNPAR:n];
[self setLength:l];
[self enCode:msg:length:dst];
}

@end


測試程式
RSCode.m

// Created by Tasuka Hsu on 03/25/2010.
// Copyleft 2010 http://Tasuka.IDV.TW.
// Original C Source code is from http://rscode.sourceforge.net

#import <Foundation/Foundation.h>
#import "RSEncodeDecode.h"

#define ML (sizeof(msg)+NPAR)

// Introduce a byte error at LOC
void byte_err(int err,int loc,unsigned char *dst)
{
//printf("Adding Error at loc %d, data %#x\n",loc,dst[loc-1]);
dst[loc-1]^=err;
}

// Pass in location of error (first byte position is labeled starting at 1, not 0), and the codeword.
void byte_erasure(int loc,unsigned char *dst,int cwsize,int *erasures)
{
//printf("Erasure at loc %d, data %#x\n",loc,dst[loc-1]);
dst[loc-1]=0;
}

int main (int argc, const char * argv[]) {
NSAutoreleasePool *pool=[[NSAutoreleasePool alloc] init];

unsigned char msg[]="Learn from other people's mistakes, you do'nt have time to make your own.";
unsigned char codeword[TABLESIZE];

int erasures[16];
int nerasures=0;

// Initialize RSCode Encode Decode Object
RSEncodeDecode *RS=[[RSEncodeDecode alloc] init];

[RS setNPAR:17];
[RS setLength:sizeof(msg)];

// Encode data into codeword, adding NPAR parity bytes
//[RS enCode:msg:sizeof(msg):codeword];
[RS enCode:msg:codeword];

printf("Encoded data is: \"%s\"\n",codeword);

// Add errors and two erasures
byte_err(0x35,3,codeword);
byte_err(0x41,6,codeword);
byte_err(0x40,10,codeword);

byte_err(0x23,17,codeword);
byte_err(0x34,19,codeword);

byte_err(0x30,53,codeword);
byte_err(0x32,65,codeword);

printf("with some errors: \"%s\"\n",codeword);

// We need to indicate the position of the erasures.
// Eraseure positions are indexed (1 based) from the end of the message...
erasures[nerasures++]=ML-16;
erasures[nerasures++]=ML-54;

// Now decode -- encoded codeword size must be passed
[RS deCode:codeword:ML];

// check if syndrome is all zeros
if([RS checkSyndrome]!=0){
[RS correctErrorsErasures:codeword:ML:nerasures:erasures];
printf("Corrected codeword: \"%s\"\n",codeword);
}
// Release RSCode Encode Decode Object
[RS release];

[pool drain];
return 0;
}


回應

發表回應

填寫下面來發表回應。









插入項目:


遊覽已上傳的圖片








回應需經過管理人員認可後才會出現在網頁上.