diff --git a/Release/qra64.exe b/Release/qra64.exe new file mode 100644 index 0000000..c17c76f Binary files /dev/null and b/Release/qra64.exe differ diff --git a/Release/qra65.exe b/Release/qra65.exe index 1c86dd7..49fb581 100644 Binary files a/Release/qra65.exe and b/Release/qra65.exe differ diff --git a/Release/qracodes-mt.exe b/Release/qracodes-mt.exe deleted file mode 100644 index e178a99..0000000 Binary files a/Release/qracodes-mt.exe and /dev/null differ diff --git a/Release/qracodes-st.exe b/Release/qracodes-st.exe new file mode 100644 index 0000000..d9f86bc Binary files /dev/null and b/Release/qracodes-st.exe differ diff --git a/Release/qracodes.exe b/Release/qracodes.exe index cbb54a4..31afb10 100644 Binary files a/Release/qracodes.exe and b/Release/qracodes.exe differ diff --git a/qracodes-mt.sln b/qra64.sln similarity index 54% rename from qracodes-mt.sln rename to qra64.sln index 704f899..d6d1655 100644 --- a/qracodes-mt.sln +++ b/qra64.sln @@ -1,7 +1,7 @@  Microsoft Visual Studio Solution File, Format Version 10.00 # Visual Studio 2008 -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "qracodes-mt", "qracodes-mt\qracodes-mt.vcproj", "{272D3997-2EE4-4D21-9F03-CE07F40B6134}" +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "qra64", "qra64\qra64.vcproj", "{8C6ACED4-B254-4959-9559-6A9F2306391F}" EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution @@ -9,10 +9,10 @@ Global Release|Win32 = Release|Win32 EndGlobalSection GlobalSection(ProjectConfigurationPlatforms) = postSolution - {272D3997-2EE4-4D21-9F03-CE07F40B6134}.Debug|Win32.ActiveCfg = Debug|Win32 - {272D3997-2EE4-4D21-9F03-CE07F40B6134}.Debug|Win32.Build.0 = Debug|Win32 - {272D3997-2EE4-4D21-9F03-CE07F40B6134}.Release|Win32.ActiveCfg = Release|Win32 - {272D3997-2EE4-4D21-9F03-CE07F40B6134}.Release|Win32.Build.0 = Release|Win32 + {8C6ACED4-B254-4959-9559-6A9F2306391F}.Debug|Win32.ActiveCfg = Debug|Win32 + {8C6ACED4-B254-4959-9559-6A9F2306391F}.Debug|Win32.Build.0 = Debug|Win32 + {8C6ACED4-B254-4959-9559-6A9F2306391F}.Release|Win32.ActiveCfg = Release|Win32 + {8C6ACED4-B254-4959-9559-6A9F2306391F}.Release|Win32.Build.0 = Release|Win32 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/qra64.suo b/qra64.suo new file mode 100644 index 0000000..872effb Binary files /dev/null and b/qra64.suo differ diff --git a/qra64/Makefile.Win b/qra64/Makefile.Win new file mode 100644 index 0000000..7bd10c2 --- /dev/null +++ b/qra64/Makefile.Win @@ -0,0 +1,30 @@ +FC = gfortran +CC = gcc +CFLAGS = -O2 -Wall -I. -D_WIN32 + +# Default rules +%.o: %.c + ${CC} ${CFLAGS} -c $< +%.o: %.f + ${FC} ${FFLAGS} -c $< +%.o: %.F + ${FC} ${FFLAGS} -c $< +%.o: %.f90 + ${FC} ${FFLAGS} -c $< +%.o: %.F90 + ${FC} ${FFLAGS} -c $< + +all: qra64.exe + +OBJS1 = main.o qra64.o +qra64.exe: $(OBJS1) + ${CC} -o qra64.exe $(OBJS1) ../qracodes/libqra64.a -lm + +OBJS2 = qra64sim.o options.o wavhdr.o +qra64sim.exe: $(OBJS2) + ${FC} -o qra64sim.exe $(OBJS2) ../qracodes/libqra64.a -lm + +.PHONY : clean + +clean: + $(RM) *.o qra64.exe qra64sim.exe diff --git a/qra64/build.sh b/qra64/build.sh new file mode 100644 index 0000000..62bac0f --- /dev/null +++ b/qra64/build.sh @@ -0,0 +1,10 @@ +g++ -march=native -pthread -O3 -I../qracodes/ \ +*.c \ +../qracodes/normrnd.c \ +../qracodes/npfwht.c \ +../qracodes/pdmath.c \ +../qracodes/qra12_63_64_irr_b.c \ +../qracodes/qra13_64_64_irr_e.c \ +../qracodes/qracodes.c \ +-lpthread -o qra64 + diff --git a/qra64/main.c b/qra64/main.c new file mode 100644 index 0000000..85b57ff --- /dev/null +++ b/qra64/main.c @@ -0,0 +1,458 @@ +/* +main.c +QRA64 mode encode/decode test + +(c) 2016 - Nico Palermo, IV3NWV + +Thanks to Andrea Montefusco IW0HDV for his help on adapting the sources +to OSs other than MS Windows + +------------------------------------------------------------------------------ +This file is part of the qracodes project, a Forward Error Control +encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. + +Files in this package: + main.c - this file + qra64.c/.h - qra64 mode encode/decoding functions + + ../qracodes/normrnd.{c,h} - random gaussian number generator + ../qracodes/npfwht.{c,h} - Fast Walsh-Hadamard Transforms + ../qracodes/pdmath.{c,h} - Elementary math on probability distributions + ../qracodes/qra12_63_64_irr_b.{c,h} - Tables for a QRA(12,63) irregular RA + code over GF(64) + ../qracodes/qra13_64_64_irr_e.{c,h} - Tables for a QRA(13,64) irregular RA + code over GF(64) + ../qracodes/qracodes.{c,h} - QRA codes encoding/decoding functions + +------------------------------------------------------------------------------- + + qracodes 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. + qracodes 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 qracodes source distribution. + If not, see . + +----------------------------------------------------------------------------- + +The code used by the QRA64 mode is the code: QRA13_64_64_IRR_E: K=13 +N=64 Q=64 irregular QRA code (defined in qra13_64_64_irr_e.{h,c}). + +This code has been designed to include a CRC as the 13th information +symbol and improve the code UER (Undetected Error Rate). The CRC +symbol is not sent along the channel (the codes are punctured) and the +resulting code is still a (12,63) code with an effective code rate of +R = 12/63. +*/ + +// OS dependent defines and includes ------------------------------------------ + +#if _WIN32 // note the underscore: without it, it's not msdn official! +// Windows (x64 and x86) +#include // required only for GetTickCount(...) +#include // _beginthread +#endif + +#if __linux__ +#include +#include + +unsigned GetTickCount(void) { + struct timespec ts; + unsigned theTick = 0U; + clock_gettime( CLOCK_REALTIME, &ts ); + theTick = ts.tv_nsec / 1000000; + theTick += ts.tv_sec * 1000; + return theTick; +} +#endif + +#if __APPLE__ +#endif + +#include +#include +#include + +#include "qra64.h" +#include "../qracodes/normrnd.h" // gaussian numbers generator + +// ---------------------------------------------------------------------------- + +// channel types +#define CHANNEL_AWGN 0 +#define CHANNEL_RAYLEIGH 1 + +void printwordd(char *msg, int *x, int size) +{ + int k; + printf("\n%s ",msg); + for (k=0;k=0) { // decoded + printf("K1JT rx: received with apcode=%d %s\n",rc, decode_type[rc]); + +// Step 2a: K1JT replies to IV3NWV (with no grid) + printf("K1JT tx: IV3NWV K1JT\n"); + encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_BLANK); + qra64_encode(codec_k1jt, y, x); + rx = mfskchannel(y,channel_type,EbNodB); + +// Step 2b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?] or [IV3NWV ?] + rc = qra64_decode(codec_iv3nwv, xdec,rx); + if (rc>=0) { // decoded + printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]); + +// Step 3a: IV3NWV replies to K1JT with a 73 + printf("IV3NWV tx: K1JT IV3NWV 73\n"); + encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV, GRID_73); + qra64_encode(codec_iv3nwv, y, x); + rx = mfskchannel(y,channel_type,EbNodB); + +// Step 3b: K1JT attempts to decode [? ? ?] or [K1JT IV3NWV ?] + rc = qra64_decode(codec_k1jt, xdec,rx); + if (rc>=0) { // decoded + printf("K1JT rx: received with apcode=%d %s\n",rc, decode_type[rc]); + +// Step 4a: K1JT replies to IV3NWV with a 73 + printf("K1JT tx: IV3NWV K1JT 73\n"); + encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_73); + qra64_encode(codec_k1jt, y, x); + rx = mfskchannel(y,channel_type,EbNodB); + +// Step 4b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?], or [IV3NWV ?] + rc = qra64_decode(codec_iv3nwv, xdec,rx); + if (rc>=0) { // decoded + printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]); + return 0; + } + } + } + } + printf("the other party did not decode\n"); + return -1; +} + +int test_proc_2(int channel_type, float EbNodB, int mode) +{ +/* +Here we simulate the decoder of K1JT after K1JT has sent a msg [IV3NWV K1JT] +and IV3NWV sends him the msg [K1JT IV3NWV JN66]. + +If mode=QRA_NOAP, K1JT decoder attempts to decode only msgs of type [? ? ?]. + +If mode=QRA_AUTOP, K1JT decoder will attempt to decode also the msgs +[K1JT IV3NWV] and [K1JT IV3NWV ?]. + +In the case a decode is successful the return code of the qra64_decode function +indicates the amount of a-priori information required to decode the received +message according to this table: + + rc=0 [? ? ?] AP0 + rc=1 [CQ ? ?] AP27 + rc=2 [CQ ? ] AP42 + rc=3 [CALL ? ?] AP29 + rc=4 [CALL ? ] AP44 + rc=5 [CALL CALL ?] AP57 + +The return code is <0 when decoding is unsuccessful + +This test simulates the situation ntx times and reports how many times +a particular type decode among the above 6 cases succeded. +*/ + + int x[QRA64_K], xdec[QRA64_K]; + int y[QRA64_N]; + float *rx; + int rc,k; + + int ndecok[6] = { 0, 0, 0, 0, 0, 0}; + int ntx = 100,ndec=0; + + qra64codec *codec_iv3nwv = qra64_init(mode,CALL_IV3NWV); // codec for IV3NWV + qra64codec *codec_k1jt = qra64_init(mode,CALL_K1JT); // codec for K1JT + +// This will enable K1JT's decoder to look for IV3NWV calls + encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT,GRID_BLANK); + qra64_encode(codec_k1jt, y, x); + printf("K1JT tx: IV3NWV K1JT\n"); + + // IV3NWV reply to K1JT + printf("IV3NWV tx: K1JT IV3NWV JN66\n"); + encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV,GRID_JN66); + qra64_encode(codec_iv3nwv, y, x); + + printf("Simulating decodes by K1JT up to AP56 ..."); + + for (k=0;k=0) + ndecok[rc]++; + } + printf("\n"); + + printf("Transimtted:%d - Decoded:\n",ntx); + for (k=0;k<6;k++) { + printf("%3d with %s\n",ndecok[k],decode_type[k]); + ndec += ndecok[k]; + } + printf("Total: %d/%d\n",ndec,ntx); + printf("\n"); + + return 0; +} + +void syntax(void) +{ + printf("\nQRA64 Mode Tests\n"); + printf("2016, Nico Palermo - IV3NWV\n\n"); + printf("---------------------------\n\n"); + printf("Syntax: qra64 [-s] [-c] [-a] [-t] [-h]\n"); + printf("Options: \n"); + printf(" -s : set simulation SNR in 2500 Hz BW (default:-27.5 dB)\n"); + printf(" -c : set channel type 0=AWGN (default) 1=Rayleigh\n"); + printf(" -a : set decode type 0=NO_AP 1=AUTO_AP (default)\n"); + printf(" -t: 0=simulate seq of msgs between IV3NWV and K1JT (default)\n"); + printf(" 1=simulate K1JT receiving K1JT IV3NWV JN66\n"); + printf(" -h: this help\n"); +} + +int main(int argc, char* argv[]) +{ + int k, rc, nok=0; + float SNRdB = -27.5f; + unsigned int channel = CHANNEL_AWGN; + unsigned int mode = QRA_AUTOAP; + unsigned int testtype=0; + int nqso = 100; + float EbNodB; + +// Parse the command line + while(--argc) { + argv++; + if (strncmp(*argv,"-h",2)==0) { + syntax(); + return 0; + } else { + if (strncmp(*argv,"-a",2)==0) { + mode = ( int)atoi((*argv)+2); + if (mode>1) { + printf("Invalid decoding mode\n"); + syntax(); + return -1; + } + } else { + if (strncmp(*argv,"-s",2)==0) { + SNRdB = (float)atof((*argv)+2); + if (SNRdB>0 || SNRdB<-40) { + printf("SNR should be in the range [-40..0]\n"); + syntax(); + return -1; + } + } else { + if (strncmp(*argv,"-t",2)==0) { + testtype = ( int)atoi((*argv)+2); + if (testtype>1) { + printf("Invalid test type\n"); + syntax(); + return -1; + } + } else { + if (strncmp(*argv,"-c",2)==0) { + channel = ( int)atoi((*argv)+2); + if (channel>CHANNEL_RAYLEIGH) { + printf("Invalid channel type\n"); + syntax(); + return -1; + } + } else { + printf("Invalid option\n"); + syntax(); + return -1; + } + } + } + } + } + } + + EbNodB = SNRdB+29.1f; + +#if defined(__linux__) || defined(__unix__) + srand48(GetTickCount()); +#endif + + if (testtype==0) { + for (k=0;k. + +----------------------------------------------------------------------------- + +Code used in this sowftware release: + +QRA13_64_64_IRR_E: K=13 N=64 Q=64 irregular QRA code (defined in +qra13_64_64_irr_e.h /.c) + +Codes with K=13 are designed to include a CRC as the 13th information symbol +and improve the code UER (Undetected Error Rate). +The CRC symbol is not sent along the channel (the codes are punctured) and the +resulting code is a (12,63) code +*/ +//---------------------------------------------------------------------------- + +#include +#include + +#include "qra64.h" +#include "../qracodes/qracodes.h" +#include "../qracodes/qra13_64_64_irr_e.h" +#include "../qracodes/pdmath.h" + +// Code parameters of the QRA64 mode +#define QRA64_CODE qra_13_64_64_irr_e +#define QRA64_NMSG 218 // Must much value indicated in QRA64_CODE.NMSG + +#define QRA64_KC (QRA64_K+1) // Information symbols (crc included) +#define QRA64_NC (QRA64_N+1) // Codeword length (as defined in the code) +#define QRA64_NITER 100 // max number of iterations per decode + +// static functions declarations ---------------------------------------------- +static int calc_crc6(const int *x, int sz); +static void ix_mask(float *dst, const float *src, const int *mask, + const int *x); +static int qra64_do_decode(int *x, const float *pix, const int *ap_mask, + const int *ap_x); + +// a-priori information masks for fields in JT65-like msgs -------------------- +#define MASK_CQQRZ 0xFFFFFFC // CQ/QRZ calls common bits +#define MASK_CALL1 0xFFFFFFF +#define MASK_CALL2 0xFFFFFFF +#define MASK_GRIDFULL 0xFFFF +#define MASK_GRIDBIT 0x8000 // b[15] is 1 for free text, 0 otherwise +// ---------------------------------------------------------------------------- + +qra64codec *qra64_init(int flags, const int mycall) +{ + + // Eb/No value for which we optimize the decoder metric + const float EbNodBMetric = 2.8f; + const float EbNoMetric = (float)pow(10,EbNodBMetric/10); + const float R = 1.0f*(QRA64_KC)/(QRA64_NC); + + qra64codec *pcodec = (qra64codec*)malloc(sizeof(qra64codec)); + + if (!pcodec) + return 0; // can't allocate memory + + pcodec->decEsNoMetric = 1.0f*QRA64_m*R*EbNoMetric; + pcodec->apflags = flags; + + if (flags!=QRA_AUTOAP) + return pcodec; + + // initialize messages and mask for decoding with a-priori information + + pcodec->apmycall = mycall; + pcodec->apsrccall = 0; + + // encode CQ/QRZ messages and masks + // NOTE: Here we handle only CQ and QRZ msgs + // 'CQ nnn', 'CQ DX' and 'DE' msgs + // will be handled by the decoder as messages with no a-priori knowledge + encodemsg_jt65(pcodec->apmsg_cqqrz, CALL_CQ, 0, GRID_BLANK); + encodemsg_jt65(pcodec->apmask_cqqrz, MASK_CQQRZ,0, MASK_GRIDBIT); // AP27 + encodemsg_jt65(pcodec->apmask_cqqrz_ooo, MASK_CQQRZ,0, MASK_GRIDFULL);// AP42 + + // encode [mycall ? x] messages and set masks + encodemsg_jt65(pcodec->apmsg_call1, mycall, 0, GRID_BLANK); + encodemsg_jt65(pcodec->apmask_call1, MASK_CALL1, 0, MASK_GRIDBIT); // AP29 + encodemsg_jt65(pcodec->apmask_call1_ooo, MASK_CALL1,0, MASK_GRIDFULL);// AP44 + + // set mask for [mycall srccall ?] messages + encodemsg_jt65(pcodec->apmask_call1_call2,MASK_CALL1,MASK_CALL2, + MASK_GRIDBIT); // AP56 + return pcodec; +} + +void qra64_encode(qra64codec *pcodec, int *y, const int *x) +{ + int encx[QRA64_KC]; // encoder input buffer + int ency[QRA64_NC]; // encoder output buffer + + int call1,call2,grid; + + memcpy(encx,x,QRA64_K*sizeof(int)); // Copy input to encoder buffer + encx[QRA64_K]=calc_crc6(encx,QRA64_K); // Compute and add crc symbol + qra_encode(&QRA64_CODE, ency, encx); // encode msg+crc using given QRA code + + // copy codeword to output puncturing the crc symbol + memcpy(y,ency,QRA64_K*sizeof(int)); // copy information symbols + memcpy(y+QRA64_K,ency+QRA64_KC,QRA64_C*sizeof(int)); // copy parity symbols + + if (pcodec->apflags!=QRA_AUTOAP) + return; + + // look if the msg sent is a std type message (bit15 of grid field = 0) + if ((x[9]&0x80)==1) + return; // no, it's a text message + + // It's a [call1 call2 grid] message + + // We assume that call2 is our call (but we don't check it) + // call1 the station callsign we are calling or indicates a general call (CQ/QRZ/etc..) + decodemsg_jt65(&call1,&call2,&grid,x); + + if ((call1>=CALL_CQ && call1<=CALL_CQ999) || call1==CALL_CQDX || + call1==CALL_DE) { + // We are making a general call; don't know who might reply (srccall) + // Reset apsrccall to 0 so decoder won't look for [mycall srccall ?] msgs + pcodec->apsrccall = 0; + } else { + // We are replying to someone named call1 + // Set apmsg_call1_call2 so decoder will try for [mycall call1 ?] msgs + pcodec->apsrccall = call1; + encodemsg_jt65(pcodec->apmsg_call1_call2, pcodec->apmycall, + pcodec->apsrccall, 0); + } +} + +int qra64_decode(qra64codec *pcodec, int *x, const float *rxen) +{ + int k; + float *srctmp, *dsttmp; + float ix[QRA64_NC*QRA64_M]; // (depunctured) intrisic information + int rc; + + if (QRA64_NMSG!=QRA64_CODE.NMSG) // sanity check + return -16; // QRA64_NMSG define is wrong + + // compute symbols intrinsic probabilities from received energy observations + qra_mfskbesselmetric(ix, rxen, QRA64_m, QRA64_N,pcodec->decEsNoMetric); + + // de-puncture observations adding a uniform distribution for the crc symbol + + // move check symbols distributions one symbol towards the end + dsttmp = PD_ROWADDR(ix,QRA64_M, QRA64_NC-1); //Point to last symbol prob dist + srctmp = dsttmp-QRA64_M; // source is the previous pd + for (k=0;k=0) return 0; // successfull decode with AP0 + + if (pcodec->apflags!=QRA_AUTOAP) return rc; // rc<0 = unsuccessful decode + + // Attempt to decode CQ calls + rc = qra64_do_decode(x,ix,pcodec->apmask_cqqrz, pcodec->apmsg_cqqrz); // AP27 + if (rc>=0) return 1; // decoded [cq/qrz ? ?] + + rc = qra64_do_decode(x, ix, pcodec->apmask_cqqrz_ooo, + pcodec->apmsg_cqqrz); // AP42 + if (rc>=0) return 2; // decoded [cq ? ooo] + + // attempt to decode calls directed to us (mycall) + rc = qra64_do_decode(x, ix, pcodec->apmask_call1, + pcodec->apmsg_call1); // AP29 + if (rc>=0) return 3; // decoded [mycall ? ?] + + rc = qra64_do_decode(x, ix, pcodec->apmask_call1_ooo, + pcodec->apmsg_call1); // AP44 + if (rc>=0) return 4; // decoded [mycall ? ooo] + + // if apsrccall is set attempt to decode [mycall srccall ?] msgs + if (pcodec->apsrccall==0) return rc; // nothing more to do + + rc = qra64_do_decode(x, ix, pcodec->apmask_call1_call2, + pcodec->apmsg_call1_call2); // AP57 + if (rc>=0) return 5; // decoded [mycall srccall ?] + + return rc; +} + +// Static functions definitions ---------------------------------------------- + +// Decode with given a-priori information +static int qra64_do_decode(int *x, const float *pix, const int *ap_mask, + const int *ap_x) +{ + int rc; + const float *ixsrc; + float ix_masked[QRA64_NC*QRA64_M]; // Masked intrinsic information + float ex[QRA64_NC*QRA64_M]; // Extrinsic information from the decoder + + float v2cmsg[QRA64_NMSG*QRA64_M]; // buffers for the decoder messages + float c2vmsg[QRA64_NMSG*QRA64_M]; + int xdec[QRA64_KC]; + + if (ap_mask==NULL) { // no a-priori information + ixsrc = pix; // intrinsic source is what passed as argument + } else { + // a-priori information provided + // mask channel observations with a-priori + ix_mask(ix_masked,pix,ap_mask,ap_x); + ixsrc = ix_masked; // intrinsic source is the masked version + } + + // run the decoding algorithm + rc = qra_extrinsic(&QRA64_CODE,ex,ixsrc,QRA64_NITER,v2cmsg,c2vmsg); + if (rc<0) + return -1; // no convergence in given iterations + + // decode + qra_mapdecode(&QRA64_CODE,xdec,ex,ixsrc); + + // verify crc + if (calc_crc6(xdec,QRA64_K)!=xdec[QRA64_K]) // crc doesn't match (detected error) + return -2; // decoding was succesfull but crc doesn't match + + // success. copy decoded message to output buffer + memcpy(x,xdec,QRA64_K*sizeof(int)); + + return 0; +} +// crc functions -------------------------------------------------------------- +// crc-6 generator polynomial +// g(x) = x^6 + a5*x^5 + ... + a1*x + a0 + +// g(x) = x^6 + x + 1 +#define CRC6_GEN_POL 0x30 // MSB=a0 LSB=a5 + +// g(x) = x^6 + x^2 + x + 1 (See: https://users.ece.cmu.edu/~koopman/crc/) +// #define CRC6_GEN_POL 0x38 // MSB=a0 LSB=a5. Simulation results are similar + +static int calc_crc6(const int *x, int sz) +{ + // todo: compute it faster using a look up table + int k,j,t,sr = 0; + for (k=0;k>1) ^ CRC6_GEN_POL; + else + sr = (sr>>1); + t>>=1; + } + } + return sr; +} + +static void ix_mask(float *dst, const float *src, const int *mask, + const int *x) +{ + // mask intrinsic information (channel observations) with a priori knowledge + + int k,kk, smask; + float *row; + + memcpy(dst,src,(QRA64_NC*QRA64_M)*sizeof(float)); + + for (k=0;k>22)&0x3F; + y[1]= (call1>>16)&0x3F; + y[2]= (call1>>10)&0x3F; + y[3]= (call1>>4)&0x3F; + y[4]= (call1<<2)&0x3F; + + y[4] |= (call2>>26)&0x3F; + y[5]= (call2>>20)&0x3F; + y[6]= (call2>>14)&0x3F; + y[7]= (call2>>8)&0x3F; + y[8]= (call2>>2)&0x3F; + y[9]= (call2<<4)&0x3F; + + y[9] |= (grid>>12)&0x3F; + y[10]= (grid>>6)&0x3F; + y[11]= (grid)&0x3F; + +} +void decodemsg_jt65(int *call1, int *call2, int *grid, const int *x) +{ + int nc1, nc2, ng; + + nc1 = x[4]>>2; + nc1 |= x[3]<<4; + nc1 |= x[2]<<10; + nc1 |= x[1]<<16; + nc1 |= x[0]<<22; + + nc2 = x[9]>>4; + nc2 |= x[8]<<2; + nc2 |= x[7]<<8; + nc2 |= x[6]<<14; + nc2 |= x[5]<<20; + nc2 |= (x[4]&0x03)<<26; + + ng = x[11]; + ng |= x[10]<<6; + ng |= (x[9]&0x0F)<<12; + + *call1 = nc1; + *call2 = nc2; + *grid = ng; +} diff --git a/qra64/qra64.h b/qra64/qra64.h new file mode 100644 index 0000000..8fd7cc5 --- /dev/null +++ b/qra64/qra64.h @@ -0,0 +1,123 @@ +// qra64.h +// Encoding/decoding functions for the QRA64 mode +// +// (c) 2016 - Nico Palermo, IV3NWV +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qracodes source distribution. +// If not, see . + +#ifndef _qra64_h_ +#define _qra64_h_ + +// qra64_init(...) initialization flags +#define QRA_NOAP 0 // don't use a-priori knowledge +#define QRA_AUTOAP 1 // use auto a-priori knowledge + +// QRA code parameters +#define QRA64_K 12 // information symbols +#define QRA64_N 63 // codeword length +#define QRA64_C 51 // (number of parity checks C=(N-K)) +#define QRA64_M 64 // code alphabet size +#define QRA64_m 6 // bits per symbol + +// packed predefined callsigns and fields as defined in JT65 +#define CALL_CQ 0xFA08319 +#define CALL_QRZ 0xFA0831A +#define CALL_CQ000 0xFA0831B +#define CALL_CQ999 0xFA08702 +#define CALL_CQDX 0x5624C39 +#define CALL_DE 0xFF641D1 +#define GRID_BLANK 0x7E91 + +typedef struct { + float decEsNoMetric; + int apflags; + int apmycall; + int apsrccall; + int apmsg_cqqrz[12]; // [cq/qrz ? blank] + int apmsg_call1[12]; // [mycall ? blank] + int apmsg_call1_call2[12]; // [mycall srccall ?] + int apmask_cqqrz[12]; + int apmask_cqqrz_ooo[12]; + int apmask_call1[12]; + int apmask_call1_ooo[12]; + int apmask_call1_call2[12]; +} qra64codec; + +#ifdef __cplusplus +extern "C" { +#endif + +qra64codec *qra64_init(int flags, const int mycall); +// QRA64 mode initialization function +// arguments: +// flags: set the decoder mode +// When flags = QRA_NOAP no a-priori information will be used by the decoder +// When flags = QRA_AUTOAP the decoder will attempt to decode with the amount +// of available a-priori information +// mycall: 28-bit packed callsign of the user (as computed by JT65) +// returns: +// Pointer to the qra64codec data structure allocated and inizialized by the function +// this handle should be passed to the encoding/decoding functions +// +// 0 if unsuccessful (can't allocate memory) +// ------------------------------------------------------------------------------------------- + +void qra64_encode(qra64codec *pcodec, int *y, const int *x); +// QRA64 mode encoder +// arguments: +// pcodec = pointer to a qra64codec data structure as returned by qra64_init +// x = pointer to the message to encode +// x must point to an array of integers (i.e. defined as int x[12]) +// y = pointer to the encoded message +// y must point to an array of integers of lenght 63 (i.e. defined as int y[63]) +// ------------------------------------------------------------------------------------------- + +int qra64_decode(qra64codec *pcodec, int *x, const float *r); +// QRA64 mode decoder +// arguments: +// pcodec = pointer to a qra64codec data structure as returned by qra64_init +// x = pointer to the array of integers where the decoded message will be stored +// x must point to an array of integers (i.e. defined as int x[12]) +// r = pointer to the received symbols energies (squared amplitudes) +// r must point to an array of QRA64_M*QRA64_N (=64*63=4032) float numbers. +// The first QRA_M entries should be the energies of the first symbol in the codeword +// The last QRA_M entries should be the energies of the last symbol in the codeword +// +// return code: +// +// The return code is <0 when decoding is unsuccessful +// -16 indicates that the definition of QRA64_NMSG does not match what required by the code +// If the decoding process is successfull the return code is accordingly to the following table +// rc=0 [? ? ?] AP0 (decoding with no a-priori) +// rc=1 [CQ ? ?] AP27 +// rc=2 [CQ ? ] AP44 +// rc=3 [CALL ? ?] AP29 +// rc=4 [CALL ? ] AP45 +// rc=5 [CALL CALL ?] AP57 +// return codes in the range 1-5 indicate the amount of a-priori information which was required +// to decode the received message and are possible only when the QRA_AUTOAP mode has been enabled. +// ------------------------------------------------------------------------------------------- + +// encode/decode std msgs in 12 symbols as done in jt65 +void encodemsg_jt65(int *y, const int call1, const int call2, const int grid); +void decodemsg_jt65(int *call1, int *call2, int *grid, const int *x); + +#ifdef __cplusplus +} +#endif + +#endif // _qra64_h_ diff --git a/qra64/qra64.vcproj b/qra64/qra64.vcproj new file mode 100644 index 0000000..5a488c9 --- /dev/null +++ b/qra64/qra64.vcproj @@ -0,0 +1,245 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/qra64/qra64_subs.c b/qra64/qra64_subs.c new file mode 100644 index 0000000..933035b --- /dev/null +++ b/qra64/qra64_subs.c @@ -0,0 +1,36 @@ +// qra64_subs.c +// Fortran interface routines for QRA64 + +#include "qra64.h" +#include + +void qra64_enc_(int x[], int y[]) +{ + int ncall=0xf70c238; //K1ABC + qra64codec *codec = qra64_init(0,ncall); //codec for ncall + qra64_encode(codec, y, x); +} + +void qra64_dec_(float r[], int* nmycall, int xdec[], int* rc) +{ +// Return codes: +// rc=-16 failed sanity check +// rc=-2 decoded, but crc check failed +// rc=-1 no decode +// rc=0 [? ? ?] AP0 (decoding with no a-priori information) +// rc=1 [CQ ? ?] AP27 +// rc=2 [CQ ? ] AP42 +// rc=3 [CALL ? ?] AP29 +// rc=4 [CALL ? ] AP44 +// rc=5 [CALL CALL ?] AP57 + + static ncall0=-1; + int ncall=*nmycall; + static qra64codec *codec; + + if(ncall!=ncall0) { + codec = qra64_init(1,ncall); //codec for ncall + ncall0=ncall; + } + *rc = qra64_decode(codec,xdec,r); +} diff --git a/qra64/qra64sim.f90 b/qra64/qra64sim.f90 new file mode 100644 index 0000000..66fa30b --- /dev/null +++ b/qra64/qra64sim.f90 @@ -0,0 +1,243 @@ +program qra64sim + +! Generate simulated QRA64 data for testing the decoder. + + use wavhdr + use packjt + use options + parameter (NMAX=54*12000) ! = 648,000 + parameter (NFFT=10*65536,NH=NFFT/2) + type(hdr) h !Header for .wav file + integer*2 iwave(NMAX) !Generated waveform + integer*4 itone(84) !Channel symbols (values 0-63) + real*4 xnoise(NMAX) !Generated random noise + real*4 dat(NMAX) !Generated real data + complex cdat(NMAX) !Generated complex waveform + complex cspread(0:NFFT-1) !Complex amplitude for Rayleigh fading + complex z + real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq + character msg*22,fname*11,csubmode*1,c,optarg*500,numbuf*32 + character msgsent*22 + logical :: display_help=.false.,seed_prngs=.true. + type (option) :: long_options(9) = [ & + option ('help',.false.,'h','Display this help message',''), & + option ('sub-mode',.true.,'M','message','MSG'), & + option ('sub-mode',.true.,'m','sub mode, default MODE=A','MODE'), & + option ('num-sigs',.true.,'n','number of signals per file, default SIGNALS=10','SIGNALS'), & + option ('doppler-spread',.true.,'d','Doppler spread, default SPREAD=0.0','SPREAD'), & + option ('time-offset',.true.,'t','Time delta, default SECONDS=0.0','SECONDS'), & + option ('num-files',.true.,'f','Number of files to generate, default FILES=1','FILES'), & + option ('no-prng-seed',.false.,'p','Do not seed PRNGs (use for reproducible tests)',''), & + option ('strength',.true.,'s','S/N in dB (2500Hz reference b/w), default SNR=0','SNR') ] + +! Default parameters: + csubmode='A' + mode65=1 + nsigs=10 + fspread=0. + xdt=0. + snrdb=0. + nfiles=1 + + do + call getopt('hM:m:n:d:t:f:ps:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.) + if( nstat .ne. 0 ) then + exit + end if + select case (c) + case ('h') + display_help = .true. + case ('M') + read (optarg(:narglen), *) msg + print*,msg + case ('m') + read (optarg(:narglen), *) csubmode + if(csubmode.eq.'A') mode65=1 + if(csubmode.eq.'B') mode65=2 + if(csubmode.eq.'C') mode65=4 + case ('n') + read (optarg(:narglen), *,err=10) nsigs + case ('d') + read (optarg(:narglen), *,err=10) fspread + case ('t') + read (optarg(:narglen), *) numbuf + if (numbuf(1:1) == '\') then + read (numbuf(2:), *,err=10) xdt + else + read (numbuf, *,err=10) xdt + end if + case ('f') + read (optarg(:narglen), *,err=10) nfiles + case ('p') + seed_prngs=.false. + case ('s') + read (optarg(:narglen), *) numbuf + if (numbuf(1:1) == '\') then + read (numbuf(2:), *,err=10) snrdb + else + read (numbuf, *,err=10) snrdb + end if + end select + cycle +10 display_help=.true. + print *, 'Optional argument format error for option -', c + end do + + if(display_help .or. nstat.lt.0 .or. nremain.ge.1) then + print *, '' + print *, 'Usage: jt65sim [OPTIONS]' + print *, '' + print *, ' Generate one or more simulated JT65 signals in .WAV file(s)' + print *, '' + print *, 'Example: jt65sim -m B -n 10 -d 0.2 -s \\-24.5 -t 0.0 -f 4' + print *, '' + print *, 'OPTIONS: NB Use \ (\\ on *nix shells) to escape -ve arguments' + print *, '' + do i = 1, size (long_options) + call long_options(i) % print (6) + end do + go to 999 + endif + + if (seed_prngs) then + call init_random_seed() !Seed Fortran RANDOM_NUMBER generator + call sgran() !Seed C rand generator (used in gran) + end if + + rms=100. + fsample=12000.d0 !Sample rate (Hz) + dt=1.d0/fsample !Sample interval (s) + twopi=8.d0*atan(1.d0) + npts=54*12000 !Total samples in .wav file + nsps=6912 + baud=12000.d0/nsps !Keying rate = 1.7361111111 + nsym=84 !Number of channel symbols + h=default_header(12000,npts) + dfsig=2000.0/nsigs !Freq spacing between sigs in file (Hz) + ichk=0 + + do ifile=1,nfiles !Loop over requested number of files + write(fname,1002) ifile !Output filename +1002 format('000000_',i4.4) + open(10,file=fname//'.wav',access='stream',status='unknown') + + xnoise=0. + cdat=0. + if(snrdb.lt.90) then + do i=1,npts + xnoise(i)=gran() !Generate gaussian noise + enddo + endif + +! msg="K1ABC W9XYZ EN37" + do isig=1,nsigs !Generate requested number of sigs + if(mod(nsigs,2).eq.0) f0=1500.0 + dfsig*(isig-0.5-nsigs/2) + if(mod(nsigs,2).eq.1) f0=1500.0 + dfsig*(isig-(nsigs+1)/2) + xsnr=snrdb + if(snrdb.eq.0.0) xsnr=-19 - isig + if(csubmode.eq.'B' .and. snrdb.eq.0.0) xsnr=-21 - isig + if(csubmode.eq.'C' .and. snrdb.eq.0.0) xsnr=-21 - isig + + call genqra64(msg,ichk,msgsent,itone,itype) + + bandwidth_ratio=2500.0/6000.0 + sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr) + if(xsnr.gt.90.0) sig=1.0 + write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,msg +1020 format(i4,i4,f10.3,2x,a1,2x,f5.1,f6.2,f5.1,1x,a22) + + phi=0.d0 + dphi=0.d0 + k=12000 + xdt*12000 !Start audio at t = xdt + 1.0 s + isym0=-99 + do i=1,npts !Add this signal into cdat() + isym=i/nsps + 1 + if(isym.gt.nsym) exit + if(isym.ne.isym0) then + freq=f0 + itone(isym)*baud*mode65 + dphi=twopi*freq*dt + isym0=isym + endif + phi=phi + dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + z=cmplx(cos(xphi),sin(xphi)) + k=k+1 + if(k.ge.1) cdat(k)=cdat(k) + sig*z + enddo + enddo + + if(fspread.ne.0) then !Apply specified Doppler spread + df=12000.0/nfft + twopi=8*atan(1.0) + cspread(0)=1.0 + cspread(NH)=0. + +! The following options were added 3/15/2016 to make the half-power tone +! widths equal to the requested Doppler spread. (Previously we effectively +! used b=1.0 and Gaussian shape, which made the tones 1.665 times wider.) +! b=2.0*sqrt(log(2.0)) !Gaussian (before 3/15/2016) +! b=2.0 !Lorenzian 3/15 - 3/27 + b=6.0 !Lorenzian 3/28 onward + + do i=1,NH + f=i*df + x=b*f/fspread + z=0. + a=0. + if(x.lt.3.0) then !Cutoff beyond x=3 +! a=sqrt(exp(-x*x)) !Gaussian + a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian + call random_number(r1) + phi1=twopi*r1 + z=a*cmplx(cos(phi1),sin(phi1)) + endif + cspread(i)=z + z=0. + if(x.lt.50.0) then + call random_number(r2) + phi2=twopi*r2 + z=a*cmplx(cos(phi2),sin(phi2)) + endif + cspread(NFFT-i)=z + enddo + + do i=0,NFFT-1 + f=i*df + if(i.gt.NH) f=(i-nfft)*df + s=real(cspread(i))**2 + aimag(cspread(i))**2 +! write(13,3000) i,f,s,cspread(i) +!3000 format(i5,f10.3,3f12.6) + enddo +! s=real(cspread(0))**2 + aimag(cspread(0))**2 +! write(13,3000) 1024,0.0,s,cspread(0) + + call four2a(cspread,NFFT,1,1,1) !Transform to time domain + + sum=0. + do i=0,NFFT-1 + p=real(cspread(i))**2 + aimag(cspread(i))**2 + sum=sum+p + enddo + avep=sum/NFFT + fac=sqrt(1.0/avep) + cspread=fac*cspread !Normalize to constant avg power + cdat=cspread(1:npts)*cdat !Apply Rayleigh fading + +! do i=0,NFFT-1 +! p=real(cspread(i))**2 + aimag(cspread(i))**2 +! write(14,3010) i,p,cspread(i) +!3010 format(i8,3f12.6) +! enddo + + endif + + dat=aimag(cdat) + xnoise !Add the generated noise + fac=32767.0/nsigs + if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts)) + if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts)) + write(10) h,iwave(1:npts) !Save the .wav file + close(10) + enddo + +999 end program qra64sim diff --git a/qra65.suo b/qra65.suo index 97d8a4f..6dc0656 100644 Binary files a/qra65.suo and b/qra65.suo differ diff --git a/qra65/build.sh b/qra65/build.sh index 44f7045..166ed13 100755 --- a/qra65/build.sh +++ b/qra65/build.sh @@ -1,10 +1,10 @@ -g++ -march=native -pthread -O3 -DFTZ_ENABLE -DQRA_DEBUG -I../qracodes-mt/ \ -*.c \ -../qracodes-mt/normrnd.c \ -../qracodes-mt/npfwht.c \ -../qracodes-mt/pdmath.c \ -../qracodes-mt/qra12_63_64_irr_b.c \ -../qracodes-mt/qra13_64_64_irr_e.c \ -../qracodes-mt/qracodes.c \ +g++ -march=native -pthread -O3 -I../qracodes/ \ +*.c \ +../qracodes/normrnd.c \ +../qracodes/npfwht.c \ +../qracodes/pdmath.c \ +../qracodes/qra12_63_64_irr_b.c \ +../qracodes/qra13_64_64_irr_e.c \ +../qracodes/qracodes.c \ -lpthread -o qra65 diff --git a/qra65/main.c b/qra65/main.c index 982874b..234b687 100644 --- a/qra65/main.c +++ b/qra65/main.c @@ -14,12 +14,12 @@ // main.c - this file // qra65.c/.h - qra65 mode encode/decoding functions // -// ../qracodes-mt/normrnd.c/.h - random gaussian number generator -// ../qracodes-mt/npfwht.c/.h - Fast Walsh-Hadamard Transforms -// ../qracodes-mt/pdmath.c/.h - Elementary math on probability distributions -// ../qracodes-mt/qra12_63_64_irr_b.c/.h - Tables for a QRA(12,63) irregular RA code over GF(64) -// ../qracodes-mt/qra13_64_64_irr_e.c/.h - Tables for a QRA(13,64) irregular RA code " " -// ../qracodes-mt/qracodes.c/.h - QRA codes encoding/decoding functions +// ../qracodes/normrnd.c/.h - random gaussian number generator +// ../qracodes/npfwht.c/.h - Fast Walsh-Hadamard Transforms +// ../qracodes/pdmath.c/.h - Elementary math on probability distributions +// ../qracodes/qra12_63_64_irr_b.c/.h - Tables for a QRA(12,63) irregular RA code over GF(64) +// ../qracodes/qra13_64_64_irr_e.c/.h - Tables for a QRA(13,64) irregular RA code " " +// ../qracodes/qracodes.c/.h - QRA codes encoding/decoding functions // // ------------------------------------------------------------------------------- // @@ -77,7 +77,7 @@ unsigned GetTickCount(void) { #include #include "qra65.h" -#include "../qracodes-mt/normrnd.h" // gaussian numbers generator +#include "../qracodes/normrnd.h" // gaussian numbers generator // ----------------------------------------------------------------------------------- diff --git a/qra65/qra65.vcproj b/qra65/qra65.vcproj index 82d2553..6407280 100644 --- a/qra65/qra65.vcproj +++ b/qra65/qra65.vcproj @@ -179,19 +179,19 @@ > @@ -209,19 +209,19 @@ UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}" > diff --git a/qracodes-mt.suo b/qracodes-mt.suo deleted file mode 100644 index 855a66f..0000000 Binary files a/qracodes-mt.suo and /dev/null differ diff --git a/qracodes-mt/build.sh b/qracodes-mt/build.sh deleted file mode 100755 index 13ebfb1..0000000 --- a/qracodes-mt/build.sh +++ /dev/null @@ -1,2 +0,0 @@ -gcc -Wall -march=native -pthread -O3 -DFTZ_ENABLE *.c -lpthread -lm -o qracodes-mt - diff --git a/qracodes-mt/main.c b/qracodes-mt/main.c deleted file mode 100644 index f518dfd..0000000 --- a/qracodes-mt/main.c +++ /dev/null @@ -1,750 +0,0 @@ -// main.c -// Word Error Rate test example for Q-ary RA codes over GF(64) -// -// (c) 2016 - Nico Palermo, IV3NWV -// -// Thanks to Andrea Montefusco IW0HDV for his help on adapting the sources -// to OSs other than MS Windows -// -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. -// -// Files in this package: -// main.c - this file -// normrnd.c/.h - random gaussian number generator -// npfwht.c/.h - Fast Walsh-Hadamard Transforms -// pdmath.c/.h - Elementary math on probability distributions -// qra12_63_64_irr_b.c/.h - Tables for a QRA(12,63) irregular RA code over GF(64) -// qra13_64_64_irr_e.c/.h - Tables for a QRA(13,64) irregular RA code " " -// qracodes.c/.h - QRA codes encoding/decoding functions -// -// ------------------------------------------------------------------------------- -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -// ----------------------------------------------------------------------------- - -// Two codes are available for simulations in this sowftware release: - -// QRA12_63_64_IRR_B: K=12 N=63 Q=64 irregular QRA code (defined in qra12_63_64_irr_b.h /.c) -// QRA13_64_64_IRR_E: K=13 N=64 Q=64 irregular QRA code (defined in qra13_64_64_irr_b.h /.c) - -// Codes with K=13 are designed to include a CRC as the 13th information symbol -// and improve the code UER (Undetected Error Rate). -// The CRC symbol is not sent along the channel (the codes are punctured) and the -// resulting code is still a (12,63) code with an effective code rate of R = 12/63. - -// ------------------------------------------------------------------------------ - -// OS dependent defines and includes -------------------------------------------- - -#if _WIN32 // note the underscore: without it, it's not msdn official! - // Windows (x64 and x86) - #include // required only for GetTickCount(...) - #include // _beginthread -#endif - -#if defined(__linux__) - -// remove unwanted macros -#define __cdecl - -// implements Windows API -#include - - unsigned int GetTickCount(void) { - struct timespec ts; - unsigned int theTick = 0U; - clock_gettime( CLOCK_REALTIME, &ts ); - theTick = ts.tv_nsec / 1000000; - theTick += ts.tv_sec * 1000; - return theTick; -} - -// Convert Windows millisecond sleep -// -// VOID WINAPI Sleep(_In_ DWORD dwMilliseconds); -// -// to Posix usleep (in microseconds) -// -// int usleep(useconds_t usec); -// -#include -#define Sleep(x) usleep(x*1000) - -#endif - -#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) -#include -#endif - -#if __APPLE__ -#endif - -#include -#include - -#include "qracodes.h" -#include "normrnd.h" // gaussian numbers generator -#include "pdmath.h" // operations on probability distributions - -// defined codes -#include "qra12_63_64_irr_b.h" -#include "qra13_64_64_irr_e.h" - -// ----------------------------------------------------------------------------------- - -#define NTHREADS_MAX 160 - -// channel types -#define CHANNEL_AWGN 0 -#define CHANNEL_RAYLEIGH 1 - -// amount of a-priori information provided to the decoder -#define AP_NONE 0 -#define AP_28 1 -#define AP_44 2 -#define AP_56 3 - -const char ap_str[4][16] = { - "None", - "28 bit", - "44 bit", - "56 bit" -}; - -const char fnameout_pfx[2][64] = { - "wer-awgn-", - "wer-rayleigh-" -}; -const char fnameout_sfx[4][64] = { - "-ap00.txt", - "-ap28.txt", - "-ap44.txt", - "-ap56.txt" -}; - -const int ap_masks_jt65[4][13] = { -// Each row must be 13 entries long (to handle puntc. codes 13,64) -// The mask of 13th symbol (crc) is alway initializated to 0 - // AP0 - no a-priori knowledge - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - // AP28 - 1st field known [cq ? ?] or [dst ? ?] - {0x3F,0x3F,0x3F,0x3F,0x3C, 0, 0, 0, 0, 0, 0, 0}, - // AP44 - 1st and 3rd fields known [cq ? 0] or [dst ? 0] - {0x3F,0x3F,0x3F,0x3F,0x3C, 0, 0, 0, 0,0x0F,0x3F,0x3F}, - // AP56 - 1st and 2nd fields known [dst src ?] - {0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x30, 0, 0} -}; - -void ix_mask(const qracode *pcode, float *r, const int *mask, const int *x); - -void printword(char *msg, int *x, int size) -{ - int k; - printf("\n%s ",msg); - for (k=0;kc msg buffer - float *qra_c2vmsg; //[qra_NMSG*qra_M]; MP decoder c->v msg buffer - float *rp; // [qra_N*qra_M]; received samples (real component) buffer - float *rq; // [qra_N*qra_M]; received samples (imag component) buffer - float *chp; //[qra_N]; channel gains (real component) buffer - float *chq; //[qra_N]; channel gains (imag component) buffer - float *r; //[qra_N*qra_M]; received samples (amplitude) buffer - float *ix; // [qra_N*qra_M]; // intrinsic information to the MP algorithm - float *ex; // [qra_N*qra_M]; // extrinsic information from the MP algorithm -#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) - pthread_t thread; -#endif -} wer_test_ds; - -typedef void( __cdecl *pwer_test_thread)(wer_test_ds*); - -// crc-6 generator polynomial -// g(x) = x^6 + a5*x^5 + ... + a1*x + a0 - -// g(x) = x^6 + x + 1 -#define CRC6_GEN_POL 0x30 // MSB=a0 LSB=a5 - -// g(x) = x^6 + x^2 + x + 1 (as suggested by Joe. See: https://users.ece.cmu.edu/~koopman/crc/) -// #define CRC6_GEN_POL 0x38 // MSB=a0 LSB=a5. Simulation results are similar - -int calc_crc6(int *x, int sz) -{ - int k,j,t,sr = 0; - for (k=0;k>1) ^ CRC6_GEN_POL; - else - sr = (sr>>1); - t>>=1; - } - } - return sr; -} - -void wer_test_thread(wer_test_ds *pdata) -{ - const qracode *pcode=pdata->pcode; - const int qra_K = pcode->K; - const int qra_N = pcode->N; - const int qra_M = pcode->M; - const int qra_m = pcode->m; - const int NSAMPLES = pcode->N*pcode->M; - - const float No = 1.0f; // noise spectral density - const float sigma = (float)sqrt(No/2.0f); // std dev of noise I/Q components - const float sigmach = (float)sqrt(1/2.0f); // std dev of channel I/Q gains - - // Eb/No value for which we optimize the bessel metric - const float EbNodBMetric = 2.8f; - const float EbNoMetric = (float)pow(10,EbNodBMetric/10); - - int k,t,diff; - - float R; - float EsNoMetric; - float EbNo, EsNo, Es, A; - int channel_type, code_type; - int nt=0; // transmitted codewords - int nerrs = 0; // total number of errors - int nerrsu = 0; // number of undetected errors - int rc; - - - // inizialize pointer to required buffers - int *x=pdata->x; // message buffer - int *y=pdata->y, *ydec=pdata->ydec; // encoded/decoded codeword buffers - float *qra_v2cmsg=pdata->qra_v2cmsg; // table of the v->c messages - float *qra_c2vmsg=pdata->qra_c2vmsg; // table of the c->v messages - float *rp=pdata->rp; // received samples (real component) - float *rq=pdata->rq; // received samples (imag component) - float *chp=pdata->chp; // channel gains (real component) - float *chq=pdata->chq; // channel gains (imag component) - float *r=pdata->r; // received samples amplitudes - float *ix=pdata->ix; // intrinsic information to the MP algorithm - float *ex=pdata->ex; // extrinsic information from the MP algorithm - - channel_type = pdata->channel_type; - code_type = pcode->type; - - // define the (true) code rate accordingly to the code type - switch(code_type) { - case QRATYPE_CRC: - R = 1.0f*(qra_K-1)/qra_N; - break; - case QRATYPE_CRCPUNCTURED: - R = 1.0f*(qra_K-1)/(qra_N-1); - break; - case QRATYPE_NORMAL: - default: - R = 1.0f*(qra_K)/(qra_N); - } - - EsNoMetric = 1.0f*qra_m*R*EbNoMetric; - - EbNo = (float)pow(10,pdata->EbNodB/10); - EsNo = 1.0f*qra_m*R*EbNo; - Es = EsNo*No; - A = (float)sqrt(Es); - - - // encode the input - if (code_type==QRATYPE_CRC || code_type==QRATYPE_CRCPUNCTURED) { - // compute the information message symbol check as the (negated) xor of all the - // information message symbols - for (k=0;k<(qra_K-1);k++) - x[k]=k%qra_M; - x[k]=calc_crc6(x,qra_K-1); - } - else - for (k=0;kstop==0) { - - // simulate the channel - // NOTE: in the case that the code is punctured, for simplicity - // we compute the channel outputs and the metric also for the crc symbol - // then we ignore its observation. - normrnd_s(rp,NSAMPLES,0,sigma); - normrnd_s(rq,NSAMPLES,0,sigma); - - if (channel_type == CHANNEL_AWGN) { - for (k=0;kdone = 1; - return; // unknown channel type - } - - // compute the squares of the amplitudes of the received samples - for (k=0;km,pcode->N,EsNoMetric); - - if (code_type==QRATYPE_CRCPUNCTURED) { - // ignore observations of the CRC symbol as it is not actually sent - // over the channel - pd_init(PD_ROWADDR(ix,qra_M,qra_K),pd_uniform(qra_m),qra_M); - } - - - if (pdata->ap_index!=0) - // mask channel observations with a priori knowledge - ix_mask(pcode,ix,ap_masks_jt65[pdata->ap_index],x); - - - // compute the extrinsic symbols probabilities with the message-passing algorithm - // stop if extrinsic information does not converges to 1 within the given number of iterations - rc = qra_extrinsic(pcode,ex,ix,100,qra_v2cmsg,qra_c2vmsg); - - if (rc>=0) { // the MP algorithm converged to Iex~1 in rc iterations - - // decode the codeword - qra_mapdecode(pcode,ydec,ex,ix); - - // look for undetected errors - if (code_type==QRATYPE_CRC || code_type==QRATYPE_CRCPUNCTURED) { - - diff = 0; - for (k=0;k<(qra_K-1);k++) - diff |= (ydec[k]!=x[k]); - t = calc_crc6(ydec,qra_K-1); - if (t!=ydec[k]) // error detected - crc doesn't matches - nerrs += 1; - else - if (diff) { // decoded message is not equal to the transmitted one but - // the crc test passed - // add as undetected error - nerrsu += 1; - nerrs += 1; - // uncomment to see what the undetected error pattern looks like - //printword("U", ydec); - } - } - else - for (k=0;knt=nt; - pdata->nerrs=nerrs; - pdata->nerrsu=nerrsu; - - } - - pdata->done=1; - - #if _WIN32 - _endthread(); - #endif -} - -#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) - -void *wer_test_pthread(void *p) -{ - wer_test_thread ((wer_test_ds *)p); - return 0; -} - -#endif - - -void ix_mask(const qracode *pcode, float *r, const int *mask, const int *x) -{ - // mask intrinsic information (channel observations) with a priori knowledge - - int k,kk, smask; - const int qra_K=pcode->K; - const int qra_M=pcode->M; - const int qra_m=pcode->m; - - for (k=0;kNTHREADS_MAX) { - printf("Error: nthreads should be <=%d\n",NTHREADS_MAX); - return -1; - } - - sprintf(fnameout,"%s%s%s", - fnameout_pfx[chtype], - pcode->name, - fnameout_sfx[ap_index]); - - fout = fopen(fnameout,"w"); - fprintf(fout,"# Channel (0=AWGN,1=Rayleigh), Eb/No (dB), Transmitted codewords, Errors, Undetected Errors, Avg dec. time (ms), WER\n"); - - printf("\nTesting the code %s over the %s channel\nSimulation data will be saved to %s\n", - pcode->name, - chtype==CHANNEL_AWGN?"AWGN":"Rayleigh", - fnameout); - fflush (stdout); - - // init fixed thread parameters and preallocate buffers - for (j=0;jK*sizeof(int)); - wt[j].y = (int*)malloc(pcode->N*sizeof(int)); - wt[j].ydec = (int*)malloc(pcode->N*sizeof(int)); - wt[j].qra_v2cmsg = (float*)malloc(pcode->NMSG*pcode->M*sizeof(float)); - wt[j].qra_c2vmsg = (float*)malloc(pcode->NMSG*pcode->M*sizeof(float)); - wt[j].rp = (float*)malloc(pcode->N*pcode->M*sizeof(float)); - wt[j].rq = (float*)malloc(pcode->N*pcode->M*sizeof(float)); - wt[j].chp = (float*)malloc(pcode->N*sizeof(float)); - wt[j].chq = (float*)malloc(pcode->N*sizeof(float)); - wt[j].r = (float*)malloc(pcode->N*pcode->M*sizeof(float)); - wt[j].ix = (float*)malloc(pcode->N*pcode->M*sizeof(float)); - wt[j].ex = (float*)malloc(pcode->N*pcode->M*sizeof(float)); - } - - - for (k=0;k=nerrstgt[k]) { - for (j=0;j] [-t] [-c] [-a] [-f[-h]\n"); - printf("Options: \n"); - printf(" -q: code to simulate. 0=qra_12_63_64_irr_b\n"); - printf(" 1=qra_13_64_64_irr_e (default)\n"); - printf(" -t : number of threads to be used for the simulation [1..24]\n"); - printf(" (default=8)\n"); - printf(" -c : channel_type. 0=AWGN 1=Rayleigh \n"); - printf(" (default=AWGN)\n"); - printf(" -a : amount of a-priori information provided to decoder. \n"); - printf(" 0= No a-priori (default)\n"); - printf(" 1= 28 bit \n"); - printf(" 2= 44 bit \n"); - printf(" 3= 56 bit \n"); - printf(" -f : name of the file containing the Eb/No values to be simulated\n"); - printf(" (default=ebnovalues.txt)\n"); - printf(" This file should contain lines in this format:\n"); - printf(" # Eb/No(dB) Target Errors\n"); - printf(" 0.1 5000\n"); - printf(" 0.6 5000\n"); - printf(" 1.1 1000\n"); - printf(" 1.6 1000\n"); - printf(" ...\n"); - printf(" (lines beginning with a # are treated as comments\n\n"); - - printf(" sizeof(unsigned int) = %lu bytes\n", (unsigned long) sizeof(unsigned int)); - printf(" sizeof(unsigned long) = %lu bytes\n", (unsigned long) sizeof(unsigned long)); - printf(" sizeof(unsigned long long) = %lu bytes\n", (unsigned long) sizeof(unsigned long long)); - printf(" sizeof(unsigned float) = %lu bytes\n", (unsigned long) sizeof(float)); - printf(" sizeof(unsigned double) = %lu bytes\n", (unsigned long) sizeof(double)); - printf(" sizeof(void *) = %lu bytes\n", (unsigned long) sizeof(void *)); - - printf("\n\n"); -} - -#define SIM_POINTS_MAX 20 - -int main(int argc, char* argv[]) -{ - - float EbNodB[SIM_POINTS_MAX]; - int nerrstgt[SIM_POINTS_MAX]; - FILE *fin; - - char fnamein[128]= "ebnovalues.txt"; - char buf[128]; - - int nitems = 0; - int code_idx = 1; - int nthreads = 8; - int ch_type = CHANNEL_AWGN; - int ap_index = AP_NONE; - - // parse command line - while(--argc) { - argv++; - if (strncmp(*argv,"-h",2)==0) { - syntax(); - return 0; - } - else - if (strncmp(*argv,"-q",2)==0) { - code_idx = (int)atoi((*argv)+2); - if (code_idx>1) { - printf("Invalid code index\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-t",2)==0) { - nthreads = (int)atoi((*argv)+2); - - printf("nthreads = %d\n",nthreads); - - if (nthreads>NTHREADS_MAX) { - printf("Invalid number of threads\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-c",2)==0) { - ch_type = (int)atoi((*argv)+2); - if (ch_type>CHANNEL_RAYLEIGH) { - printf("Invalid channel type\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-a",2)==0) { - ap_index = (int)atoi((*argv)+2); - if (ap_index>AP_56) { - printf("Invalid a-priori information index\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-f",2)==0) { - strncpy(fnamein,(*argv)+2,127); - } - else - if (strncmp(*argv,"-h",2)==0) { - syntax(); - return -1; - } - else { - printf("Invalid option\n"); - syntax(); - return -1; - } - } - - // parse points to be simulated from the input file - fin = fopen(fnamein,"r"); - if (!fin) { - printf("Can't open file: %s\n",fnamein); - syntax(); - } - - while (fgets(buf,128,fin)!=0) - if (*buf=='#' || *buf=='\n' ) - continue; - else - if (nitems==SIM_POINTS_MAX) - break; - else - if (sscanf(buf,"%f %u",&EbNodB[nitems],&nerrstgt[nitems])!=2) { - printf("Invalid input file format\n"); - syntax(); - return -1; - } - else - nitems++; - - fclose(fin); - - if (nitems==0) { - printf("No Eb/No point specified in file %s\n",fnamein); - syntax(); - return -1; - } - - printf("\nQ-ary Repeat-Accumulate Code Word Error Rate Simulator\n"); - printf("2016, Nico Palermo - IV3NWV\n\n"); - - printf("Nthreads = %d\n",nthreads); - printf("Channel = %s\n",ch_type==CHANNEL_AWGN?"AWGN":"Rayleigh"); - printf("Codename = %s\n",codetotest[code_idx]->name); - printf("A-priori = %s\n",ap_str[ap_index]); - printf("Eb/No input file = %s\n\n",fnamein); - - wer_test_proc(codetotest[code_idx], nthreads, ch_type, ap_index, EbNodB, nerrstgt, nitems); - - return 0; -} - diff --git a/qracodes-mt/qracodes.h b/qracodes-mt/qracodes.h deleted file mode 100644 index fb74b2e..0000000 --- a/qracodes-mt/qracodes.h +++ /dev/null @@ -1,79 +0,0 @@ -// qracodes.h -// Q-ary RA codes encoding/decoding functions -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -#ifndef _qracodes_h_ -#define _qracodes_h_ - -// type of codes -#define QRATYPE_NORMAL 0x00 // normal code -#define QRATYPE_CRC 0x01 // code with crc - last information symbol is a CRC -#define QRATYPE_CRCPUNCTURED 0x02 // the CRC symbol is punctured (not sent along the channel) - - -typedef struct { - // code parameters - const int K; // number of information symbols - const int N; // codeword length in symbols - const int m; // bits/symbol - const int M; // Symbol alphabet cardinality (2^m) - const int a; // code grouping factor - const int NC; // number of check symbols (N-K) - const int V; // number of variables in the code graph (N) - const int C; // number of factors in the code graph (N +(N-K)+1) - const int NMSG; // number of msgs in the code graph - const int MAXVDEG; // maximum variable degree - const int MAXCDEG; // maximum factor degree - const int type; // see QRATYPE_xx defines - const float R; // code rate (K/N) - const char name[64]; // code name - // tables used by the encoder - const int *acc_input_idx; - const int *acc_input_wlog; - const int *gflog; - const int *gfexp; - // tables used by the decoder ------------------------- - const int *msgw; - const int *vdeg; - const int *cdeg; - const int *v2cmidx; - const int *c2vmidx; - const int *gfpmat; -} qracode; -// Uncomment the header file of the code which needs to be tested - -//#include "qra12_63_64_irr_b.h" // irregular code (12,63) over GF(64) -//#include "qra13_64_64_irr_e.h" // irregular code with good performance and best UER protection at AP56 -//#include "qra13_64_64_reg_a.h" // regular code with good UER but perf. inferior to that of code qra12_63_64_irr_b - -#ifdef __cplusplus -extern "C" { -#endif - -int qra_encode(const qracode *pcode, int *y, const int *x); -void qra_mfskbesselmetric(float *pix, const float *rsq, const int m, const int N, float EsNoMetric); -int qra_extrinsic(const qracode *pcode, float *pex, const float *pix, int maxiter,float *qra_v2cmsg,float *qra_c2vmsg); -void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix); - -#ifdef __cplusplus -} -#endif - -#endif // _qracodes_h_ diff --git a/qracodes-st.sln b/qracodes-st.sln new file mode 100644 index 0000000..a90eb89 --- /dev/null +++ b/qracodes-st.sln @@ -0,0 +1,20 @@ + +Microsoft Visual Studio Solution File, Format Version 10.00 +# Visual Studio 2008 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "qracodes-st", "qracodes-st\qracodes-st.vcproj", "{9E921A72-CE23-46C1-88A6-AC4D5BF5F091}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Win32 = Debug|Win32 + Release|Win32 = Release|Win32 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {9E921A72-CE23-46C1-88A6-AC4D5BF5F091}.Debug|Win32.ActiveCfg = Debug|Win32 + {9E921A72-CE23-46C1-88A6-AC4D5BF5F091}.Debug|Win32.Build.0 = Debug|Win32 + {9E921A72-CE23-46C1-88A6-AC4D5BF5F091}.Release|Win32.ActiveCfg = Release|Win32 + {9E921A72-CE23-46C1-88A6-AC4D5BF5F091}.Release|Win32.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/qracodes-st.suo b/qracodes-st.suo new file mode 100644 index 0000000..ab5b3d8 Binary files /dev/null and b/qracodes-st.suo differ diff --git a/qracodes-st/build.sh b/qracodes-st/build.sh new file mode 100644 index 0000000..f8dceb5 --- /dev/null +++ b/qracodes-st/build.sh @@ -0,0 +1,3 @@ +g++ -march=native -O3 -DFTZ_ENABLE -DQRA_DEBUG -DTEST_WER_AWGN *.c -o qracodesawgn && ./qracodesawgn +g++ -march=native -O3 -DFTZ_ENABLE -DQRA_DEBUG -DTEST_WER_RAYLEIGH *.c -o qracodesrayleigh && ./qracodesrayleigh + diff --git a/qracodes-st/main.c b/qracodes-st/main.c new file mode 100644 index 0000000..200a28f --- /dev/null +++ b/qracodes-st/main.c @@ -0,0 +1,270 @@ +// main.c +// Word Error Rate test for the Q-ary RA code (12,63) over GF(64) +// First, single threaded simulator. +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// Files in the package: +// main.c - this file +// normrnd.c/.h - random gaussian number generator and header +// npfwht.c/.h - Fast Walsh-Hadamard Transforms +// pdmath.c/.h - Elementary math on probability distributions +// qra12_63_64_irr_b.c/.h - Tables and defines for a P(12,63) IRA code over GF(64) +// qracodes.c/.h - QRA (Q-ary Repeat & Accumulate) codes encoding/decoding functions +// +// ------------------------------------------------------------------------------- +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#if _WIN32 // note the underscore: without it, it's not msdn official! + // Windows (x64 and x86) + #include // required only for GetTickCount(...) +#endif + +#if __linux__ +#include +#include + +unsigned GetTickCount() { + struct timespec ts; + unsigned theTick = 0U; + clock_gettime( CLOCK_REALTIME, &ts ); + theTick = ts.tv_nsec / 1000000; + theTick += ts.tv_sec * 1000; + return theTick; +} +#endif + +#if __APPLE__ +#endif + +#include + +#ifdef FTZ_ENABLE // shouldn't be required anymore - kept here just as a remind in the case we need it + // round off errors in the qra_extrinsics fixed with biasing +#include // used to set the Flush to Zero flag in the FPU +#endif + +#include "normrnd.h" // gaussian number generators +#include "qra12_63_64_irr_b.h" // tables and defines for the P(12,63) Q-ary RA code over GF(64) +#include "qracodes.h" // encoding/decoding functions + +#define NSAMPLES (qra_N*qra_M) // samples per message (codeword length * codeword alphabet size) + +#define CHANNEL_AWGN 0 +#define CHANNEL_RAYLEIGH 1 + +int wer_test(int channel_type, float EbNodB, int desirederrs, FILE *fout) +{ + int k; + +// uint x[qra_K]= {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}; + uint x[qra_K]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,}; + uint y[qra_N],ydec[qra_N]; + + float rp[qra_N*qra_M]; // received samples (real component) + float rq[qra_N*qra_M]; // received samples (imag component) + float chp[qra_N]; // channel gains (real component) + float chq[qra_N]; // channel gains (imag component) + float r[qra_N*qra_M]; // r = sqrt(rp^2+rq^2) + + float No = 1.0f; // noise spectral density + float sigma = (float)sqrt(No/2.0f); // std dev of noise I/Q components + float sigmach = (float)sqrt(1/2.0f); // std dev of channel I/Q gains + float R = 1.0f*qra_K/qra_N; // code rate + + float EbNo = (float)pow(10,EbNodB/10); + float EsNo = 1.0f*qra_m*R*EbNo; + float Es = EsNo*No; + float A = (float)sqrt(Es); + + // Eb/No value for which we optimize the bessel metric + const float EbNodBMetric = 2.8f; + const float EbNoMetric = (float)pow(10,EbNodBMetric/10); + const float EsNoMetric = 1.0f*qra_m*R*EbNoMetric; + + float ex[qra_N*qra_M]; // extrinsic information from the MP algorithm + float avgt; // average time per codeword encoding/channel sim/decoding + int nt=0; // transmitted codewords + int nerrs = 0; // total number of errors + int nerrsu = 0; // number of undetected errors + int nep = -1; + + unsigned int cini,cend; // time tick counters + + int rc; + + // encode the input + qra_encode(y,x); + + + cini = GetTickCount(); + + while (1) { + + // simulate the channel + normrnd_s(rp,NSAMPLES,0,sigma); + normrnd_s(rq,NSAMPLES,0,sigma); + + if (channel_type == CHANNEL_AWGN) { + for (k=0;k=0) { // the MP algorithm converged to Iex~1 in rc iterations + + // decode the codeword + qra_mapdecode(ydec,ex,r); + + // look for undetected errors + for (k=0;k. - - -#include "normrnd.h" - -#if _WIN32 // note the underscore: without it, it's not msdn official! - // Windows (x64 and x86) - #include // required only for GetTickCount(...) - #define K_RAND_MAX UINT_MAX -#elif __unix__ // all unices, not all compilers - #include - #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF - #define K_RAND_MAX 0x7FFFFFFF // that's the max number generated by lrand48 - // Unix -#elif __linux__ - #include - #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF - #define K_RAND_MAX 0x7FFFFFFF // that's the max number generated by lrand48 - // linux -#elif __APPLE__ - // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though... -#endif - - -// use MS rand_s(...) function -void normrnd_s(float *dst, int nitems, float mean, float stdev) -{ - unsigned int r; - float phi=0, u=0; - int set = 0; - - while (nitems--) - if (set==1) { - *dst++ = (float)sin(phi)*u*stdev+mean; - set = 0; - } - else { - rand_s((unsigned int*)&r); phi = (M_2PI/(1.0f+K_RAND_MAX))*r; - rand_s((unsigned int*)&r); u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+K_RAND_MAX))*(1.0f+r) ) ); - *dst++ = (float)cos(phi)*u*stdev+mean; - set=1; - } -} - -/* NOT USED -// use MS rand() function -void normrnd(float *dst, int nitems, float mean, float stdev) -{ - float phi=0, u=0; - int set = 0; - - while (nitems--) - if (set==1) { - *dst++ = (float)sin(phi)*u*stdev+mean; - set = 0; - } - else { - phi = (M_2PI/(1.0f+RAND_MAX))*rand(); - u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+RAND_MAX))*(1.0f+rand()) ) ); - *dst++ = (float)cos(phi)*u*stdev+mean; - set=1; - } -} -*/ +// normrnd.c +// functions to generate gaussian distributed numbers +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + + +#include "normrnd.h" + +#if _WIN32 // note the underscore: without it, it's not msdn official! + // Windows (x64 and x86) + #include // required only for GetTickCount(...) + #define K_RAND_MAX UINT_MAX +#elif __unix__ // all unices, not all compilers + #include + #include + #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF + #define K_RAND_MAX 0x7FFFFFFF // that's the max number generated by lrand48 + // Unix +#elif __linux__ + #include + #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF + #define K_RAND_MAX 0x7FFFFFFF // that's the max number generated by lrand48 + // linux +#elif __APPLE__ + // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though... +#endif + +// use MS rand_s(...) function +void normrnd_s(float *dst, int nitems, float mean, float stdev) +{ + unsigned int r; + float phi, u; + int set = 0; + + while (nitems--) + if (set==1) { + *dst++ = (float)sin(phi)*u*stdev+mean; + set = 0; + } + else { + rand_s(&r); phi = (M_2PI/(1.0f+K_RAND_MAX))*r; + rand_s(&r); u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+K_RAND_MAX))*(1.0f+r) ) ); + *dst++ = (float)cos(phi)*u*stdev+mean; + set=1; + } +} + + +// use MS rand() function +void normrnd(float *dst, int nitems, float mean, float stdev) +{ + float phi, u; + int set = 0; + + while (nitems--) + if (set==1) { + *dst++ = (float)sin(phi)*u*stdev+mean; + set = 0; + } + else { + phi = (M_2PI/(1.0f+RAND_MAX))*rand(); + u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+RAND_MAX))*(1.0f+rand()) ) ); + *dst++ = (float)cos(phi)*u*stdev+mean; + set=1; + } +} diff --git a/qracodes-mt/normrnd.h b/qracodes-st/normrnd.h similarity index 93% rename from qracodes-mt/normrnd.h rename to qracodes-st/normrnd.h index 04bbf01..8471821 100644 --- a/qracodes-mt/normrnd.h +++ b/qracodes-st/normrnd.h @@ -1,51 +1,49 @@ -// normrnd.h -// Functions to generate gaussian distributed numbers -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -#ifndef _normrnd_h_ -#define _normrnd_h_ - -#define _CRT_RAND_S -#include - -#define _USE_MATH_DEFINES -#include -#define M_2PI (2.0f*(float)M_PI) - -#ifdef __cplusplus -extern "C" { -#endif - -void normrnd_s(float *dst, int nitems, float mean, float stdev); -// generate a random array of numbers with a gaussian distribution of given mean and stdev -// use MS rand_s(...) function - -/* not used -void normrnd(float *dst, int nitems, float mean, float stdev); -// generate a random array of numbers with a gaussian distribution of given mean and stdev -// use MS rand() function -*/ - -#ifdef __cplusplus -} -#endif - -#endif // _normrnd_h_ - +// normrnd.h +// Functions to generate gaussian distributed numbers +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#ifndef _normrnd_h_ +#define _normrnd_h_ + +#define _CRT_RAND_S +#include + +#define _USE_MATH_DEFINES +#include +#define M_2PI (2.0f*(float)M_PI) + +#ifdef __cplusplus +extern "C" { +#endif + +void normrnd_s(float *dst, int nitems, float mean, float stdev); +// generate a random array of numbers with a gaussian distribution of given mean and stdev +// use MS rand_s(...) function + +void normrnd(float *dst, int nitems, float mean, float stdev); +// generate a random array of numbers with a gaussian distribution of given mean and stdev +// use MS rand() function + +#ifdef __cplusplus +} +#endif + +#endif // _normrnd_h_ + diff --git a/qracodes-mt/npfwht.c b/qracodes-st/npfwht.c similarity index 97% rename from qracodes-mt/npfwht.c rename to qracodes-st/npfwht.c index 0f8fa42..70df5ab 100644 --- a/qracodes-mt/npfwht.c +++ b/qracodes-st/npfwht.c @@ -1,216 +1,216 @@ -// npfwht.c -// Basic implementation of the Fast Walsh-Hadamard Transforms -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -#include "npfwht.h" - -#define WHBFY(dst,src,base,offs,dist) { dst[base+offs]=src[base+offs]+src[base+offs+dist]; dst[base+offs+dist]=src[base+offs]-src[base+offs+dist]; } - -typedef void (*pnp_fwht)(float*,float*); - -static void np_fwht2(float *dst, float *src); - -static void np_fwht1(float *dst, float *src); -static void np_fwht2(float *dst, float *src); -static void np_fwht4(float *dst, float *src); -static void np_fwht8(float *dst, float *src); -static void np_fwht16(float *dst, float *src); -static void np_fwht32(float *dst, float *src); -static void np_fwht64(float *dst, float *src); - -static pnp_fwht np_fwht_tab[7] = { - np_fwht1, - np_fwht2, - np_fwht4, - np_fwht8, - np_fwht16, - np_fwht32, - np_fwht64 -}; - -void np_fwht(int nlogdim, float *dst, float *src) -{ - np_fwht_tab[nlogdim](dst,src); -} - -static void np_fwht1(float *dst, float *src) -{ - dst[0] = src[0]; -} - - -static void np_fwht2(float *dst, float *src) -{ - float t[2]; - - WHBFY(t,src,0,0,1); - dst[0]= t[0]; - dst[1]= t[1]; -} - -static void np_fwht4(float *dst, float *src) -{ - float t[4]; - - // group 1 - WHBFY(t,src,0,0,2); WHBFY(t,src,0,1,2); - // group 2 - WHBFY(dst,t,0,0,1); WHBFY(dst,t,2,0,1); -}; - - -static void np_fwht8(float *dst, float *src) -{ - float t[16]; - float *t1=t, *t2=t+8; - - // group 1 - WHBFY(t1,src,0,0,4); WHBFY(t1,src,0,1,4); WHBFY(t1,src,0,2,4); WHBFY(t1,src,0,3,4); - // group 2 - WHBFY(t2,t1,0,0,2); WHBFY(t2,t1,0,1,2); WHBFY(t2,t1,4,0,2); WHBFY(t2,t1,4,1,2); - // group 3 - WHBFY(dst,t2,0,0,1); WHBFY(dst,t2,2,0,1); WHBFY(dst,t2,4,0,1); WHBFY(dst,t2,6,0,1); -}; - - -static void np_fwht16(float *dst, float *src) -{ - float t[32]; - float *t1=t, *t2=t+16; - - // group 1 - WHBFY(t1,src,0,0,8); WHBFY(t1,src,0,1,8); WHBFY(t1,src,0,2,8); WHBFY(t1,src,0,3,8); - WHBFY(t1,src,0,4,8); WHBFY(t1,src,0,5,8); WHBFY(t1,src,0,6,8); WHBFY(t1,src,0,7,8); - // group 2 - WHBFY(t2,t1,0,0,4); WHBFY(t2,t1,0,1,4); WHBFY(t2,t1,0,2,4); WHBFY(t2,t1,0,3,4); - WHBFY(t2,t1,8,0,4); WHBFY(t2,t1,8,1,4); WHBFY(t2,t1,8,2,4); WHBFY(t2,t1,8,3,4); - // group 3 - WHBFY(t1,t2,0,0,2); WHBFY(t1,t2,0,1,2); WHBFY(t1,t2,4,0,2); WHBFY(t1,t2,4,1,2); - WHBFY(t1,t2,8,0,2); WHBFY(t1,t2,8,1,2); WHBFY(t1,t2,12,0,2); WHBFY(t1,t2,12,1,2); - // group 4 - WHBFY(dst,t1,0,0,1); WHBFY(dst,t1,2,0,1); WHBFY(dst,t1,4,0,1); WHBFY(dst,t1,6,0,1); - WHBFY(dst,t1,8,0,1); WHBFY(dst,t1,10,0,1); WHBFY(dst,t1,12,0,1); WHBFY(dst,t1,14,0,1); - -} - -static void np_fwht32(float *dst, float *src) -{ - float t[64]; - float *t1=t, *t2=t+32; - - // group 1 - WHBFY(t1,src,0,0,16); WHBFY(t1,src,0,1,16); WHBFY(t1,src,0,2,16); WHBFY(t1,src,0,3,16); - WHBFY(t1,src,0,4,16); WHBFY(t1,src,0,5,16); WHBFY(t1,src,0,6,16); WHBFY(t1,src,0,7,16); - WHBFY(t1,src,0,8,16); WHBFY(t1,src,0,9,16); WHBFY(t1,src,0,10,16); WHBFY(t1,src,0,11,16); - WHBFY(t1,src,0,12,16); WHBFY(t1,src,0,13,16); WHBFY(t1,src,0,14,16); WHBFY(t1,src,0,15,16); - - // group 2 - WHBFY(t2,t1,0,0,8); WHBFY(t2,t1,0,1,8); WHBFY(t2,t1,0,2,8); WHBFY(t2,t1,0,3,8); - WHBFY(t2,t1,0,4,8); WHBFY(t2,t1,0,5,8); WHBFY(t2,t1,0,6,8); WHBFY(t2,t1,0,7,8); - WHBFY(t2,t1,16,0,8); WHBFY(t2,t1,16,1,8); WHBFY(t2,t1,16,2,8); WHBFY(t2,t1,16,3,8); - WHBFY(t2,t1,16,4,8); WHBFY(t2,t1,16,5,8); WHBFY(t2,t1,16,6,8); WHBFY(t2,t1,16,7,8); - - // group 3 - WHBFY(t1,t2,0,0,4); WHBFY(t1,t2,0,1,4); WHBFY(t1,t2,0,2,4); WHBFY(t1,t2,0,3,4); - WHBFY(t1,t2,8,0,4); WHBFY(t1,t2,8,1,4); WHBFY(t1,t2,8,2,4); WHBFY(t1,t2,8,3,4); - WHBFY(t1,t2,16,0,4); WHBFY(t1,t2,16,1,4); WHBFY(t1,t2,16,2,4); WHBFY(t1,t2,16,3,4); - WHBFY(t1,t2,24,0,4); WHBFY(t1,t2,24,1,4); WHBFY(t1,t2,24,2,4); WHBFY(t1,t2,24,3,4); - - // group 4 - WHBFY(t2,t1,0,0,2); WHBFY(t2,t1,0,1,2); WHBFY(t2,t1,4,0,2); WHBFY(t2,t1,4,1,2); - WHBFY(t2,t1,8,0,2); WHBFY(t2,t1,8,1,2); WHBFY(t2,t1,12,0,2); WHBFY(t2,t1,12,1,2); - WHBFY(t2,t1,16,0,2); WHBFY(t2,t1,16,1,2); WHBFY(t2,t1,20,0,2); WHBFY(t2,t1,20,1,2); - WHBFY(t2,t1,24,0,2); WHBFY(t2,t1,24,1,2); WHBFY(t2,t1,28,0,2); WHBFY(t2,t1,28,1,2); - - // group 5 - WHBFY(dst,t2,0,0,1); WHBFY(dst,t2,2,0,1); WHBFY(dst,t2,4,0,1); WHBFY(dst,t2,6,0,1); - WHBFY(dst,t2,8,0,1); WHBFY(dst,t2,10,0,1); WHBFY(dst,t2,12,0,1); WHBFY(dst,t2,14,0,1); - WHBFY(dst,t2,16,0,1); WHBFY(dst,t2,18,0,1); WHBFY(dst,t2,20,0,1); WHBFY(dst,t2,22,0,1); - WHBFY(dst,t2,24,0,1); WHBFY(dst,t2,26,0,1); WHBFY(dst,t2,28,0,1); WHBFY(dst,t2,30,0,1); - -} - -static void np_fwht64(float *dst, float *src) -{ - float t[128]; - float *t1=t, *t2=t+64; - - - // group 1 - WHBFY(t1,src,0,0,32); WHBFY(t1,src,0,1,32); WHBFY(t1,src,0,2,32); WHBFY(t1,src,0,3,32); - WHBFY(t1,src,0,4,32); WHBFY(t1,src,0,5,32); WHBFY(t1,src,0,6,32); WHBFY(t1,src,0,7,32); - WHBFY(t1,src,0,8,32); WHBFY(t1,src,0,9,32); WHBFY(t1,src,0,10,32); WHBFY(t1,src,0,11,32); - WHBFY(t1,src,0,12,32); WHBFY(t1,src,0,13,32); WHBFY(t1,src,0,14,32); WHBFY(t1,src,0,15,32); - WHBFY(t1,src,0,16,32); WHBFY(t1,src,0,17,32); WHBFY(t1,src,0,18,32); WHBFY(t1,src,0,19,32); - WHBFY(t1,src,0,20,32); WHBFY(t1,src,0,21,32); WHBFY(t1,src,0,22,32); WHBFY(t1,src,0,23,32); - WHBFY(t1,src,0,24,32); WHBFY(t1,src,0,25,32); WHBFY(t1,src,0,26,32); WHBFY(t1,src,0,27,32); - WHBFY(t1,src,0,28,32); WHBFY(t1,src,0,29,32); WHBFY(t1,src,0,30,32); WHBFY(t1,src,0,31,32); - - // group 2 - WHBFY(t2,t1,0,0,16); WHBFY(t2,t1,0,1,16); WHBFY(t2,t1,0,2,16); WHBFY(t2,t1,0,3,16); - WHBFY(t2,t1,0,4,16); WHBFY(t2,t1,0,5,16); WHBFY(t2,t1,0,6,16); WHBFY(t2,t1,0,7,16); - WHBFY(t2,t1,0,8,16); WHBFY(t2,t1,0,9,16); WHBFY(t2,t1,0,10,16); WHBFY(t2,t1,0,11,16); - WHBFY(t2,t1,0,12,16); WHBFY(t2,t1,0,13,16); WHBFY(t2,t1,0,14,16); WHBFY(t2,t1,0,15,16); - - WHBFY(t2,t1,32,0,16); WHBFY(t2,t1,32,1,16); WHBFY(t2,t1,32,2,16); WHBFY(t2,t1,32,3,16); - WHBFY(t2,t1,32,4,16); WHBFY(t2,t1,32,5,16); WHBFY(t2,t1,32,6,16); WHBFY(t2,t1,32,7,16); - WHBFY(t2,t1,32,8,16); WHBFY(t2,t1,32,9,16); WHBFY(t2,t1,32,10,16); WHBFY(t2,t1,32,11,16); - WHBFY(t2,t1,32,12,16); WHBFY(t2,t1,32,13,16); WHBFY(t2,t1,32,14,16); WHBFY(t2,t1,32,15,16); - - // group 3 - WHBFY(t1,t2,0,0,8); WHBFY(t1,t2,0,1,8); WHBFY(t1,t2,0,2,8); WHBFY(t1,t2,0,3,8); - WHBFY(t1,t2,0,4,8); WHBFY(t1,t2,0,5,8); WHBFY(t1,t2,0,6,8); WHBFY(t1,t2,0,7,8); - WHBFY(t1,t2,16,0,8); WHBFY(t1,t2,16,1,8); WHBFY(t1,t2,16,2,8); WHBFY(t1,t2,16,3,8); - WHBFY(t1,t2,16,4,8); WHBFY(t1,t2,16,5,8); WHBFY(t1,t2,16,6,8); WHBFY(t1,t2,16,7,8); - WHBFY(t1,t2,32,0,8); WHBFY(t1,t2,32,1,8); WHBFY(t1,t2,32,2,8); WHBFY(t1,t2,32,3,8); - WHBFY(t1,t2,32,4,8); WHBFY(t1,t2,32,5,8); WHBFY(t1,t2,32,6,8); WHBFY(t1,t2,32,7,8); - WHBFY(t1,t2,48,0,8); WHBFY(t1,t2,48,1,8); WHBFY(t1,t2,48,2,8); WHBFY(t1,t2,48,3,8); - WHBFY(t1,t2,48,4,8); WHBFY(t1,t2,48,5,8); WHBFY(t1,t2,48,6,8); WHBFY(t1,t2,48,7,8); - - // group 4 - WHBFY(t2,t1,0,0,4); WHBFY(t2,t1,0,1,4); WHBFY(t2,t1,0,2,4); WHBFY(t2,t1,0,3,4); - WHBFY(t2,t1,8,0,4); WHBFY(t2,t1,8,1,4); WHBFY(t2,t1,8,2,4); WHBFY(t2,t1,8,3,4); - WHBFY(t2,t1,16,0,4); WHBFY(t2,t1,16,1,4); WHBFY(t2,t1,16,2,4); WHBFY(t2,t1,16,3,4); - WHBFY(t2,t1,24,0,4); WHBFY(t2,t1,24,1,4); WHBFY(t2,t1,24,2,4); WHBFY(t2,t1,24,3,4); - WHBFY(t2,t1,32,0,4); WHBFY(t2,t1,32,1,4); WHBFY(t2,t1,32,2,4); WHBFY(t2,t1,32,3,4); - WHBFY(t2,t1,40,0,4); WHBFY(t2,t1,40,1,4); WHBFY(t2,t1,40,2,4); WHBFY(t2,t1,40,3,4); - WHBFY(t2,t1,48,0,4); WHBFY(t2,t1,48,1,4); WHBFY(t2,t1,48,2,4); WHBFY(t2,t1,48,3,4); - WHBFY(t2,t1,56,0,4); WHBFY(t2,t1,56,1,4); WHBFY(t2,t1,56,2,4); WHBFY(t2,t1,56,3,4); - - // group 5 - WHBFY(t1,t2,0,0,2); WHBFY(t1,t2,0,1,2); WHBFY(t1,t2,4,0,2); WHBFY(t1,t2,4,1,2); - WHBFY(t1,t2,8,0,2); WHBFY(t1,t2,8,1,2); WHBFY(t1,t2,12,0,2); WHBFY(t1,t2,12,1,2); - WHBFY(t1,t2,16,0,2); WHBFY(t1,t2,16,1,2); WHBFY(t1,t2,20,0,2); WHBFY(t1,t2,20,1,2); - WHBFY(t1,t2,24,0,2); WHBFY(t1,t2,24,1,2); WHBFY(t1,t2,28,0,2); WHBFY(t1,t2,28,1,2); - WHBFY(t1,t2,32,0,2); WHBFY(t1,t2,32,1,2); WHBFY(t1,t2,36,0,2); WHBFY(t1,t2,36,1,2); - WHBFY(t1,t2,40,0,2); WHBFY(t1,t2,40,1,2); WHBFY(t1,t2,44,0,2); WHBFY(t1,t2,44,1,2); - WHBFY(t1,t2,48,0,2); WHBFY(t1,t2,48,1,2); WHBFY(t1,t2,52,0,2); WHBFY(t1,t2,52,1,2); - WHBFY(t1,t2,56,0,2); WHBFY(t1,t2,56,1,2); WHBFY(t1,t2,60,0,2); WHBFY(t1,t2,60,1,2); - - // group 6 - WHBFY(dst,t1,0,0,1); WHBFY(dst,t1,2,0,1); WHBFY(dst,t1,4,0,1); WHBFY(dst,t1,6,0,1); - WHBFY(dst,t1,8,0,1); WHBFY(dst,t1,10,0,1); WHBFY(dst,t1,12,0,1); WHBFY(dst,t1,14,0,1); - WHBFY(dst,t1,16,0,1); WHBFY(dst,t1,18,0,1); WHBFY(dst,t1,20,0,1); WHBFY(dst,t1,22,0,1); - WHBFY(dst,t1,24,0,1); WHBFY(dst,t1,26,0,1); WHBFY(dst,t1,28,0,1); WHBFY(dst,t1,30,0,1); - WHBFY(dst,t1,32,0,1); WHBFY(dst,t1,34,0,1); WHBFY(dst,t1,36,0,1); WHBFY(dst,t1,38,0,1); - WHBFY(dst,t1,40,0,1); WHBFY(dst,t1,42,0,1); WHBFY(dst,t1,44,0,1); WHBFY(dst,t1,46,0,1); - WHBFY(dst,t1,48,0,1); WHBFY(dst,t1,50,0,1); WHBFY(dst,t1,52,0,1); WHBFY(dst,t1,54,0,1); - WHBFY(dst,t1,56,0,1); WHBFY(dst,t1,58,0,1); WHBFY(dst,t1,60,0,1); WHBFY(dst,t1,62,0,1); +// npfwht.c +// Basic implementation of the Fast Walsh-Hadamard Transforms +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#include "npfwht.h" + +#define WHBFY(dst,src,base,offs,dist) { dst[base+offs]=src[base+offs]+src[base+offs+dist]; dst[base+offs+dist]=src[base+offs]-src[base+offs+dist]; } + +typedef void (*pnp_fwht)(float*,float*); + +static void np_fwht2(float *dst, float *src); + +static void np_fwht1(float *dst, float *src); +static void np_fwht2(float *dst, float *src); +static void np_fwht4(float *dst, float *src); +static void np_fwht8(float *dst, float *src); +static void np_fwht16(float *dst, float *src); +static void np_fwht32(float *dst, float *src); +static void np_fwht64(float *dst, float *src); + +static pnp_fwht np_fwht_tab[7] = { + np_fwht1, + np_fwht2, + np_fwht4, + np_fwht8, + np_fwht16, + np_fwht32, + np_fwht64 +}; + +void np_fwht(int nlogdim, float *dst, float *src) +{ + np_fwht_tab[nlogdim](dst,src); +} + +static void np_fwht1(float *dst, float *src) +{ + dst[0] = src[0]; +} + + +static void np_fwht2(float *dst, float *src) +{ + float t[2]; + + WHBFY(t,src,0,0,1); + dst[0]= t[0]; + dst[1]= t[1]; +} + +static void np_fwht4(float *dst, float *src) +{ + float t[4]; + + // group 1 + WHBFY(t,src,0,0,2); WHBFY(t,src,0,1,2); + // group 2 + WHBFY(dst,t,0,0,1); WHBFY(dst,t,2,0,1); +}; + + +static void np_fwht8(float *dst, float *src) +{ + float t[16]; + float *t1=t, *t2=t+8; + + // group 1 + WHBFY(t1,src,0,0,4); WHBFY(t1,src,0,1,4); WHBFY(t1,src,0,2,4); WHBFY(t1,src,0,3,4); + // group 2 + WHBFY(t2,t1,0,0,2); WHBFY(t2,t1,0,1,2); WHBFY(t2,t1,4,0,2); WHBFY(t2,t1,4,1,2); + // group 3 + WHBFY(dst,t2,0,0,1); WHBFY(dst,t2,2,0,1); WHBFY(dst,t2,4,0,1); WHBFY(dst,t2,6,0,1); +}; + + +static void np_fwht16(float *dst, float *src) +{ + float t[32]; + float *t1=t, *t2=t+16; + + // group 1 + WHBFY(t1,src,0,0,8); WHBFY(t1,src,0,1,8); WHBFY(t1,src,0,2,8); WHBFY(t1,src,0,3,8); + WHBFY(t1,src,0,4,8); WHBFY(t1,src,0,5,8); WHBFY(t1,src,0,6,8); WHBFY(t1,src,0,7,8); + // group 2 + WHBFY(t2,t1,0,0,4); WHBFY(t2,t1,0,1,4); WHBFY(t2,t1,0,2,4); WHBFY(t2,t1,0,3,4); + WHBFY(t2,t1,8,0,4); WHBFY(t2,t1,8,1,4); WHBFY(t2,t1,8,2,4); WHBFY(t2,t1,8,3,4); + // group 3 + WHBFY(t1,t2,0,0,2); WHBFY(t1,t2,0,1,2); WHBFY(t1,t2,4,0,2); WHBFY(t1,t2,4,1,2); + WHBFY(t1,t2,8,0,2); WHBFY(t1,t2,8,1,2); WHBFY(t1,t2,12,0,2); WHBFY(t1,t2,12,1,2); + // group 4 + WHBFY(dst,t1,0,0,1); WHBFY(dst,t1,2,0,1); WHBFY(dst,t1,4,0,1); WHBFY(dst,t1,6,0,1); + WHBFY(dst,t1,8,0,1); WHBFY(dst,t1,10,0,1); WHBFY(dst,t1,12,0,1); WHBFY(dst,t1,14,0,1); + +} + +static void np_fwht32(float *dst, float *src) +{ + float t[64]; + float *t1=t, *t2=t+32; + + // group 1 + WHBFY(t1,src,0,0,16); WHBFY(t1,src,0,1,16); WHBFY(t1,src,0,2,16); WHBFY(t1,src,0,3,16); + WHBFY(t1,src,0,4,16); WHBFY(t1,src,0,5,16); WHBFY(t1,src,0,6,16); WHBFY(t1,src,0,7,16); + WHBFY(t1,src,0,8,16); WHBFY(t1,src,0,9,16); WHBFY(t1,src,0,10,16); WHBFY(t1,src,0,11,16); + WHBFY(t1,src,0,12,16); WHBFY(t1,src,0,13,16); WHBFY(t1,src,0,14,16); WHBFY(t1,src,0,15,16); + + // group 2 + WHBFY(t2,t1,0,0,8); WHBFY(t2,t1,0,1,8); WHBFY(t2,t1,0,2,8); WHBFY(t2,t1,0,3,8); + WHBFY(t2,t1,0,4,8); WHBFY(t2,t1,0,5,8); WHBFY(t2,t1,0,6,8); WHBFY(t2,t1,0,7,8); + WHBFY(t2,t1,16,0,8); WHBFY(t2,t1,16,1,8); WHBFY(t2,t1,16,2,8); WHBFY(t2,t1,16,3,8); + WHBFY(t2,t1,16,4,8); WHBFY(t2,t1,16,5,8); WHBFY(t2,t1,16,6,8); WHBFY(t2,t1,16,7,8); + + // group 3 + WHBFY(t1,t2,0,0,4); WHBFY(t1,t2,0,1,4); WHBFY(t1,t2,0,2,4); WHBFY(t1,t2,0,3,4); + WHBFY(t1,t2,8,0,4); WHBFY(t1,t2,8,1,4); WHBFY(t1,t2,8,2,4); WHBFY(t1,t2,8,3,4); + WHBFY(t1,t2,16,0,4); WHBFY(t1,t2,16,1,4); WHBFY(t1,t2,16,2,4); WHBFY(t1,t2,16,3,4); + WHBFY(t1,t2,24,0,4); WHBFY(t1,t2,24,1,4); WHBFY(t1,t2,24,2,4); WHBFY(t1,t2,24,3,4); + + // group 4 + WHBFY(t2,t1,0,0,2); WHBFY(t2,t1,0,1,2); WHBFY(t2,t1,4,0,2); WHBFY(t2,t1,4,1,2); + WHBFY(t2,t1,8,0,2); WHBFY(t2,t1,8,1,2); WHBFY(t2,t1,12,0,2); WHBFY(t2,t1,12,1,2); + WHBFY(t2,t1,16,0,2); WHBFY(t2,t1,16,1,2); WHBFY(t2,t1,20,0,2); WHBFY(t2,t1,20,1,2); + WHBFY(t2,t1,24,0,2); WHBFY(t2,t1,24,1,2); WHBFY(t2,t1,28,0,2); WHBFY(t2,t1,28,1,2); + + // group 5 + WHBFY(dst,t2,0,0,1); WHBFY(dst,t2,2,0,1); WHBFY(dst,t2,4,0,1); WHBFY(dst,t2,6,0,1); + WHBFY(dst,t2,8,0,1); WHBFY(dst,t2,10,0,1); WHBFY(dst,t2,12,0,1); WHBFY(dst,t2,14,0,1); + WHBFY(dst,t2,16,0,1); WHBFY(dst,t2,18,0,1); WHBFY(dst,t2,20,0,1); WHBFY(dst,t2,22,0,1); + WHBFY(dst,t2,24,0,1); WHBFY(dst,t2,26,0,1); WHBFY(dst,t2,28,0,1); WHBFY(dst,t2,30,0,1); + +} + +static void np_fwht64(float *dst, float *src) +{ + float t[128]; + float *t1=t, *t2=t+64; + + + // group 1 + WHBFY(t1,src,0,0,32); WHBFY(t1,src,0,1,32); WHBFY(t1,src,0,2,32); WHBFY(t1,src,0,3,32); + WHBFY(t1,src,0,4,32); WHBFY(t1,src,0,5,32); WHBFY(t1,src,0,6,32); WHBFY(t1,src,0,7,32); + WHBFY(t1,src,0,8,32); WHBFY(t1,src,0,9,32); WHBFY(t1,src,0,10,32); WHBFY(t1,src,0,11,32); + WHBFY(t1,src,0,12,32); WHBFY(t1,src,0,13,32); WHBFY(t1,src,0,14,32); WHBFY(t1,src,0,15,32); + WHBFY(t1,src,0,16,32); WHBFY(t1,src,0,17,32); WHBFY(t1,src,0,18,32); WHBFY(t1,src,0,19,32); + WHBFY(t1,src,0,20,32); WHBFY(t1,src,0,21,32); WHBFY(t1,src,0,22,32); WHBFY(t1,src,0,23,32); + WHBFY(t1,src,0,24,32); WHBFY(t1,src,0,25,32); WHBFY(t1,src,0,26,32); WHBFY(t1,src,0,27,32); + WHBFY(t1,src,0,28,32); WHBFY(t1,src,0,29,32); WHBFY(t1,src,0,30,32); WHBFY(t1,src,0,31,32); + + // group 2 + WHBFY(t2,t1,0,0,16); WHBFY(t2,t1,0,1,16); WHBFY(t2,t1,0,2,16); WHBFY(t2,t1,0,3,16); + WHBFY(t2,t1,0,4,16); WHBFY(t2,t1,0,5,16); WHBFY(t2,t1,0,6,16); WHBFY(t2,t1,0,7,16); + WHBFY(t2,t1,0,8,16); WHBFY(t2,t1,0,9,16); WHBFY(t2,t1,0,10,16); WHBFY(t2,t1,0,11,16); + WHBFY(t2,t1,0,12,16); WHBFY(t2,t1,0,13,16); WHBFY(t2,t1,0,14,16); WHBFY(t2,t1,0,15,16); + + WHBFY(t2,t1,32,0,16); WHBFY(t2,t1,32,1,16); WHBFY(t2,t1,32,2,16); WHBFY(t2,t1,32,3,16); + WHBFY(t2,t1,32,4,16); WHBFY(t2,t1,32,5,16); WHBFY(t2,t1,32,6,16); WHBFY(t2,t1,32,7,16); + WHBFY(t2,t1,32,8,16); WHBFY(t2,t1,32,9,16); WHBFY(t2,t1,32,10,16); WHBFY(t2,t1,32,11,16); + WHBFY(t2,t1,32,12,16); WHBFY(t2,t1,32,13,16); WHBFY(t2,t1,32,14,16); WHBFY(t2,t1,32,15,16); + + // group 3 + WHBFY(t1,t2,0,0,8); WHBFY(t1,t2,0,1,8); WHBFY(t1,t2,0,2,8); WHBFY(t1,t2,0,3,8); + WHBFY(t1,t2,0,4,8); WHBFY(t1,t2,0,5,8); WHBFY(t1,t2,0,6,8); WHBFY(t1,t2,0,7,8); + WHBFY(t1,t2,16,0,8); WHBFY(t1,t2,16,1,8); WHBFY(t1,t2,16,2,8); WHBFY(t1,t2,16,3,8); + WHBFY(t1,t2,16,4,8); WHBFY(t1,t2,16,5,8); WHBFY(t1,t2,16,6,8); WHBFY(t1,t2,16,7,8); + WHBFY(t1,t2,32,0,8); WHBFY(t1,t2,32,1,8); WHBFY(t1,t2,32,2,8); WHBFY(t1,t2,32,3,8); + WHBFY(t1,t2,32,4,8); WHBFY(t1,t2,32,5,8); WHBFY(t1,t2,32,6,8); WHBFY(t1,t2,32,7,8); + WHBFY(t1,t2,48,0,8); WHBFY(t1,t2,48,1,8); WHBFY(t1,t2,48,2,8); WHBFY(t1,t2,48,3,8); + WHBFY(t1,t2,48,4,8); WHBFY(t1,t2,48,5,8); WHBFY(t1,t2,48,6,8); WHBFY(t1,t2,48,7,8); + + // group 4 + WHBFY(t2,t1,0,0,4); WHBFY(t2,t1,0,1,4); WHBFY(t2,t1,0,2,4); WHBFY(t2,t1,0,3,4); + WHBFY(t2,t1,8,0,4); WHBFY(t2,t1,8,1,4); WHBFY(t2,t1,8,2,4); WHBFY(t2,t1,8,3,4); + WHBFY(t2,t1,16,0,4); WHBFY(t2,t1,16,1,4); WHBFY(t2,t1,16,2,4); WHBFY(t2,t1,16,3,4); + WHBFY(t2,t1,24,0,4); WHBFY(t2,t1,24,1,4); WHBFY(t2,t1,24,2,4); WHBFY(t2,t1,24,3,4); + WHBFY(t2,t1,32,0,4); WHBFY(t2,t1,32,1,4); WHBFY(t2,t1,32,2,4); WHBFY(t2,t1,32,3,4); + WHBFY(t2,t1,40,0,4); WHBFY(t2,t1,40,1,4); WHBFY(t2,t1,40,2,4); WHBFY(t2,t1,40,3,4); + WHBFY(t2,t1,48,0,4); WHBFY(t2,t1,48,1,4); WHBFY(t2,t1,48,2,4); WHBFY(t2,t1,48,3,4); + WHBFY(t2,t1,56,0,4); WHBFY(t2,t1,56,1,4); WHBFY(t2,t1,56,2,4); WHBFY(t2,t1,56,3,4); + + // group 5 + WHBFY(t1,t2,0,0,2); WHBFY(t1,t2,0,1,2); WHBFY(t1,t2,4,0,2); WHBFY(t1,t2,4,1,2); + WHBFY(t1,t2,8,0,2); WHBFY(t1,t2,8,1,2); WHBFY(t1,t2,12,0,2); WHBFY(t1,t2,12,1,2); + WHBFY(t1,t2,16,0,2); WHBFY(t1,t2,16,1,2); WHBFY(t1,t2,20,0,2); WHBFY(t1,t2,20,1,2); + WHBFY(t1,t2,24,0,2); WHBFY(t1,t2,24,1,2); WHBFY(t1,t2,28,0,2); WHBFY(t1,t2,28,1,2); + WHBFY(t1,t2,32,0,2); WHBFY(t1,t2,32,1,2); WHBFY(t1,t2,36,0,2); WHBFY(t1,t2,36,1,2); + WHBFY(t1,t2,40,0,2); WHBFY(t1,t2,40,1,2); WHBFY(t1,t2,44,0,2); WHBFY(t1,t2,44,1,2); + WHBFY(t1,t2,48,0,2); WHBFY(t1,t2,48,1,2); WHBFY(t1,t2,52,0,2); WHBFY(t1,t2,52,1,2); + WHBFY(t1,t2,56,0,2); WHBFY(t1,t2,56,1,2); WHBFY(t1,t2,60,0,2); WHBFY(t1,t2,60,1,2); + + // group 6 + WHBFY(dst,t1,0,0,1); WHBFY(dst,t1,2,0,1); WHBFY(dst,t1,4,0,1); WHBFY(dst,t1,6,0,1); + WHBFY(dst,t1,8,0,1); WHBFY(dst,t1,10,0,1); WHBFY(dst,t1,12,0,1); WHBFY(dst,t1,14,0,1); + WHBFY(dst,t1,16,0,1); WHBFY(dst,t1,18,0,1); WHBFY(dst,t1,20,0,1); WHBFY(dst,t1,22,0,1); + WHBFY(dst,t1,24,0,1); WHBFY(dst,t1,26,0,1); WHBFY(dst,t1,28,0,1); WHBFY(dst,t1,30,0,1); + WHBFY(dst,t1,32,0,1); WHBFY(dst,t1,34,0,1); WHBFY(dst,t1,36,0,1); WHBFY(dst,t1,38,0,1); + WHBFY(dst,t1,40,0,1); WHBFY(dst,t1,42,0,1); WHBFY(dst,t1,44,0,1); WHBFY(dst,t1,46,0,1); + WHBFY(dst,t1,48,0,1); WHBFY(dst,t1,50,0,1); WHBFY(dst,t1,52,0,1); WHBFY(dst,t1,54,0,1); + WHBFY(dst,t1,56,0,1); WHBFY(dst,t1,58,0,1); WHBFY(dst,t1,60,0,1); WHBFY(dst,t1,62,0,1); } \ No newline at end of file diff --git a/qracodes-mt/npfwht.h b/qracodes-st/npfwht.h similarity index 94% rename from qracodes-mt/npfwht.h rename to qracodes-st/npfwht.h index 212b0af..e31ebd5 100644 --- a/qracodes-mt/npfwht.h +++ b/qracodes-st/npfwht.h @@ -1,45 +1,45 @@ -// np_fwht.h -// Basic implementation of the Fast Walsh-Hadamard Transforms -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -#ifndef _npfwht_h_ -#define _npfwht_h_ - -#ifdef __cplusplus -extern "C" { -#endif - -void np_fwht(int nlogdim, float *dst, float *src); -// Compute the Walsh-Hadamard transform of the given data up to a -// 64-dimensional transform -// -// Input parameters: -// nlogdim: log2 of the transform size. Must be in the range [0..6] -// src : pointer to the input data buffer. -// dst : pointer to the output data buffer. -// -// src and dst must point to preallocated data buffers of size 2^nlogdim*sizeof(float) -// src and dst buffers can overlap - -#ifdef __cplusplus -} -#endif - -#endif // _npfwht_ +// np_fwht.h +// Basic implementation of the Fast Walsh-Hadamard Transforms +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#ifndef _npfwht_h_ +#define _npfwht_h_ + +#ifdef __cplusplus +extern "C" { +#endif + +void np_fwht(int nlogdim, float *dst, float *src); +// Compute the Walsh-Hadamard transform of the given data up to a +// 64-dimensional transform +// +// Input parameters: +// nlogdim: log2 of the transform size. Must be in the range [0..6] +// src : pointer to the input data buffer. +// dst : pointer to the output data buffer. +// +// src and dst must point to preallocated data buffers of size 2^nlogdim*sizeof(float) +// src and dst buffers can overlap + +#ifdef __cplusplus +} +#endif + +#endif // _npfwht_ diff --git a/qracodes-mt/pdmath.c b/qracodes-st/pdmath.c similarity index 94% rename from qracodes-mt/pdmath.c rename to qracodes-st/pdmath.c index 2c55e41..ba15c2e 100644 --- a/qracodes-mt/pdmath.c +++ b/qracodes-st/pdmath.c @@ -1,385 +1,385 @@ -// pdmath.c -// Elementary math on probability distributions -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - -#include "pdmath.h" - -typedef const float *ppd_uniform; -typedef void (*ppd_imul)(float*,const float*); -typedef float (*ppd_norm)(float*); - -// define vector size in function of its logarithm in base 2 -static const int pd_log2dim[7] = { - 1,2,4,8,16,32,64 -}; - -// define uniform distributions of given size -static const float pd_uniform1[1] = { - 1. -}; -static const float pd_uniform2[2] = { - 1./2., 1./2. -}; -static const float pd_uniform4[4] = { - 1./4., 1./4.,1./4., 1./4. -}; -static const float pd_uniform8[8] = { - 1./8., 1./8.,1./8., 1./8.,1./8., 1./8.,1./8., 1./8. -}; -static const float pd_uniform16[16] = { - 1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16., - 1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16. -}; -static const float pd_uniform32[32] = { - 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., - 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., - 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., - 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32. -}; -static const float pd_uniform64[64] = { - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., - 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64. - -}; - -static const ppd_uniform pd_uniform_tab[7] = { - pd_uniform1, - pd_uniform2, - pd_uniform4, - pd_uniform8, - pd_uniform16, - pd_uniform32, - pd_uniform64 -}; - -// returns a pointer to the uniform distribution of the given logsize -const float *pd_uniform(int nlogdim) -{ - return pd_uniform_tab[nlogdim]; -} - -// in-place multiplication functions -// compute dst = dst*src for any element of the distrib - -static void pd_imul1(float *dst, const float *src) -{ - dst[0] *= src[0]; -} - -static void pd_imul2(float *dst, const float *src) -{ - dst[0] *= src[0]; dst[1] *= src[1]; -} -static void pd_imul4(float *dst, const float *src) -{ - dst[0] *= src[0]; dst[1] *= src[1]; - dst[2] *= src[2]; dst[3] *= src[3]; -} -static void pd_imul8(float *dst, const float *src) -{ - dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3]; - dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7]; -} -static void pd_imul16(float *dst, const float *src) -{ - dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3]; - dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7]; - dst[8] *= src[8]; dst[9] *= src[9]; dst[10]*= src[10]; dst[11]*= src[11]; - dst[12]*= src[12]; dst[13]*= src[13]; dst[14]*= src[14]; dst[15]*= src[15]; -} -static void pd_imul32(float *dst, const float *src) -{ - pd_imul16(dst,src); - pd_imul16(dst+16,src+16); -} -static void pd_imul64(float *dst, const float *src) -{ - pd_imul16(dst, src); - pd_imul16(dst+16, src+16); - pd_imul16(dst+32, src+32); - pd_imul16(dst+48, src+48); -} - -static const ppd_imul pd_imul_tab[7] = { - pd_imul1, - pd_imul2, - pd_imul4, - pd_imul8, - pd_imul16, - pd_imul32, - pd_imul64 -}; - -// in place multiplication -// compute dst = dst*src for any element of the distrib give their log2 size -// arguments must be pointers to array of floats of the given size -void pd_imul(float *dst, const float *src, int nlogdim) -{ - pd_imul_tab[nlogdim](dst,src); -} - -static float pd_norm1(float *ppd) -{ - float t = ppd[0]; - ppd[0] = 1.f; - return t; -} - -static float pd_norm2(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; - - if (t<=0) { - pd_init(ppd,pd_uniform(1),pd_log2dim[1]); - return t; - } - - to = t; - t = 1.f/t; - ppd[0] *=t; ppd[1] *=t; - return to; - -} - -static float pd_norm4(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; - - if (t<=0) { - pd_init(ppd,pd_uniform(2),pd_log2dim[2]); - return t; - } - - to = t; - t = 1.f/t; - ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; - return to; -} - -static float pd_norm8(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; - t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; - - if (t<=0) { - pd_init(ppd,pd_uniform(3),pd_log2dim[3]); - return t; - } - - to = t; - t = 1.f/t; - ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; - ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; - return to; -} -static float pd_norm16(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; - t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; - t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; - t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; - - if (t<=0) { - pd_init(ppd,pd_uniform(4),pd_log2dim[4]); - return t; - } - - to = t; - t = 1.f/t; - ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; - ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; - ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; - ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; - - return to; -} -static float pd_norm32(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; - t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; - t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; - t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; - t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19]; - t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23]; - t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27]; - t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31]; - - if (t<=0) { - pd_init(ppd,pd_uniform(5),pd_log2dim[5]); - return t; - } - - to = t; - t = 1.f/t; - ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; - ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; - ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; - ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; - ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t; - ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t; - ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t; - ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t; - - return to; -} - -static float pd_norm64(float *ppd) -{ - float t,to; - - t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; - t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; - t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; - t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; - t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19]; - t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23]; - t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27]; - t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31]; - - t +=ppd[32]; t +=ppd[33]; t +=ppd[34]; t +=ppd[35]; - t +=ppd[36]; t +=ppd[37]; t +=ppd[38]; t +=ppd[39]; - t +=ppd[40]; t +=ppd[41]; t +=ppd[42]; t +=ppd[43]; - t +=ppd[44]; t +=ppd[45]; t +=ppd[46]; t +=ppd[47]; - t +=ppd[48]; t +=ppd[49]; t +=ppd[50]; t +=ppd[51]; - t +=ppd[52]; t +=ppd[53]; t +=ppd[54]; t +=ppd[55]; - t +=ppd[56]; t +=ppd[57]; t +=ppd[58]; t +=ppd[59]; - t +=ppd[60]; t +=ppd[61]; t +=ppd[62]; t +=ppd[63]; - - if (t<=0) { - pd_init(ppd,pd_uniform(6),pd_log2dim[6]); - return t; - } - - to = t; - t = 1.0f/t; - ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; - ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; - ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; - ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; - ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t; - ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t; - ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t; - ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t; - - ppd[32] *=t; ppd[33] *=t; ppd[34] *=t; ppd[35] *=t; - ppd[36] *=t; ppd[37] *=t; ppd[38] *=t; ppd[39] *=t; - ppd[40] *=t; ppd[41] *=t; ppd[42] *=t; ppd[43] *=t; - ppd[44] *=t; ppd[45] *=t; ppd[46] *=t; ppd[47] *=t; - ppd[48] *=t; ppd[49] *=t; ppd[50] *=t; ppd[51] *=t; - ppd[52] *=t; ppd[53] *=t; ppd[54] *=t; ppd[55] *=t; - ppd[56] *=t; ppd[57] *=t; ppd[58] *=t; ppd[59] *=t; - ppd[60] *=t; ppd[61] *=t; ppd[62] *=t; ppd[63] *=t; - - return to; -} - - -static const ppd_norm pd_norm_tab[7] = { - pd_norm1, - pd_norm2, - pd_norm4, - pd_norm8, - pd_norm16, - pd_norm32, - pd_norm64 -}; - -float pd_norm(float *pd, int nlogdim) -{ - return pd_norm_tab[nlogdim](pd); -} - -void pd_memset(float *dst, const float *src, int ndim, int nitems) -{ - int size = PD_SIZE(ndim); - while(nitems--) { - memcpy(dst,src,size); - dst +=ndim; - } -} - -void pd_fwdperm(float *dst, float *src, const int *perm, int ndim) -{ - // TODO: non-loop implementation - while (ndim--) - dst[ndim] = src[perm[ndim]]; -} - -void pd_bwdperm(float *dst, float *src, const int *perm, int ndim) -{ - // TODO: non-loop implementation - while (ndim--) - dst[perm[ndim]] = src[ndim]; -} - -float pd_max(float *src, int ndim) -{ - // TODO: faster implementation - - float cmax=0; // we assume that prob distributions are always positive - float cval; - - while (ndim--) { - cval = src[ndim]; - if (cval>=cmax) { - cmax = cval; - } - } - - return cmax; -} - -int pd_argmax(float *pmax, float *src, int ndim) -{ - // TODO: faster implementation - - float cmax=0; // we assume that prob distributions are always positive - float cval; - int idxmax=-1; // indicates that all pd elements are <0 - - while (ndim--) { - cval = src[ndim]; - if (cval>=cmax) { - cmax = cval; - idxmax = ndim; - } - } - - if (pmax) - *pmax = cmax; - - return idxmax; -} +// pdmath.c +// Elementary math on probability distributions +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#include "pdmath.h" + +typedef const float *ppd_uniform; +typedef void (*ppd_imul)(float*,const float*); +typedef float (*ppd_norm)(float*); + +// define vector size in function of its logarithm in base 2 +static const int pd_log2dim[7] = { + 1,2,4,8,16,32,64 +}; + +// define uniform distributions of given size +static const float pd_uniform1[1] = { + 1. +}; +static const float pd_uniform2[2] = { + 1./2., 1./2. +}; +static const float pd_uniform4[4] = { + 1./4., 1./4.,1./4., 1./4. +}; +static const float pd_uniform8[8] = { + 1./8., 1./8.,1./8., 1./8.,1./8., 1./8.,1./8., 1./8. +}; +static const float pd_uniform16[16] = { + 1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16., + 1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16. +}; +static const float pd_uniform32[32] = { + 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., + 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., + 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32., + 1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32. +}; +static const float pd_uniform64[64] = { + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64., + 1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64. + +}; + +static const ppd_uniform pd_uniform_tab[7] = { + pd_uniform1, + pd_uniform2, + pd_uniform4, + pd_uniform8, + pd_uniform16, + pd_uniform32, + pd_uniform64 +}; + +// returns a pointer to the uniform distribution of the given logsize +const float *pd_uniform(int nlogdim) +{ + return pd_uniform_tab[nlogdim]; +} + +// in-place multiplication functions +// compute dst = dst*src for any element of the distrib + +static void pd_imul1(float *dst, const float *src) +{ + dst[0] *= src[0]; +} + +static void pd_imul2(float *dst, const float *src) +{ + dst[0] *= src[0]; dst[1] *= src[1]; +} +static void pd_imul4(float *dst, const float *src) +{ + dst[0] *= src[0]; dst[1] *= src[1]; + dst[2] *= src[2]; dst[3] *= src[3]; +} +static void pd_imul8(float *dst, const float *src) +{ + dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3]; + dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7]; +} +static void pd_imul16(float *dst, const float *src) +{ + dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3]; + dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7]; + dst[8] *= src[8]; dst[9] *= src[9]; dst[10]*= src[10]; dst[11]*= src[11]; + dst[12]*= src[12]; dst[13]*= src[13]; dst[14]*= src[14]; dst[15]*= src[15]; +} +static void pd_imul32(float *dst, const float *src) +{ + pd_imul16(dst,src); + pd_imul16(dst+16,src+16); +} +static void pd_imul64(float *dst, const float *src) +{ + pd_imul16(dst, src); + pd_imul16(dst+16, src+16); + pd_imul16(dst+32, src+32); + pd_imul16(dst+48, src+48); +} + +static const ppd_imul pd_imul_tab[7] = { + pd_imul1, + pd_imul2, + pd_imul4, + pd_imul8, + pd_imul16, + pd_imul32, + pd_imul64 +}; + +// in place multiplication +// compute dst = dst*src for any element of the distrib give their log2 size +// arguments must be pointers to array of floats of the given size +void pd_imul(float *dst, const float *src, int nlogdim) +{ + pd_imul_tab[nlogdim](dst,src); +} + +static float pd_norm1(float *ppd) +{ + float t = ppd[0]; + ppd[0] = 1.f; + return t; +} + +static float pd_norm2(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; + + if (t<=0) { + pd_init(ppd,pd_uniform(1),pd_log2dim[1]); + return t; + } + + to = t; + t = 1.f/t; + ppd[0] *=t; ppd[1] *=t; + return to; + +} + +static float pd_norm4(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; + + if (t<=0) { + pd_init(ppd,pd_uniform(2),pd_log2dim[2]); + return t; + } + + to = t; + t = 1.f/t; + ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; + return to; +} + +static float pd_norm8(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; + t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; + + if (t<=0) { + pd_init(ppd,pd_uniform(3),pd_log2dim[3]); + return t; + } + + to = t; + t = 1.f/t; + ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; + ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; + return to; +} +static float pd_norm16(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; + t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; + t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; + t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; + + if (t<=0) { + pd_init(ppd,pd_uniform(4),pd_log2dim[4]); + return t; + } + + to = t; + t = 1.f/t; + ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; + ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; + ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; + ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; + + return to; +} +static float pd_norm32(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; + t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; + t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; + t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; + t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19]; + t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23]; + t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27]; + t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31]; + + if (t<=0) { + pd_init(ppd,pd_uniform(5),pd_log2dim[5]); + return t; + } + + to = t; + t = 1.f/t; + ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; + ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; + ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; + ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; + ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t; + ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t; + ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t; + ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t; + + return to; +} + +static float pd_norm64(float *ppd) +{ + float t,to; + + t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3]; + t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7]; + t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11]; + t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15]; + t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19]; + t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23]; + t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27]; + t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31]; + + t +=ppd[32]; t +=ppd[33]; t +=ppd[34]; t +=ppd[35]; + t +=ppd[36]; t +=ppd[37]; t +=ppd[38]; t +=ppd[39]; + t +=ppd[40]; t +=ppd[41]; t +=ppd[42]; t +=ppd[43]; + t +=ppd[44]; t +=ppd[45]; t +=ppd[46]; t +=ppd[47]; + t +=ppd[48]; t +=ppd[49]; t +=ppd[50]; t +=ppd[51]; + t +=ppd[52]; t +=ppd[53]; t +=ppd[54]; t +=ppd[55]; + t +=ppd[56]; t +=ppd[57]; t +=ppd[58]; t +=ppd[59]; + t +=ppd[60]; t +=ppd[61]; t +=ppd[62]; t +=ppd[63]; + + if (t<=0) { + pd_init(ppd,pd_uniform(6),pd_log2dim[6]); + return t; + } + + to = t; + t = 1.0f/t; + ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t; + ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t; + ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t; + ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t; + ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t; + ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t; + ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t; + ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t; + + ppd[32] *=t; ppd[33] *=t; ppd[34] *=t; ppd[35] *=t; + ppd[36] *=t; ppd[37] *=t; ppd[38] *=t; ppd[39] *=t; + ppd[40] *=t; ppd[41] *=t; ppd[42] *=t; ppd[43] *=t; + ppd[44] *=t; ppd[45] *=t; ppd[46] *=t; ppd[47] *=t; + ppd[48] *=t; ppd[49] *=t; ppd[50] *=t; ppd[51] *=t; + ppd[52] *=t; ppd[53] *=t; ppd[54] *=t; ppd[55] *=t; + ppd[56] *=t; ppd[57] *=t; ppd[58] *=t; ppd[59] *=t; + ppd[60] *=t; ppd[61] *=t; ppd[62] *=t; ppd[63] *=t; + + return to; +} + + +static const ppd_norm pd_norm_tab[7] = { + pd_norm1, + pd_norm2, + pd_norm4, + pd_norm8, + pd_norm16, + pd_norm32, + pd_norm64 +}; + +float pd_norm(float *pd, int nlogdim) +{ + return pd_norm_tab[nlogdim](pd); +} + +void pd_memset(float *dst, const float *src, int ndim, int nitems) +{ + int size = PD_SIZE(ndim); + while(nitems--) { + memcpy(dst,src,size); + dst +=ndim; + } +} + +void pd_fwdperm(float *dst, float *src, const unsigned int *perm, int ndim) +{ + // TODO: non-loop implementation + while (ndim--) + dst[ndim] = src[perm[ndim]]; +} + +void pd_bwdperm(float *dst, float *src, const unsigned int *perm, int ndim) +{ + // TODO: non-loop implementation + while (ndim--) + dst[perm[ndim]] = src[ndim]; +} + +float pd_max(float *src, int ndim) +{ + // TODO: faster implementation + + float cmax=0; // we assume that prob distributions are always positive + float cval; + + while (ndim--) { + cval = src[ndim]; + if (cval>=cmax) { + cmax = cval; + } + } + + return cmax; +} + +int pd_argmax(float *pmax, float *src, int ndim) +{ + // TODO: faster implementation + + float cmax=0; // we assume that prob distributions are always positive + float cval; + int idxmax=-1; // indicates that all pd elements are <0 + + while (ndim--) { + cval = src[ndim]; + if (cval>=cmax) { + cmax = cval; + idxmax = ndim; + } + } + + if (pmax) + *pmax = cmax; + + return idxmax; +} diff --git a/qracodes-mt/pdmath.h b/qracodes-st/pdmath.h similarity index 92% rename from qracodes-mt/pdmath.h rename to qracodes-st/pdmath.h index bd441e0..269501d 100644 --- a/qracodes-mt/pdmath.h +++ b/qracodes-st/pdmath.h @@ -1,85 +1,85 @@ -// pdmath.h -// Elementary math on probability distributions -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ -// This file is part of the qracodes project, a Forward Error Control -// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. -// -// qracodes 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. -// qracodes 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 qracodes source distribution. -// If not, see . - - -#ifndef _pdmath_h_ -#define _pdmath_h_ - -#include - -#ifdef __cplusplus -extern "C" { -#endif - -#define PD_NDIM(nlogdim) ((1<<(nlogdim)) -#define PD_SIZE(ndim) ((ndim)*sizeof(float)) -#define PD_ROWADDR(fp,ndim,idx) (fp+((ndim)*(idx))) - -const float *pd_uniform(int nlogdim); -// Returns a pointer to a (constant) uniform distribution of the given log2 size - -#define pd_init(dst,src,ndim) memcpy(dst,src,PD_SIZE(ndim)) -// Distribution copy - -void pd_memset(float *dst, const float *src, int ndim, int nitems); -// Copy the distribution pointed by src to the array of distributions dst -// src is a pointer to the input distribution (a vector of size ndim) -// dst is a pointer to a linear array of distributions (a vector of size ndim*nitems) - -void pd_imul(float *dst, const float *src, int nlogdim); -// In place multiplication -// Compute dst = dst*src for any element of the distrib give their log2 size -// src and dst arguments must be pointers to array of floats of the given size - -float pd_norm(float *pd, int nlogdim); -// In place normalizazion -// Normalizes the input vector so that the sum of its components are one -// pd must be a pointer to an array of floats of the given size. -// If the norm of the input vector is non-positive the vector components -// are replaced with a uniform distribution -// Returns the norm of the distribution prior to the normalization - -void pd_fwdperm(float *dst, float *src, const int *perm, int ndim); -// Forward permutation of a distribution -// Computes dst[k] = src[perm[k]] for every element in the distribution -// perm must be a pointer to an array of integers of length ndim - -void pd_bwdperm(float *dst, float *src, const int *perm, int ndim); -// Backward permutation of a distribution -// Computes dst[perm[k]] = src[k] for every element in the distribution -// perm must be a pointer to an array of integers of length ndim - -float pd_max(float *src, int ndim); -// Return the maximum of the elements of the given distribution -// Assumes that the input vector is a probability distribution and that each element in the -// distribution is non negative - -int pd_argmax(float *pmax, float *src, int ndim); -// Return the index of the maximum element of the given distribution -// The maximum is stored in the variable pointed by pmax if pmax is not null -// Same note of pd_max applies. -// Return -1 if all the elements in the distribution are negative - -#ifdef __cplusplus -} -#endif - -#endif // _pdmath_h_ +// pdmath.h +// Elementary math on probability distributions +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (repeat and accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + + +#ifndef _pdmath_h_ +#define _pdmath_h_ + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +#define PD_NDIM(nlogdim) ((1<<(nlogdim)) +#define PD_SIZE(ndim) ((ndim)*sizeof(float)) +#define PD_ROWADDR(fp,ndim,idx) (fp+((ndim)*(idx))) + +const float *pd_uniform(int nlogdim); +// Returns a pointer to a (constant) uniform distribution of the given log2 size + +#define pd_init(dst,src,ndim) memcpy(dst,src,PD_SIZE(ndim)) +// Distribution copy + +void pd_memset(float *dst, const float *src, int ndim, int nitems); +// Copy the distribution pointed by src to the array of distributions dst +// src is a pointer to the input distribution (a vector of size ndim) +// dst is a pointer to a linear array of distributions (a vector of size ndim*nitems) + +void pd_imul(float *dst, const float *src, int nlogdim); +// In place multiplication +// Compute dst = dst*src for any element of the distrib give their log2 size +// src and dst arguments must be pointers to array of floats of the given size + +float pd_norm(float *pd, int nlogdim); +// In place normalizazion +// Normalizes the input vector so that the sum of its components are one +// pd must be a pointer to an array of floats of the given size. +// If the norm of the input vector is non-positive the vector components +// are replaced with a uniform distribution +// Returns the norm of the distribution prior to the normalization + +void pd_fwdperm(float *dst, float *src, const unsigned int *perm, int ndim); +// Forward permutation of a distribution +// Computes dst[k] = src[perm[k]] for every element in the distribution +// perm must be a pointer to an array of integers of length ndim + +void pd_bwdperm(float *dst, float *src, const unsigned int *perm, int ndim); +// Backward permutation of a distribution +// Computes dst[perm[k]] = src[k] for every element in the distribution +// perm must be a pointer to an array of integers of length ndim + +float pd_max(float *src, int ndim); +// Return the maximum of the elements of the given distribution +// Assumes that the input vector is a probability distribution and that each element in the +// distribution is non negative + +int pd_argmax(float *pmax, float *src, int ndim); +// Return the index of the maximum element of the given distribution +// The maximum is stored in the variable pointed by pmax if pmax is not null +// Same note of pd_max applies. +// Return -1 if all the elements in the distribution are negative + +#ifdef __cplusplus +} +#endif + +#endif // _pdmath_h_ diff --git a/qracodes-mt/qra12_63_64_irr_b.c b/qracodes-st/qra12_63_64_irr_b.c similarity index 95% rename from qracodes-mt/qra12_63_64_irr_b.c rename to qracodes-st/qra12_63_64_irr_b.c index 1da2e03..5eed61f 100644 --- a/qracodes-mt/qra12_63_64_irr_b.c +++ b/qracodes-st/qra12_63_64_irr_b.c @@ -1,10 +1,10 @@ // qra12_63_64_irr_b.c -// Encoding/Decoding tables for Q-ary RA code (12,63) over GF(64) +// Code tables and defines for the Q-ary irregular RA code (12,63) over GF(64) // Code Name: qra12_63_64_irr_b -// (12,63) RA Code over GF(64) - RF=333344455567 - -// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy - +// Systematic symbols repetition factors: 333344455567 +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. // @@ -16,18 +16,17 @@ // 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 qracodes source distribution. +// along with qra_codes source distribution. // If not, see . #include "qra12_63_64_irr_b.h" - +/* #define qra_K 12 // number of information symbols #define qra_N 63 // codeword length in symbols #define qra_m 6 // bits/symbol #define qra_M 64 // Symbol alphabet cardinality -#define qra_a 1 // grouping factor #define qra_NC 51 // number of check symbols (N-K) // Defines used by the message passing decoder -------- @@ -35,14 +34,14 @@ #define qra_V 63 // number of variables in the code graph (N) #define qra_C 115 // number of factors in the code graph (N +(N-K)+1) #define qra_NMSG 217 // number of msgs in the code graph -#define qra_MAXVDEG 8 // maximum variable degree +#define qra_MAXVDEG 8 // maximum variable degree (intrinsic check included) #define qra_MAXCDEG 3 // maximum factor degree -#define qra_R 0.19048f // code rate (K/N) -#define CODE_NAME "qra_12_63_64_irr_b" - +#define qra_R 0.19048 // code rate (K/N) +*/ +// Tables used by the encoder ------------------------- // table of the systematic symbols indexes in the accumulator chain -static const int qra_acc_input_idx[qra_NC+1] = { +const unsigned int qra_acc_input_idx[qra_NC+1] = { 3, 11, 0, 1, 7, 8, 6, 5, 10, 4, 11, 9, 0, 2, 6, 7, 8, 4, 11, 5, 10, 2, 1, 9, 3, 8, 4, 11, 5, 7, @@ -52,7 +51,7 @@ static const int qra_acc_input_idx[qra_NC+1] = { }; // table of the systematic symbols weight logarithms over GF(M) -static const int qra_acc_input_wlog[qra_NC+1] = { +const unsigned int qra_acc_input_wlog[qra_NC+1] = { 39, 0, 34, 16, 25, 0, 34, 48, 19, 13, 29, 56, 0, 5, 39, 42, 31, 0, 10, 0, 57, 62, 33, 43, 0, 14, 22, 48, 28, 20, @@ -62,7 +61,7 @@ static const int qra_acc_input_wlog[qra_NC+1] = { }; // table of the logarithms of the elements of GF(M) (log(0) never used) -static const int qra_log[qra_M] = { +const unsigned int qra_log[qra_M] = { -1, 0, 1, 6, 2, 12, 7, 26, 3, 32, 13, 35, 8, 48, 27, 18, 4, 24, 33, 16, 14, 52, 36, 54, 9, 45, 49, 38, 28, 41, @@ -73,7 +72,7 @@ static const int qra_log[qra_M] = { }; // table of GF(M) elements given their logarithm -static const int qra_exp[qra_M-1] = { +const unsigned int qra_exp[qra_M-1] = { 1, 2, 4, 8, 16, 32, 3, 6, 12, 24, 48, 35, 5, 10, 20, 40, 19, 38, 15, 30, 60, 59, 53, 41, 17, 34, 7, 14, 28, 56, @@ -83,8 +82,10 @@ static const int qra_exp[qra_M-1] = { 57, 49, 33 }; +// Tables used by the decoder ------------------------- + // table of the messages weight logarithms over GF(M) -static const int qra_msgw[qra_NMSG] = { +const unsigned int qra_msgw[qra_NMSG] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, @@ -110,7 +111,7 @@ static const int qra_msgw[qra_NMSG] = { }; // table of the degrees of the variable nodes -static const int qra_vdeg[qra_V] = { +const unsigned int qra_vdeg[qra_V] = { 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, @@ -121,7 +122,7 @@ static const int qra_vdeg[qra_V] = { }; // table of the degrees of the factor nodes -static const int qra_cdeg[qra_C] = { +const unsigned int qra_cdeg[qra_C] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, @@ -137,7 +138,7 @@ static const int qra_cdeg[qra_C] = { }; // table (uncompressed) of the v->c message indexes (-1=unused entry) -static const int qra_v2cmidx[qra_V*qra_MAXVDEG] = { +const unsigned int qra_v2cmidx[qra_V*qra_MAXVDEG] = { 0, 65, 75, 101, -1, -1, -1, -1, 1, 66, 85, 110, -1, -1, -1, -1, 2, 76, 84, 106, -1, -1, -1, -1, @@ -204,7 +205,7 @@ static const int qra_v2cmidx[qra_V*qra_MAXVDEG] = { }; // table (uncompressed) of the c->v message indexes (-1=unused entry) -static const int qra_c2vmidx[qra_C*qra_MAXCDEG] = { +const unsigned int qra_c2vmidx[qra_C*qra_MAXCDEG] = { 0, -1, -1, 1, -1, -1, 2, -1, -1, 3, -1, -1, 4, -1, -1, 5, -1, -1, 6, -1, -1, 7, -1, -1, 8, -1, -1, 9, -1, -1, 10, -1, -1, 11, -1, -1, @@ -237,7 +238,7 @@ static const int qra_c2vmidx[qra_C*qra_MAXCDEG] = { }; // permutation matrix to compute Prob(x*alfa^logw) -static const int qra_pmat[qra_M*qra_M] = { +const unsigned int qra_pmat[qra_M*qra_M] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, @@ -492,43 +493,5 @@ static const int qra_pmat[qra_M*qra_M] = { 35, 33, 39, 37, 43, 41, 47, 45, 51, 49, 55, 53, 59, 57, 63, 61 }; -const qracode qra_12_63_64_irr_b = { - qra_K, - qra_N, - qra_m, - qra_M, - qra_a, - qra_NC, - qra_V, - qra_C, - qra_NMSG, - qra_MAXVDEG, - qra_MAXCDEG, - QRATYPE_NORMAL, - qra_R, - CODE_NAME, - qra_acc_input_idx, - qra_acc_input_wlog, - qra_log, - qra_exp, - qra_msgw, - qra_vdeg, - qra_cdeg, - qra_v2cmidx, - qra_c2vmidx, - qra_pmat -}; -#undef qra_K -#undef qra_N -#undef qra_m -#undef qra_M -#undef qra_a -#undef qra_NC -#undef qra_V -#undef qra_C -#undef qra_NMSG -#undef qra_MAXVDEG -#undef qra_MAXCDEG -#undef qra_R -#undef CODE_NAME + diff --git a/qracodes-st/qra12_63_64_irr_b.h b/qracodes-st/qra12_63_64_irr_b.h new file mode 100644 index 0000000..662078e --- /dev/null +++ b/qracodes-st/qra12_63_64_irr_b.h @@ -0,0 +1,64 @@ +// qra12_63_64_irr_b.h +// Code tables and defines for the Q-ary irregular RA code (12,63) over GF(64) +// Code Name: qra12_63_64_irr_b +// Systematic symbols repetition factors: 333344455567 +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ +// This file is part of the qracodes project, a Forward Error Control +// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. +// +// qracodes 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. +// qracodes 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 qra_codes source distribution. +// If not, see . + +#ifndef _qra12_64_64_irr_b_h +#define _qra12_64_64_irr_b_h + + +#define qra_K 12 // number of information symbols +#define qra_N 63 // codeword length in symbols +#define qra_m 6 // bits/symbol +#define qra_M 64 // Symbol alphabet cardinality +#define qra_NC 51 // number of check symbols (N-K) + +// Defines used by the message passing decoder -------- + +#define qra_V 63 // number of variables in the code graph (N) +#define qra_C 115 // number of factors in the code graph (N +(N-K)+1) +#define qra_NMSG 217 // number of msgs in the code graph +#define qra_MAXVDEG 8 // maximum variable degree (intrinsic check included) +#define qra_MAXCDEG 3 // maximum factor degree +#define qra_R 0.19048 // code rate (K/N) + +#ifdef __cplusplus +extern "C" { +#endif + +// Tables used by the encoder ------------------------- +extern const unsigned int qra_acc_input_idx[qra_NC+1]; +extern const unsigned int qra_acc_input_wlog[qra_NC+1]; +extern const unsigned int qra_log[qra_M]; +extern const unsigned int qra_exp[qra_M-1]; +// Tables used by the decoder ------------------------- +extern const unsigned int qra_msgw[qra_NMSG]; +extern const unsigned int qra_vdeg[qra_V]; +extern const unsigned int qra_cdeg[qra_C]; +extern const unsigned int qra_v2cmidx[qra_V*qra_MAXVDEG]; +extern const unsigned int qra_c2vmidx[qra_C*qra_MAXCDEG]; +extern const unsigned int qra_pmat[qra_M*qra_M]; + +#ifdef __cplusplus +} +#endif + +#endif // _qra12_64_64_irr_b_h diff --git a/qracodes-mt/qracodes-mt.vcproj b/qracodes-st/qracodes-st.vcproj similarity index 87% rename from qracodes-mt/qracodes-mt.vcproj rename to qracodes-st/qracodes-st.vcproj index a98d825..c1949e7 100644 --- a/qracodes-mt/qracodes-mt.vcproj +++ b/qracodes-st/qracodes-st.vcproj @@ -2,9 +2,10 @@ @@ -40,9 +41,13 @@ @@ -57,7 +62,9 @@ /> - - @@ -222,10 +220,6 @@ RelativePath=".\qra12_63_64_irr_b.h" > - - diff --git a/qracodes-st/qracodes-st.vcproj.NICOW7.User.user b/qracodes-st/qracodes-st.vcproj.NICOW7.User.user new file mode 100644 index 0000000..89376bb --- /dev/null +++ b/qracodes-st/qracodes-st.vcproj.NICOW7.User.user @@ -0,0 +1,65 @@ + + + + + + + + + + + diff --git a/qracodes-mt/qracodes.c b/qracodes-st/qracodes.c similarity index 73% rename from qracodes-mt/qracodes.c rename to qracodes-st/qracodes.c index fc80dcc..ae96ac0 100644 --- a/qracodes-mt/qracodes.c +++ b/qracodes-st/qracodes.c @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qracodes source distribution. +// along with qra_codes source distribution. // If not, see . #include @@ -26,97 +26,55 @@ #include "pdmath.h" #include "qracodes.h" +#include "qra12_63_64_irr_b.h" // tables for the P(12,63) Q-ary RA code over GF(64) -int qra_encode(const qracode *pcode, int *y, const int *x) +// #define QRA_DEBUG + +// message tables used by the qra_extrinsic +// table of the v->c messages +static float qra_v2cmsg[qra_NMSG*qra_M]; +// table of the c->v messages +static float qra_c2vmsg[qra_NMSG*qra_M]; + +int qra_encode(uint *y, const uint *x) { - int k,j,kk,jj; - int t, chk = 0; - - const int K = pcode->K; - const int M = pcode->M; - const int NC= pcode->NC; - const int a = pcode->a; - const int *acc_input_idx = pcode->acc_input_idx; - const int *acc_input_wlog = pcode->acc_input_wlog; - const int *gflog = pcode->gflog; - const int *gfexp = pcode->gfexp; + int k; + uint t, chk = 0; // copy the systematic symbols to destination - memcpy(y,x,K*sizeof(int)); + memcpy(y,x,qra_K*sizeof(uint)); - y = y+K; // point to check symbols + y = y+qra_K; // point to check symbols // compute the code check symbols as a weighted accumulation of a permutated // sequence of the (repeated) systematic input symbols: // chk(k+1) = x(idx(k))*alfa^(logw(k)) + chk(k) - // (all operations performed over GF(M)) - - if (a==1) { // grouping factor = 1 - for (k=0;k 1 - for (k=0;kM; - const int qra_m = pcode->m; - const int qra_V = pcode->V; - const int qra_MAXVDEG = pcode->MAXVDEG; - const int *qra_vdeg = pcode->vdeg; - const int qra_C = pcode->C; - const int qra_MAXCDEG = pcode->MAXCDEG; - const int *qra_cdeg = pcode->cdeg; - const int *qra_v2cmidx = pcode->v2cmidx; - const int *qra_c2vmidx = pcode->c2vmidx; - const int *qra_pmat = pcode->gfpmat; - const int *qra_msgw = pcode->msgw; - -// float msgout[qra_M]; // buffer to store temporary results - float msgout[QRACODE_MAX_M]; // we use a fixed size in order to avoid mallocs - - float totex; // total extrinsic information + // TODO: pass the function pointers to message tables + // so that multiple instances of the function can be run on different threads (and/or CPU cores) + + float msgout[qra_M]; // buffer to store temporary results + float totex; // total extrinsic information int nit; // current iteration - int nv; // current variable - int nc; // current check + int nv; // current variable + int nc; // current check int k,kk; // loop indexes int ndeg; // current node degree @@ -254,12 +187,8 @@ int qra_extrinsic(const qracode *pcode, int rc = -1; // rc>=0 extrinsic converged to 1 at iteration rc (rc=0..maxiter-1) // rc=-1 no convergence in the given number of iterations // rc=-2 error in the code tables (code checks degrees must be >1) - // rc=-3 M is larger than QRACODE_MAX_M - - if (qra_M>QRACODE_MAX_M) - return -3; // message initialization ------------------------------------------------------- @@ -446,7 +375,7 @@ int qra_extrinsic(const qracode *pcode, return rc; } -void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix) +void qra_mapdecode(uint *xdec, float *pex, const float *pix) { // Maximum a posteriori probability decoding. // Given the intrinsic information (pix) and extrinsic information (pex) (computed with qra_extrinsic(...)) @@ -456,11 +385,7 @@ void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix // Returns: // xdec[k] = decoded (information) symbols k=[0..qra_K-1] -// Note: pex is destroyed and overwritten with mapp - - const int qra_M = pcode->M; - const int qra_m = pcode->m; - const int qra_K = pcode->K; +// Note: pex is destroyed and overwritten with pmap int k; diff --git a/qracodes-mt/qra12_63_64_irr_b.h b/qracodes-st/qracodes.h similarity index 60% rename from qracodes-mt/qra12_63_64_irr_b.h rename to qracodes-st/qracodes.h index bf70b36..8503252 100644 --- a/qracodes-mt/qra12_63_64_irr_b.h +++ b/qracodes-st/qracodes.h @@ -1,10 +1,8 @@ -// qra12_63_64_irr_b.h -// Code tables and defines for Q-ary RA code (12,63) over GF(64) -// Code Name: qra12_63_64_irr_b -// (12,63) RA Code over GF(64) - RF=333344455567 - -// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy - +// qracodes.h +// Q-ary RA codes encoding/decoding functions +// +// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// ------------------------------------------------------------------------------ // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. // @@ -16,24 +14,27 @@ // 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 qracodes source distribution. +// along with qra_codes source distribution. // If not, see . -#ifndef _qra12_63_64_irr_b_h -#define _qra12_63_64_irr_b_h +#ifndef _qracodes_h_ +#define _qracodes_h_ -#include "qracodes.h" +typedef unsigned int uint; #ifdef __cplusplus extern "C" { #endif -extern const qracode qra_12_63_64_irr_b; +int qra_encode(uint *y, const uint *x); +void qra_mfskbesselmetric(float *pix, float EsNoMetric); +int qra_extrinsic(float *pex, const float *pix, int maxiter); +void qra_mapdecode(uint *xdec, float *pex, const float *pix); #ifdef __cplusplus } #endif -#endif // _qra12_63_64_irr_b_h +#endif // _qracodes_h_ diff --git a/qracodes.sln b/qracodes.sln index fd2d6fc..4614cb7 100644 --- a/qracodes.sln +++ b/qracodes.sln @@ -1,7 +1,7 @@  Microsoft Visual Studio Solution File, Format Version 10.00 # Visual Studio 2008 -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "qracodes", "qracodes\qracodes.vcproj", "{8540B348-B8F2-4104-84CB-6286AC19158B}" +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "qracodes", "qracodes\qracodes.vcproj", "{933A58F6-199B-4723-ACFE-3013E6DD9D0A}" EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution @@ -9,10 +9,10 @@ Global Release|Win32 = Release|Win32 EndGlobalSection GlobalSection(ProjectConfigurationPlatforms) = postSolution - {8540B348-B8F2-4104-84CB-6286AC19158B}.Debug|Win32.ActiveCfg = Debug|Win32 - {8540B348-B8F2-4104-84CB-6286AC19158B}.Debug|Win32.Build.0 = Debug|Win32 - {8540B348-B8F2-4104-84CB-6286AC19158B}.Release|Win32.ActiveCfg = Release|Win32 - {8540B348-B8F2-4104-84CB-6286AC19158B}.Release|Win32.Build.0 = Release|Win32 + {933A58F6-199B-4723-ACFE-3013E6DD9D0A}.Debug|Win32.ActiveCfg = Debug|Win32 + {933A58F6-199B-4723-ACFE-3013E6DD9D0A}.Debug|Win32.Build.0 = Debug|Win32 + {933A58F6-199B-4723-ACFE-3013E6DD9D0A}.Release|Win32.ActiveCfg = Release|Win32 + {933A58F6-199B-4723-ACFE-3013E6DD9D0A}.Release|Win32.Build.0 = Release|Win32 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/qracodes.suo b/qracodes.suo index d49a136..e2128b3 100644 Binary files a/qracodes.suo and b/qracodes.suo differ diff --git a/qracodes/build.sh b/qracodes/build.sh index f8dceb5..2792635 100644 --- a/qracodes/build.sh +++ b/qracodes/build.sh @@ -1,3 +1,2 @@ -g++ -march=native -O3 -DFTZ_ENABLE -DQRA_DEBUG -DTEST_WER_AWGN *.c -o qracodesawgn && ./qracodesawgn -g++ -march=native -O3 -DFTZ_ENABLE -DQRA_DEBUG -DTEST_WER_RAYLEIGH *.c -o qracodesrayleigh && ./qracodesrayleigh +gcc -Wall -march=native -pthread -O3 *.c -lpthread -lm -o qracodes diff --git a/qracodes-mt/ebnovalues.txt b/qracodes/ebnovalues.txt similarity index 93% rename from qracodes-mt/ebnovalues.txt rename to qracodes/ebnovalues.txt index a26d7ef..7dba138 100644 --- a/qracodes-mt/ebnovalues.txt +++ b/qracodes/ebnovalues.txt @@ -1,15 +1,15 @@ -# Eb/No Values to be used during the code simulation -# Each line of this file indicates the Eb/No value to be simulated (in dB) -# and the number of errors to be detected by the decoder -1.1 1000 -1.6 1000 -2.1 1000 -2.6 1000 -3.1 1000 -3.6 1000 -4.1 1000 -4.6 1000 -5.1 500 -5.6 200 -6.1 100 +# Eb/No Values to be used during the code simulation +# Each line of this file indicates the Eb/No value to be simulated (in dB) +# and the number of errors to be detected by the decoder +1.1 1000 +1.6 1000 +2.1 1000 +2.6 1000 +3.1 1000 +3.6 1000 +4.1 1000 +4.6 1000 +5.1 500 +5.6 200 +6.1 100 6.6 50 \ No newline at end of file diff --git a/qracodes-mt/ebnovaluesfast.txt b/qracodes/ebnovaluesfast.txt similarity index 95% rename from qracodes-mt/ebnovaluesfast.txt rename to qracodes/ebnovaluesfast.txt index 84d030a..b057b31 100644 --- a/qracodes-mt/ebnovaluesfast.txt +++ b/qracodes/ebnovaluesfast.txt @@ -1,11 +1,11 @@ -# Eb/No Values to be used during the code simulation -# Each line of this file indicates the Eb/No value to be simulated (in dB) -# and the number of errors to be detected by the decoder -1.1 500 -1.6 500 -2.1 500 -2.6 500 -3.1 500 -3.6 500 -4.1 200 -4.6 50 +# Eb/No Values to be used during the code simulation +# Each line of this file indicates the Eb/No value to be simulated (in dB) +# and the number of errors to be detected by the decoder +1.1 500 +1.6 500 +2.1 500 +2.6 500 +3.1 500 +3.6 500 +4.1 200 +4.6 50 diff --git a/qracodes/main.c b/qracodes/main.c index 3a4488c..554736e 100644 --- a/qracodes/main.c +++ b/qracodes/main.c @@ -1,18 +1,24 @@ // main.c -// Word Error Rate test for the Q-ary RA code (12,63) over GF(64) +// Word Error Rate test example for Q-ary RA codes over GF(64) +// Multi-threaded simulator version + +// (c) 2016 - Nico Palermo, IV3NWV // -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// Thanks to Andrea Montefusco IW0HDV for his help on adapting the sources +// to OSs other than MS Windows +// // ------------------------------------------------------------------------------ // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. // -// Files in the package: +// Files in this package: // main.c - this file -// normrnd.c/.h - random gaussian number generator and header +// normrnd.c/.h - random gaussian number generator // npfwht.c/.h - Fast Walsh-Hadamard Transforms // pdmath.c/.h - Elementary math on probability distributions -// qra12_63_64_irr_b.c/.h - Tables and defines for a P(12,63) IRA code over GF(64) -// qracodes.c/.h - QRA (Q-ary Repeat & Accumulate) codes encoding/decoding functions +// qra12_63_64_irr_b.c/.h - Tables for a QRA(12,63) irregular RA code over GF(64) +// qra13_64_64_irr_e.c/.h - Tables for a QRA(13,64) irregular RA code " " +// qracodes.c/.h - QRA codes encoding/decoding functions // // ------------------------------------------------------------------------------- // @@ -26,96 +32,274 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . +// ----------------------------------------------------------------------------- + +// Two codes are available for simulations in this sowftware release: + +// QRA12_63_64_IRR_B: K=12 N=63 Q=64 irregular QRA code (defined in qra12_63_64_irr_b.h /.c) +// QRA13_64_64_IRR_E: K=13 N=64 Q=64 irregular QRA code (defined in qra13_64_64_irr_b.h /.c) + +// Codes with K=13 are designed to include a CRC as the 13th information symbol +// and improve the code UER (Undetected Error Rate). +// The CRC symbol is not sent along the channel (the codes are punctured) and the +// resulting code is still a (12,63) code with an effective code rate of R = 12/63. + +// ------------------------------------------------------------------------------ + +// OS dependent defines and includes -------------------------------------------- + #if _WIN32 // note the underscore: without it, it's not msdn official! // Windows (x64 and x86) #include // required only for GetTickCount(...) + #include // _beginthread #endif -#if __linux__ -#include +#if defined(__linux__) + +// remove unwanted macros +#define __cdecl + +// implements Windows API #include -unsigned GetTickCount() { + unsigned int GetTickCount(void) { struct timespec ts; - unsigned theTick = 0U; + unsigned int theTick = 0U; clock_gettime( CLOCK_REALTIME, &ts ); theTick = ts.tv_nsec / 1000000; theTick += ts.tv_sec * 1000; return theTick; } + +// Convert Windows millisecond sleep +// +// VOID WINAPI Sleep(_In_ DWORD dwMilliseconds); +// +// to Posix usleep (in microseconds) +// +// int usleep(useconds_t usec); +// +#include +#define Sleep(x) usleep(x*1000) + +#endif + +#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) +#include #endif #if __APPLE__ #endif +#include #include -#ifdef FTZ_ENABLE // shouldn't be required anymore - kept here just as a remind in the case we need it - // round off errors in the qra_extrinsics fixed with biasing -#include // used to set the Flush to Zero flag in the FPU -#endif +#include "qracodes.h" +#include "normrnd.h" // gaussian numbers generator +#include "pdmath.h" // operations on probability distributions + +// defined codes +#include "qra12_63_64_irr_b.h" +#include "qra13_64_64_irr_e.h" -#include "normrnd.h" // gaussian number generators -#include "qra12_63_64_irr_b.h" // tables and defines for the P(12,63) Q-ary RA code over GF(64) -#include "qracodes.h" // encoding/decoding functions +// ----------------------------------------------------------------------------------- -#define NSAMPLES (qra_N*qra_M) // samples per message (codeword length * codeword alphabet size) +#define NTHREADS_MAX 160 +// channel types #define CHANNEL_AWGN 0 #define CHANNEL_RAYLEIGH 1 -int wer_test(int channel_type, float EbNodB, int desirederrs, FILE *fout) +// amount of a-priori information provided to the decoder +#define AP_NONE 0 +#define AP_28 1 +#define AP_44 2 +#define AP_56 3 + +const char ap_str[4][16] = { + "None", + "28 bit", + "44 bit", + "56 bit" +}; + +const char fnameout_pfx[2][64] = { + "wer-awgn-", + "wer-rayleigh-" +}; +const char fnameout_sfx[4][64] = { + "-ap00.txt", + "-ap28.txt", + "-ap44.txt", + "-ap56.txt" +}; + +const int ap_masks_jt65[4][13] = { +// Each row must be 13 entries long (to handle puntc. codes 13,64) +// The mask of 13th symbol (crc) is alway initializated to 0 + // AP0 - no a-priori knowledge + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + // AP28 - 1st field known [cq ? ?] or [dst ? ?] + {0x3F,0x3F,0x3F,0x3F,0x3C, 0, 0, 0, 0, 0, 0, 0}, + // AP44 - 1st and 3rd fields known [cq ? 0] or [dst ? 0] + {0x3F,0x3F,0x3F,0x3F,0x3C, 0, 0, 0, 0,0x0F,0x3F,0x3F}, + // AP56 - 1st and 2nd fields known [dst src ?] + {0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x3F,0x30, 0, 0} +}; + +void ix_mask(const qracode *pcode, float *r, const int *mask, const int *x); + +void printword(char *msg, int *x, int size) { int k; + printf("\n%s ",msg); + for (k=0;kc msg buffer + float *qra_c2vmsg; //[qra_NMSG*qra_M]; MP decoder c->v msg buffer + float *rp; // [qra_N*qra_M]; received samples (real component) buffer + float *rq; // [qra_N*qra_M]; received samples (imag component) buffer + float *chp; //[qra_N]; channel gains (real component) buffer + float *chq; //[qra_N]; channel gains (imag component) buffer + float *r; //[qra_N*qra_M]; received samples (amplitude) buffer + float *ix; // [qra_N*qra_M]; // intrinsic information to the MP algorithm + float *ex; // [qra_N*qra_M]; // extrinsic information from the MP algorithm +#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) + pthread_t thread; +#endif +} wer_test_ds; + +typedef void( __cdecl *pwer_test_thread)(wer_test_ds*); -// uint x[qra_K]= {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}; - uint x[qra_K]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,}; - uint y[qra_N],ydec[qra_N]; +// crc-6 generator polynomial +// g(x) = x^6 + a5*x^5 + ... + a1*x + a0 - float rp[qra_N*qra_M]; // received samples (real component) - float rq[qra_N*qra_M]; // received samples (imag component) - float chp[qra_N]; // channel gains (real component) - float chq[qra_N]; // channel gains (imag component) - float r[qra_N*qra_M]; // r = sqrt(rp^2+rq^2) +// g(x) = x^6 + x + 1 +#define CRC6_GEN_POL 0x30 // MSB=a0 LSB=a5 - float No = 1.0f; // noise spectral density - float sigma = (float)sqrt(No/2.0f); // std dev of noise I/Q components - float sigmach = (float)sqrt(1/2.0f); // std dev of channel I/Q gains - float R = 1.0f*qra_K/qra_N; // code rate +// g(x) = x^6 + x^2 + x + 1 (as suggested by Joe. See: https://users.ece.cmu.edu/~koopman/crc/) +// #define CRC6_GEN_POL 0x38 // MSB=a0 LSB=a5. Simulation results are similar + +int calc_crc6(int *x, int sz) +{ + int k,j,t,sr = 0; + for (k=0;k>1) ^ CRC6_GEN_POL; + else + sr = (sr>>1); + t>>=1; + } + } + return sr; +} - float EbNo = (float)pow(10,EbNodB/10); - float EsNo = 1.0f*qra_m*R*EbNo; - float Es = EsNo*No; - float A = (float)sqrt(Es); +void wer_test_thread(wer_test_ds *pdata) +{ + const qracode *pcode=pdata->pcode; + const int qra_K = pcode->K; + const int qra_N = pcode->N; + const int qra_M = pcode->M; + const int qra_m = pcode->m; + const int NSAMPLES = pcode->N*pcode->M; + + const float No = 1.0f; // noise spectral density + const float sigma = (float)sqrt(No/2.0f); // std dev of noise I/Q components + const float sigmach = (float)sqrt(1/2.0f); // std dev of channel I/Q gains // Eb/No value for which we optimize the bessel metric - const float EbNodBMetric = 2.8f; + const float EbNodBMetric = 2.8f; const float EbNoMetric = (float)pow(10,EbNodBMetric/10); - const float EsNoMetric = 1.0f*qra_m*R*EbNoMetric; - float ex[qra_N*qra_M]; // extrinsic information from the MP algorithm - float avgt; // average time per codeword encoding/channel sim/decoding + int k,t,diff; + + float R; + float EsNoMetric; + float EbNo, EsNo, Es, A; + int channel_type, code_type; int nt=0; // transmitted codewords int nerrs = 0; // total number of errors int nerrsu = 0; // number of undetected errors - int nep = -1; + int rc; - unsigned int cini,cend; // time tick counters - int rc; + // inizialize pointer to required buffers + int *x=pdata->x; // message buffer + int *y=pdata->y, *ydec=pdata->ydec; // encoded/decoded codeword buffers + float *qra_v2cmsg=pdata->qra_v2cmsg; // table of the v->c messages + float *qra_c2vmsg=pdata->qra_c2vmsg; // table of the c->v messages + float *rp=pdata->rp; // received samples (real component) + float *rq=pdata->rq; // received samples (imag component) + float *chp=pdata->chp; // channel gains (real component) + float *chq=pdata->chq; // channel gains (imag component) + float *r=pdata->r; // received samples amplitudes + float *ix=pdata->ix; // intrinsic information to the MP algorithm + float *ex=pdata->ex; // extrinsic information from the MP algorithm + + channel_type = pdata->channel_type; + code_type = pcode->type; + + // define the (true) code rate accordingly to the code type + switch(code_type) { + case QRATYPE_CRC: + R = 1.0f*(qra_K-1)/qra_N; + break; + case QRATYPE_CRCPUNCTURED: + R = 1.0f*(qra_K-1)/(qra_N-1); + break; + case QRATYPE_NORMAL: + default: + R = 1.0f*(qra_K)/(qra_N); + } - // encode the input - qra_encode(y,x); + EsNoMetric = 1.0f*qra_m*R*EbNoMetric; + + EbNo = (float)pow(10,pdata->EbNodB/10); + EsNo = 1.0f*qra_m*R*EbNo; + Es = EsNo*No; + A = (float)sqrt(Es); - cini = GetTickCount(); + // encode the input + if (code_type==QRATYPE_CRC || code_type==QRATYPE_CRCPUNCTURED) { + // compute the information message symbol check as the (negated) xor of all the + // information message symbols + for (k=0;k<(qra_K-1);k++) + x[k]=k%qra_M; + x[k]=calc_crc6(x,qra_K-1); + } + else + for (k=0;kstop==0) { // simulate the channel + // NOTE: in the case that the code is punctured, for simplicity + // we compute the channel outputs and the metric also for the crc symbol + // then we ignore its observation. normrnd_s(rp,NSAMPLES,0,sigma); normrnd_s(rq,NSAMPLES,0,sigma); @@ -131,138 +315,436 @@ int wer_test(int channel_type, float EbNodB, int desirederrs, FILE *fout) rq[k*qra_M+y[k]]+=A*chq[k]; } } - else - return -1; // unknown channel type + else { + pdata->done = 1; + return; // unknown channel type + } // compute the squares of the amplitudes of the received samples for (k=0;km,pcode->N,EsNoMetric); + + if (code_type==QRATYPE_CRCPUNCTURED) { + // ignore observations of the CRC symbol as it is not actually sent + // over the channel + pd_init(PD_ROWADDR(ix,qra_M,qra_K),pd_uniform(qra_m),qra_M); + } + + + if (pdata->ap_index!=0) + // mask channel observations with a priori knowledge + ix_mask(pcode,ix,ap_masks_jt65[pdata->ap_index],x); + // compute the extrinsic symbols probabilities with the message-passing algorithm - // stop if extrinsic information does not converges to 1 within 50 iterations - rc = qra_extrinsic(ex,r,50); + // stop if extrinsic information does not converges to 1 within the given number of iterations + rc = qra_extrinsic(pcode,ex,ix,100,qra_v2cmsg,qra_c2vmsg); if (rc>=0) { // the MP algorithm converged to Iex~1 in rc iterations // decode the codeword - qra_mapdecode(ydec,ex,r); + qra_mapdecode(pcode,ydec,ex,ix); // look for undetected errors - for (k=0;knt=nt; + pdata->nerrs=nerrs; + pdata->nerrsu=nerrsu; } - cend = GetTickCount(); - printf("\n"); fflush (stdout); + pdata->done=1; - avgt = 1.0f*(cend-cini)/nt; // average time per decode(ms) - printf("Elapsed time=%5.1fs Avg time per word=%6.2fms ntx=%6d errs=%5d errsu=%3d pe=%.2e\n",0.001f*(cend-cini),avgt, nt, nerrs, nerrsu, 1.0f*nerrs/nt); - fflush (stdout); - fprintf(fout,"%01d %.2f %d %d %d %.2f\n",channel_type, EbNodB,nt,nerrs,nerrsu, avgt); + #if _WIN32 + _endthread(); + #endif +} + +#if defined(__linux__) || ( defined(__MINGW32__) || defined (__MIGW64__) ) +void *wer_test_pthread(void *p) +{ + wer_test_thread ((wer_test_ds *)p); return 0; } -void wer_test_awgn(void) +#endif + + +void ix_mask(const qracode *pcode, float *r, const int *mask, const int *x) { - int k,nn; + // mask intrinsic information (channel observations) with a priori knowledge + + int k,kk, smask; + const int qra_K=pcode->K; + const int qra_M=pcode->M; + const int qra_m=pcode->m; + + for (k=0;kNTHREADS_MAX) { + printf("Error: nthreads should be <=%d\n",NTHREADS_MAX); + return -1; + } + + sprintf(fnameout,"%s%s%s", + fnameout_pfx[chtype], + pcode->name, + fnameout_sfx[ap_index]); - nn = sizeof(EbNodB)/sizeof(float); + fout = fopen(fnameout,"w"); + fprintf(fout,"# Channel (0=AWGN,1=Rayleigh), Eb/No (dB), Transmitted codewords, Errors, Undetected Errors, Avg dec. time (ms), WER\n"); - printf("\nTesting the code over the AWGN channel\n"); - fout = fopen("wer-awgn.txt","w"); + printf("\nTesting the code %s over the %s channel\nSimulation data will be saved to %s\n", + pcode->name, + chtype==CHANNEL_AWGN?"AWGN":"Rayleigh", + fnameout); + fflush (stdout); - for (k=0;kK*sizeof(int)); + wt[j].y = (int*)malloc(pcode->N*sizeof(int)); + wt[j].ydec = (int*)malloc(pcode->N*sizeof(int)); + wt[j].qra_v2cmsg = (float*)malloc(pcode->NMSG*pcode->M*sizeof(float)); + wt[j].qra_c2vmsg = (float*)malloc(pcode->NMSG*pcode->M*sizeof(float)); + wt[j].rp = (float*)malloc(pcode->N*pcode->M*sizeof(float)); + wt[j].rq = (float*)malloc(pcode->N*pcode->M*sizeof(float)); + wt[j].chp = (float*)malloc(pcode->N*sizeof(float)); + wt[j].chq = (float*)malloc(pcode->N*sizeof(float)); + wt[j].r = (float*)malloc(pcode->N*pcode->M*sizeof(float)); + wt[j].ix = (float*)malloc(pcode->N*pcode->M*sizeof(float)); + wt[j].ex = (float*)malloc(pcode->N*pcode->M*sizeof(float)); } - fclose(fout); -} + for (k=0;k=nerrstgt[k]) { + for (j=0;j] [-t] [-c] [-a] [-f[-h]\n"); + printf("Options: \n"); + printf(" -q: code to simulate. 0=qra_12_63_64_irr_b\n"); + printf(" 1=qra_13_64_64_irr_e (default)\n"); + printf(" -t : number of threads to be used for the simulation [1..24]\n"); + printf(" (default=8)\n"); + printf(" -c : channel_type. 0=AWGN 1=Rayleigh \n"); + printf(" (default=AWGN)\n"); + printf(" -a : amount of a-priori information provided to decoder. \n"); + printf(" 0= No a-priori (default)\n"); + printf(" 1= 28 bit \n"); + printf(" 2= 44 bit \n"); + printf(" 3= 56 bit \n"); + printf(" -f : name of the file containing the Eb/No values to be simulated\n"); + printf(" (default=ebnovalues.txt)\n"); + printf(" This file should contain lines in this format:\n"); + printf(" # Eb/No(dB) Target Errors\n"); + printf(" 0.1 5000\n"); + printf(" 0.6 5000\n"); + printf(" 1.1 1000\n"); + printf(" 1.6 1000\n"); + printf(" ...\n"); + printf(" (lines beginning with a # are treated as comments\n\n"); + + printf(" sizeof(unsigned int) = %lu bytes\n", (unsigned long) sizeof(unsigned int)); + printf(" sizeof(unsigned long) = %lu bytes\n", (unsigned long) sizeof(unsigned long)); + printf(" sizeof(unsigned long long) = %lu bytes\n", (unsigned long) sizeof(unsigned long long)); + printf(" sizeof(unsigned float) = %lu bytes\n", (unsigned long) sizeof(float)); + printf(" sizeof(unsigned double) = %lu bytes\n", (unsigned long) sizeof(double)); + printf(" sizeof(void *) = %lu bytes\n", (unsigned long) sizeof(void *)); + + printf("\n\n"); +} + +#define SIM_POINTS_MAX 20 int main(int argc, char* argv[]) { -#ifdef FTZ_ENABLE - // set the Flush-to-Zero flag in the FPU -#if _WIN32 - int ftzmode; - _MM_GET_FLUSH_ZERO_MODE(ftzmode); -#elif __linux__ - int ftzmode = - _MM_GET_FLUSH_ZERO_MODE(); -#else -#error Platform does not support _MM_GET_FLUSH_ZERO_MODE -#endif - _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); -#endif -#ifdef TEST_WER_AWGN - wer_test_awgn(); -#endif + float EbNodB[SIM_POINTS_MAX]; + int nerrstgt[SIM_POINTS_MAX]; + FILE *fin; + + char fnamein[128]= "ebnovalues.txt"; + char buf[128]; + + int nitems = 0; + int code_idx = 1; + int nthreads = 8; + int ch_type = CHANNEL_AWGN; + int ap_index = AP_NONE; + + // parse command line + while(--argc) { + argv++; + if (strncmp(*argv,"-h",2)==0) { + syntax(); + return 0; + } + else + if (strncmp(*argv,"-q",2)==0) { + code_idx = (int)atoi((*argv)+2); + if (code_idx>1) { + printf("Invalid code index\n"); + syntax(); + return -1; + } + } + else + if (strncmp(*argv,"-t",2)==0) { + nthreads = (int)atoi((*argv)+2); + + printf("nthreads = %d\n",nthreads); + + if (nthreads>NTHREADS_MAX) { + printf("Invalid number of threads\n"); + syntax(); + return -1; + } + } + else + if (strncmp(*argv,"-c",2)==0) { + ch_type = (int)atoi((*argv)+2); + if (ch_type>CHANNEL_RAYLEIGH) { + printf("Invalid channel type\n"); + syntax(); + return -1; + } + } + else + if (strncmp(*argv,"-a",2)==0) { + ap_index = (int)atoi((*argv)+2); + if (ap_index>AP_56) { + printf("Invalid a-priori information index\n"); + syntax(); + return -1; + } + } + else + if (strncmp(*argv,"-f",2)==0) { + strncpy(fnamein,(*argv)+2,127); + } + else + if (strncmp(*argv,"-h",2)==0) { + syntax(); + return -1; + } + else { + printf("Invalid option\n"); + syntax(); + return -1; + } + } -#ifdef TEST_WER_RAYLEIGH - wer_test_rayleigh(); -#endif + // parse points to be simulated from the input file + fin = fopen(fnamein,"r"); + if (!fin) { + printf("Can't open file: %s\n",fnamein); + syntax(); + } -#ifdef FTZ_ENABLE - // restore flush_zero_mode - _MM_SET_FLUSH_ZERO_MODE(ftzmode); -#endif + while (fgets(buf,128,fin)!=0) + if (*buf=='#' || *buf=='\n' ) + continue; + else + if (nitems==SIM_POINTS_MAX) + break; + else + if (sscanf(buf,"%f %u",&EbNodB[nitems],&nerrstgt[nitems])!=2) { + printf("Invalid input file format\n"); + syntax(); + return -1; + } + else + nitems++; + + fclose(fin); + + if (nitems==0) { + printf("No Eb/No point specified in file %s\n",fnamein); + syntax(); + return -1; + } + + printf("\nQ-ary Repeat-Accumulate Code Word Error Rate Simulator\n"); + printf("2016, Nico Palermo - IV3NWV\n\n"); + + printf("Nthreads = %d\n",nthreads); + printf("Channel = %s\n",ch_type==CHANNEL_AWGN?"AWGN":"Rayleigh"); + printf("Codename = %s\n",codetotest[code_idx]->name); + printf("A-priori = %s\n",ap_str[ap_index]); + printf("Eb/No input file = %s\n\n",fnamein); + + wer_test_proc(codetotest[code_idx], nthreads, ch_type, ap_index, EbNodB, nerrstgt, nitems); return 0; } diff --git a/qracodes/normrnd.c b/qracodes/normrnd.c index c952ce9..e2eac45 100644 --- a/qracodes/normrnd.c +++ b/qracodes/normrnd.c @@ -2,6 +2,10 @@ // functions to generate gaussian distributed numbers // // (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy +// +// Credits to Andrea Montefusco - IW0HDV for his help on adapting the sources +// to OSs other than MS Windows +// // ------------------------------------------------------------------------------ // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. @@ -16,7 +20,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . @@ -25,10 +29,9 @@ #if _WIN32 // note the underscore: without it, it's not msdn official! // Windows (x64 and x86) #include // required only for GetTickCount(...) - #define K_RAND_MAX UINT_MAX + #define K_RAND_MAX UINT_MAX #elif __unix__ // all unices, not all compilers #include - #include #define rand_s(x) (*x)=(unsigned int)lrand48() // returns unsigned integers in the range 0..0x7FFFFFFF #define K_RAND_MAX 0x7FFFFFFF // that's the max number generated by lrand48 // Unix @@ -41,11 +44,12 @@ // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though... #endif + // use MS rand_s(...) function void normrnd_s(float *dst, int nitems, float mean, float stdev) { unsigned int r; - float phi, u; + float phi=0, u=0; int set = 0; while (nitems--) @@ -54,18 +58,18 @@ void normrnd_s(float *dst, int nitems, float mean, float stdev) set = 0; } else { - rand_s(&r); phi = (M_2PI/(1.0f+K_RAND_MAX))*r; - rand_s(&r); u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+K_RAND_MAX))*(1.0f+r) ) ); + rand_s((unsigned int*)&r); phi = (M_2PI/(1.0f+K_RAND_MAX))*r; + rand_s((unsigned int*)&r); u = (float)sqrt(-2.0f* log( (1.0f/(1.0f+K_RAND_MAX))*(1.0f+r) ) ); *dst++ = (float)cos(phi)*u*stdev+mean; set=1; } } - +/* NOT USED // use MS rand() function void normrnd(float *dst, int nitems, float mean, float stdev) { - float phi, u; + float phi=0, u=0; int set = 0; while (nitems--) @@ -80,3 +84,4 @@ void normrnd(float *dst, int nitems, float mean, float stdev) set=1; } } +*/ diff --git a/qracodes/normrnd.h b/qracodes/normrnd.h index 3b59904..04bbf01 100644 --- a/qracodes/normrnd.h +++ b/qracodes/normrnd.h @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #ifndef _normrnd_h_ @@ -37,9 +37,11 @@ void normrnd_s(float *dst, int nitems, float mean, float stdev); // generate a random array of numbers with a gaussian distribution of given mean and stdev // use MS rand_s(...) function +/* not used void normrnd(float *dst, int nitems, float mean, float stdev); // generate a random array of numbers with a gaussian distribution of given mean and stdev // use MS rand() function +*/ #ifdef __cplusplus } diff --git a/qracodes/npfwht.c b/qracodes/npfwht.c index 13018ba..0f8fa42 100644 --- a/qracodes/npfwht.c +++ b/qracodes/npfwht.c @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #include "npfwht.h" diff --git a/qracodes/npfwht.h b/qracodes/npfwht.h index bddd8b0..212b0af 100644 --- a/qracodes/npfwht.h +++ b/qracodes/npfwht.h @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #ifndef _npfwht_h_ diff --git a/qracodes/pdmath.c b/qracodes/pdmath.c index 0ebad49..2c55e41 100644 --- a/qracodes/pdmath.c +++ b/qracodes/pdmath.c @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #include "pdmath.h" @@ -331,14 +331,14 @@ void pd_memset(float *dst, const float *src, int ndim, int nitems) } } -void pd_fwdperm(float *dst, float *src, const unsigned int *perm, int ndim) +void pd_fwdperm(float *dst, float *src, const int *perm, int ndim) { // TODO: non-loop implementation while (ndim--) dst[ndim] = src[perm[ndim]]; } -void pd_bwdperm(float *dst, float *src, const unsigned int *perm, int ndim) +void pd_bwdperm(float *dst, float *src, const int *perm, int ndim) { // TODO: non-loop implementation while (ndim--) diff --git a/qracodes/pdmath.h b/qracodes/pdmath.h index b082cf8..bd441e0 100644 --- a/qracodes/pdmath.h +++ b/qracodes/pdmath.h @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . @@ -57,12 +57,12 @@ float pd_norm(float *pd, int nlogdim); // are replaced with a uniform distribution // Returns the norm of the distribution prior to the normalization -void pd_fwdperm(float *dst, float *src, const unsigned int *perm, int ndim); +void pd_fwdperm(float *dst, float *src, const int *perm, int ndim); // Forward permutation of a distribution // Computes dst[k] = src[perm[k]] for every element in the distribution // perm must be a pointer to an array of integers of length ndim -void pd_bwdperm(float *dst, float *src, const unsigned int *perm, int ndim); +void pd_bwdperm(float *dst, float *src, const int *perm, int ndim); // Backward permutation of a distribution // Computes dst[perm[k]] = src[k] for every element in the distribution // perm must be a pointer to an array of integers of length ndim diff --git a/qracodes/qra12_63_64_irr_b.c b/qracodes/qra12_63_64_irr_b.c index 8e3f48a..d9bdf19 100644 --- a/qracodes/qra12_63_64_irr_b.c +++ b/qracodes/qra12_63_64_irr_b.c @@ -1,10 +1,10 @@ // qra12_63_64_irr_b.c -// Code tables and defines for the Q-ary irregular RA code (12,63) over GF(64) +// Encoding/Decoding tables for Q-ary RA code (12,63) over GF(64) // Code Name: qra12_63_64_irr_b -// Systematic symbols repetition factors: 333344455567 -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ +// (12,63) RA Code over GF(64) - RF=333344455567 + +// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy + // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. // @@ -16,17 +16,18 @@ // 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 qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #include "qra12_63_64_irr_b.h" -/* + #define qra_K 12 // number of information symbols #define qra_N 63 // codeword length in symbols #define qra_m 6 // bits/symbol #define qra_M 64 // Symbol alphabet cardinality +#define qra_a 1 // grouping factor #define qra_NC 51 // number of check symbols (N-K) // Defines used by the message passing decoder -------- @@ -34,14 +35,14 @@ #define qra_V 63 // number of variables in the code graph (N) #define qra_C 115 // number of factors in the code graph (N +(N-K)+1) #define qra_NMSG 217 // number of msgs in the code graph -#define qra_MAXVDEG 8 // maximum variable degree (intrinsic check included) +#define qra_MAXVDEG 8 // maximum variable degree #define qra_MAXCDEG 3 // maximum factor degree -#define qra_R 0.19048 // code rate (K/N) -*/ -// Tables used by the encoder ------------------------- +#define qra_R 0.19048f // code rate (K/N) +#define CODE_NAME "qra_12_63_64_irr_b" + // table of the systematic symbols indexes in the accumulator chain -const unsigned int qra_acc_input_idx[qra_NC+1] = { +static const int qra_acc_input_idx[qra_NC+1] = { 3, 11, 0, 1, 7, 8, 6, 5, 10, 4, 11, 9, 0, 2, 6, 7, 8, 4, 11, 5, 10, 2, 1, 9, 3, 8, 4, 11, 5, 7, @@ -51,7 +52,7 @@ const unsigned int qra_acc_input_idx[qra_NC+1] = { }; // table of the systematic symbols weight logarithms over GF(M) -const unsigned int qra_acc_input_wlog[qra_NC+1] = { +static const int qra_acc_input_wlog[qra_NC+1] = { 39, 0, 34, 16, 25, 0, 34, 48, 19, 13, 29, 56, 0, 5, 39, 42, 31, 0, 10, 0, 57, 62, 33, 43, 0, 14, 22, 48, 28, 20, @@ -61,7 +62,7 @@ const unsigned int qra_acc_input_wlog[qra_NC+1] = { }; // table of the logarithms of the elements of GF(M) (log(0) never used) -const unsigned int qra_log[qra_M] = { +static const int qra_log[qra_M] = { -1, 0, 1, 6, 2, 12, 7, 26, 3, 32, 13, 35, 8, 48, 27, 18, 4, 24, 33, 16, 14, 52, 36, 54, 9, 45, 49, 38, 28, 41, @@ -72,7 +73,7 @@ const unsigned int qra_log[qra_M] = { }; // table of GF(M) elements given their logarithm -const unsigned int qra_exp[qra_M-1] = { +static const int qra_exp[qra_M-1] = { 1, 2, 4, 8, 16, 32, 3, 6, 12, 24, 48, 35, 5, 10, 20, 40, 19, 38, 15, 30, 60, 59, 53, 41, 17, 34, 7, 14, 28, 56, @@ -82,10 +83,8 @@ const unsigned int qra_exp[qra_M-1] = { 57, 49, 33 }; -// Tables used by the decoder ------------------------- - // table of the messages weight logarithms over GF(M) -const unsigned int qra_msgw[qra_NMSG] = { +static const int qra_msgw[qra_NMSG] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, @@ -111,7 +110,7 @@ const unsigned int qra_msgw[qra_NMSG] = { }; // table of the degrees of the variable nodes -const unsigned int qra_vdeg[qra_V] = { +static const int qra_vdeg[qra_V] = { 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, @@ -122,7 +121,7 @@ const unsigned int qra_vdeg[qra_V] = { }; // table of the degrees of the factor nodes -const unsigned int qra_cdeg[qra_C] = { +static const int qra_cdeg[qra_C] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, @@ -138,7 +137,7 @@ const unsigned int qra_cdeg[qra_C] = { }; // table (uncompressed) of the v->c message indexes (-1=unused entry) -const unsigned int qra_v2cmidx[qra_V*qra_MAXVDEG] = { +static const int qra_v2cmidx[qra_V*qra_MAXVDEG] = { 0, 65, 75, 101, -1, -1, -1, -1, 1, 66, 85, 110, -1, -1, -1, -1, 2, 76, 84, 106, -1, -1, -1, -1, @@ -205,7 +204,7 @@ const unsigned int qra_v2cmidx[qra_V*qra_MAXVDEG] = { }; // table (uncompressed) of the c->v message indexes (-1=unused entry) -const unsigned int qra_c2vmidx[qra_C*qra_MAXCDEG] = { +static const int qra_c2vmidx[qra_C*qra_MAXCDEG] = { 0, -1, -1, 1, -1, -1, 2, -1, -1, 3, -1, -1, 4, -1, -1, 5, -1, -1, 6, -1, -1, 7, -1, -1, 8, -1, -1, 9, -1, -1, 10, -1, -1, 11, -1, -1, @@ -238,7 +237,7 @@ const unsigned int qra_c2vmidx[qra_C*qra_MAXCDEG] = { }; // permutation matrix to compute Prob(x*alfa^logw) -const unsigned int qra_pmat[qra_M*qra_M] = { +static const int qra_pmat[qra_M*qra_M] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, @@ -493,5 +492,43 @@ const unsigned int qra_pmat[qra_M*qra_M] = { 35, 33, 39, 37, 43, 41, 47, 45, 51, 49, 55, 53, 59, 57, 63, 61 }; +const qracode qra_12_63_64_irr_b = { + qra_K, + qra_N, + qra_m, + qra_M, + qra_a, + qra_NC, + qra_V, + qra_C, + qra_NMSG, + qra_MAXVDEG, + qra_MAXCDEG, + QRATYPE_NORMAL, + qra_R, + CODE_NAME, + qra_acc_input_idx, + qra_acc_input_wlog, + qra_log, + qra_exp, + qra_msgw, + qra_vdeg, + qra_cdeg, + qra_v2cmidx, + qra_c2vmidx, + qra_pmat +}; - +#undef qra_K +#undef qra_N +#undef qra_m +#undef qra_M +#undef qra_a +#undef qra_NC +#undef qra_V +#undef qra_C +#undef qra_NMSG +#undef qra_MAXVDEG +#undef qra_MAXCDEG +#undef qra_R +#undef CODE_NAME diff --git a/qracodes/qra12_63_64_irr_b.h b/qracodes/qra12_63_64_irr_b.h index b8edaee..3641e4d 100644 --- a/qracodes/qra12_63_64_irr_b.h +++ b/qracodes/qra12_63_64_irr_b.h @@ -1,10 +1,10 @@ // qra12_63_64_irr_b.h -// Code tables and defines for the Q-ary irregular RA code (12,63) over GF(64) +// Code tables and defines for Q-ary RA code (12,63) over GF(64) // Code Name: qra12_63_64_irr_b -// Systematic symbols repetition factors: 333344455567 -// -// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy -// ------------------------------------------------------------------------------ +// (12,63) RA Code over GF(64) - RF=333344455567 + +// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy + // This file is part of the qracodes project, a Forward Error Control // encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. // @@ -16,49 +16,24 @@ // 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 qra_codes source distribution. +// along with qracodes source distribution. // If not, see . -#ifndef _qra12_64_64_irr_b_h -#define _qra12_64_64_irr_b_h - - -#define qra_K 12 // number of information symbols -#define qra_N 63 // codeword length in symbols -#define qra_m 6 // bits/symbol -#define qra_M 64 // Symbol alphabet cardinality -#define qra_NC 51 // number of check symbols (N-K) - -// Defines used by the message passing decoder -------- +#ifndef _qra12_63_64_irr_b_h +#define _qra12_63_64_irr_b_h -#define qra_V 63 // number of variables in the code graph (N) -#define qra_C 115 // number of factors in the code graph (N +(N-K)+1) -#define qra_NMSG 217 // number of msgs in the code graph -#define qra_MAXVDEG 8 // maximum variable degree (intrinsic check included) -#define qra_MAXCDEG 3 // maximum factor degree -#define qra_R 0.19048 // code rate (K/N) +#include "qracodes.h" #ifdef __cplusplus extern "C" { #endif -// Tables used by the encoder ------------------------- -extern const unsigned int qra_acc_input_idx[qra_NC+1]; -extern const unsigned int qra_acc_input_wlog[qra_NC+1]; -extern const unsigned int qra_log[qra_M]; -extern const unsigned int qra_exp[qra_M-1]; -// Tables used by the decoder ------------------------- -extern const unsigned int qra_msgw[qra_NMSG]; -extern const unsigned int qra_vdeg[qra_V]; -extern const unsigned int qra_cdeg[qra_C]; -extern const unsigned int qra_v2cmidx[qra_V*qra_MAXVDEG]; -extern const unsigned int qra_c2vmidx[qra_C*qra_MAXCDEG]; -extern const unsigned int qra_pmat[qra_M*qra_M]; +extern const qracode qra_12_63_64_irr_b; #ifdef __cplusplus } #endif -#endif // _qra12_64_64_irr_b_h +#endif // _qra12_63_64_irr_b_h diff --git a/qracodes-mt/qra13_64_64_irr_e.c b/qracodes/qra13_64_64_irr_e.c similarity index 100% rename from qracodes-mt/qra13_64_64_irr_e.c rename to qracodes/qra13_64_64_irr_e.c diff --git a/qracodes-mt/qra13_64_64_irr_e.h b/qracodes/qra13_64_64_irr_e.h similarity index 100% rename from qracodes-mt/qra13_64_64_irr_e.h rename to qracodes/qra13_64_64_irr_e.h diff --git a/qracodes/qracodes.c b/qracodes/qracodes.c index b04bfb1..4b9fe34 100644 --- a/qracodes/qracodes.c +++ b/qracodes/qracodes.c @@ -16,7 +16,7 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #include @@ -26,55 +26,97 @@ #include "pdmath.h" #include "qracodes.h" -#include "qra12_63_64_irr_b.h" // tables for the P(12,63) Q-ary RA code over GF(64) -// #define QRA_DEBUG - -// message tables used by the qra_extrinsic -// table of the v->c messages -static float qra_v2cmsg[qra_NMSG*qra_M]; -// table of the c->v messages -static float qra_c2vmsg[qra_NMSG*qra_M]; - -int qra_encode(uint *y, const uint *x) +int qra_encode(const qracode *pcode, int *y, const int *x) { - int k; - uint t, chk = 0; + int k,j,kk,jj; + int t, chk = 0; + + const int K = pcode->K; + const int M = pcode->M; + const int NC= pcode->NC; + const int a = pcode->a; + const int *acc_input_idx = pcode->acc_input_idx; + const int *acc_input_wlog = pcode->acc_input_wlog; + const int *gflog = pcode->gflog; + const int *gfexp = pcode->gfexp; // copy the systematic symbols to destination - memcpy(y,x,qra_K*sizeof(uint)); + memcpy(y,x,K*sizeof(int)); - y = y+qra_K; // point to check symbols + y = y+K; // point to check symbols // compute the code check symbols as a weighted accumulation of a permutated // sequence of the (repeated) systematic input symbols: // chk(k+1) = x(idx(k))*alfa^(logw(k)) + chk(k) - // (all operations performed over GF(qra_M)) - for (k=0;k 1 + for (k=0;kM; + const int qra_m = pcode->m; + const int qra_V = pcode->V; + const int qra_MAXVDEG = pcode->MAXVDEG; + const int *qra_vdeg = pcode->vdeg; + const int qra_C = pcode->C; + const int qra_MAXCDEG = pcode->MAXCDEG; + const int *qra_cdeg = pcode->cdeg; + const int *qra_v2cmidx = pcode->v2cmidx; + const int *qra_c2vmidx = pcode->c2vmidx; + const int *qra_pmat = pcode->gfpmat; + const int *qra_msgw = pcode->msgw; + +// float msgout[qra_M]; // buffer to store temporary results + float msgout[QRACODE_MAX_M]; // we use a fixed size in order to avoid mallocs + + float totex; // total extrinsic information int nit; // current iteration - int nv; // current variable - int nc; // current check + int nv; // current variable + int nc; // current check int k,kk; // loop indexes int ndeg; // current node degree @@ -187,8 +254,12 @@ int qra_extrinsic(float *pex, const float *pix, int maxiter) int rc = -1; // rc>=0 extrinsic converged to 1 at iteration rc (rc=0..maxiter-1) // rc=-1 no convergence in the given number of iterations // rc=-2 error in the code tables (code checks degrees must be >1) + // rc=-3 M is larger than QRACODE_MAX_M + + if (qra_M>QRACODE_MAX_M) + return -3; // message initialization ------------------------------------------------------- @@ -375,7 +446,7 @@ int qra_extrinsic(float *pex, const float *pix, int maxiter) return rc; } -void qra_mapdecode(uint *xdec, float *pex, const float *pix) +void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix) { // Maximum a posteriori probability decoding. // Given the intrinsic information (pix) and extrinsic information (pex) (computed with qra_extrinsic(...)) @@ -385,7 +456,11 @@ void qra_mapdecode(uint *xdec, float *pex, const float *pix) // Returns: // xdec[k] = decoded (information) symbols k=[0..qra_K-1] -// Note: pex is destroyed and overwritten with pmap +// Note: pex is destroyed and overwritten with mapp + + const int qra_M = pcode->M; + const int qra_m = pcode->m; + const int qra_K = pcode->K; int k; diff --git a/qracodes/qracodes.h b/qracodes/qracodes.h index f87d0b0..fb74b2e 100644 --- a/qracodes/qracodes.h +++ b/qracodes/qracodes.h @@ -16,22 +16,61 @@ // GNU General Public License for more details. // You should have received a copy of the GNU General Public License -// along with qra_codes source distribution. +// along with qracodes source distribution. // If not, see . #ifndef _qracodes_h_ #define _qracodes_h_ -typedef unsigned int uint; +// type of codes +#define QRATYPE_NORMAL 0x00 // normal code +#define QRATYPE_CRC 0x01 // code with crc - last information symbol is a CRC +#define QRATYPE_CRCPUNCTURED 0x02 // the CRC symbol is punctured (not sent along the channel) + + +typedef struct { + // code parameters + const int K; // number of information symbols + const int N; // codeword length in symbols + const int m; // bits/symbol + const int M; // Symbol alphabet cardinality (2^m) + const int a; // code grouping factor + const int NC; // number of check symbols (N-K) + const int V; // number of variables in the code graph (N) + const int C; // number of factors in the code graph (N +(N-K)+1) + const int NMSG; // number of msgs in the code graph + const int MAXVDEG; // maximum variable degree + const int MAXCDEG; // maximum factor degree + const int type; // see QRATYPE_xx defines + const float R; // code rate (K/N) + const char name[64]; // code name + // tables used by the encoder + const int *acc_input_idx; + const int *acc_input_wlog; + const int *gflog; + const int *gfexp; + // tables used by the decoder ------------------------- + const int *msgw; + const int *vdeg; + const int *cdeg; + const int *v2cmidx; + const int *c2vmidx; + const int *gfpmat; +} qracode; +// Uncomment the header file of the code which needs to be tested + +//#include "qra12_63_64_irr_b.h" // irregular code (12,63) over GF(64) +//#include "qra13_64_64_irr_e.h" // irregular code with good performance and best UER protection at AP56 +//#include "qra13_64_64_reg_a.h" // regular code with good UER but perf. inferior to that of code qra12_63_64_irr_b #ifdef __cplusplus extern "C" { #endif -int qra_encode(uint *y, const uint *x); -void qra_mfskbesselmetric(float *pix, float EsNoMetric); -int qra_extrinsic(float *pex, const float *pix, int maxiter); -void qra_mapdecode(uint *xdec, float *pex, const float *pix); +int qra_encode(const qracode *pcode, int *y, const int *x); +void qra_mfskbesselmetric(float *pix, const float *rsq, const int m, const int N, float EsNoMetric); +int qra_extrinsic(const qracode *pcode, float *pex, const float *pix, int maxiter,float *qra_v2cmsg,float *qra_c2vmsg); +void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix); #ifdef __cplusplus } diff --git a/qracodes/qracodes.vcproj b/qracodes/qracodes.vcproj index 445f06e..991be6e 100644 --- a/qracodes/qracodes.vcproj +++ b/qracodes/qracodes.vcproj @@ -3,9 +3,8 @@ ProjectType="Visual C++" Version="9,00" Name="qracodes" - ProjectGUID="{8540B348-B8F2-4104-84CB-6286AC19158B}" + ProjectGUID="{933A58F6-199B-4723-ACFE-3013E6DD9D0A}" RootNamespace="qracodes" - Keyword="Win32Proj" TargetFrameworkVersion="196613" > @@ -41,13 +40,9 @@ @@ -62,9 +57,7 @@ /> + + @@ -220,6 +222,10 @@ RelativePath=".\qra12_63_64_irr_b.h" > + +