WSJT-X/lib/qra/qra65/main.c

// main.c 
// QRA65 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
//    qra65.c/.h     - qra65 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 "     "
//    ../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 <http://www.gnu.org/licenses/>.

// -----------------------------------------------------------------------------

// The code used by the QRA65 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 <windows.h>   // required only for GetTickCount(...)
	#include <process.h>   // _beginthread
#endif

#if __linux__
#include <unistd.h>
#include <time.h>

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 <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "qra65.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<size;k++)
		printf("%2d ",x[k]);
	printf("\n");
}
void printwordh(char *msg, int *x, int size)
{
	int k;
	printf("\n%s ",msg);
	for (k=0;k<size;k++)
		printf("%02hx ",x[k]);
	printf("\n");
}

#define NSAMPLES (QRA65_N*QRA65_M)

static float rp[NSAMPLES];
static float rq[NSAMPLES];
static float chp[NSAMPLES];
static float chq[NSAMPLES];
static float r[NSAMPLES];

float *mfskchannel(int *x, int channel_type, float EbNodB)
{
	// simulate a MFSK channel (AWGN or Rayleigh)
	// x is a pointer to the transmitted codeword (an array of QRA65_N integers in the range 0..63)
	//
	// the function returns the received symbol energies (squared amplitudes)
	// as an array of (QRA65_M*QRA65_N) float numbers
	// the first QRA65_M entries of this array are the energies of the first symbol in the codeword
	// the second QRA65_M entries are those of the second symbol
	// and so on up to the last codeword symbol

	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
	const float R = 1.0f*QRA65_K/QRA65_N;	

	float EbNo = (float)pow(10,EbNodB/10);
	float EsNo = 1.0f*QRA65_m*R*EbNo;
	float Es = EsNo*No;
	float A = (float)sqrt(Es);
	int k;

	normrnd_s(rp,NSAMPLES,0,sigma);
	normrnd_s(rq,NSAMPLES,0,sigma);

	if (channel_type == CHANNEL_AWGN) 
			for (k=0;k<QRA65_N;k++) 
				rp[k*QRA65_M+x[k]]+=A;
	else 
		if (channel_type == CHANNEL_RAYLEIGH) {
			normrnd_s(chp,QRA65_N,0,sigmach);
			normrnd_s(chq,QRA65_N,0,sigmach);
			for (k=0;k<QRA65_N;k++) {
				rp[k*QRA65_M+x[k]]+=A*chp[k];
				rq[k*QRA65_M+x[k]]+=A*chq[k];
				}
			}
		else {
			return 0;	// unknown channel type
			}

	// compute the squares of the amplitudes of the received samples
	for (k=0;k<NSAMPLES;k++) 
		r[k] = rp[k]*rp[k] + rq[k]*rq[k];

	return r;

}

// These defines are some packed fields as computed by JT65 
#define CALL_IV3NWV		0x7F85AE7	
#define CALL_K1JT		0xF70DDD7
#define GRID_JN66		0x3AE4		// JN66
#define GRID_73 		0x7ED0		// 73

char decode_type[6][32] = {
	"[?    ?    ?] AP0",
	"[CQ   ?    ?] AP27",
	"[CQ   ?     ] AP44",
	"[CALL ?    ?] AP29",
	"[CALL ?     ] AP45",
	"[CALL CALL ?] AP57"
};

int test_proc_1(int channel_type, float EbNodB, int mode)
{
	// Here we simulate the following (dummy) QSO:
	//
	// 1) CQ IV3NWV
	// 2)                 IV3NWV K1JT
	// 3) K1JT IV3NWV 73
	// 4)                 IV3NWV K1JT 73
	//
	// No message repetition is attempted
	// The QSO is counted as successfull if IV3NWV received the last message
	// When mode=QRA_AUTOAP each decoder attempts to decode the message sent
	// by the other station using the a-priori information derived by what
	// has been already decoded in a previous phase of the QSO if
	// decoding with no a-priori information has not been successful.
	// I.e. at step 1) K1JT's decoder first attempts to decode msgs of type [? ? ?] and if this
	// attempt fails, it attempts to decode [CQ/QRZ ? ?] or [[CQ/QRZ ?] msgs
	// At step 2) if IV3NWV's decoder is unable to decode K1JT's without AP it attempts to decode
	// messages of the type [IV3NWV ? ?] and [IV3NWV ?].
	// At step 3) K1JT's decoder attempts to decode [? ? ?] and [K1JT IV3NWV ?] (this last decode
	// type has been enabled by K1JT's encoder at step 2)
	// At step 4) IV3NWV's decoder attempts to decode [? ? ?] and [IV3NWV K1JT ?] (this last decode
	// type has been enabled by IV3NWV's encoder at step 3)

	// At each step the simulation reports if a decode was successful.
	// In this case it also reports the type of decode (see table decode_type above)

	// When mode=QRA_NOAP, only [? ? ?] decodes are attempted and no a-priori information
	// is used by the decoder

	// The function returns 0 if all of the four messages have been decoded by their recipients
	// (with no retries) and -1 if any of them could not be decoded


	int x[QRA65_K], xdec[QRA65_K];
	int y[QRA65_N];
	float *rx;
	int rc;

	// each simulated station must use its own codec
	// as it might work with different a-priori information
	qra65codec *codec_iv3nwv = qra65_init(mode,CALL_IV3NWV);	// codec for IV3NWV
	qra65codec *codec_k1jt   = qra65_init(mode,CALL_K1JT);		// codec for K1JT


	// IV3NWV makes a CQ call (with no grid)
	printf("IV3NWV tx: CQ IV3NWV\n");
	encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_BLANK);
	qra65_encode(codec_iv3nwv, y, x);
	rx = mfskchannel(y,channel_type,EbNodB);

	// K1JT attempts to decode
	rc = qra65_decode(codec_k1jt, xdec,rx);
	if (rc>=0) { // decoded
		printf("K1JT   rx: received with apcode=%d %s\n",rc, decode_type[rc]);
		// K1JT replies to IV3NWV (with no grid)
		printf("K1JT   tx: IV3NWV K1JT\n");
		encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_BLANK);
		qra65_encode(codec_k1jt, y, x);
		rx = mfskchannel(y,channel_type,EbNodB);

		// IV3NWV attempts to decode
		rc = qra65_decode(codec_iv3nwv, xdec,rx);
		if (rc>=0) { // decoded
			printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]);
			// IV3NWV replies to K1JT with a 73
			printf("IV3NWV tx: K1JT   IV3NWV 73\n");
			encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV, GRID_73);
			qra65_encode(codec_iv3nwv, y, x);
			rx = mfskchannel(y,channel_type,EbNodB);

			// K1JT attempts to decode
			rc = qra65_decode(codec_k1jt, xdec,rx);
			if (rc>=0) { // decoded
				printf("K1JT   rx: received with apcode=%d %s\n",rc, decode_type[rc]);
				// K1JT replies to IV3NWV with a 73
				printf("K1JT   tx: IV3NWV K1JT   73\n");
				encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_73);
				qra65_encode(codec_k1jt, y, x);
				rx = mfskchannel(y,channel_type,EbNodB);

				// IV3NWV attempts to decode
				rc = qra65_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 qra65_decode function
	// indicates the amount of a-priori information required to decode the received message
	// accordingly to this table:
	//	rc=0    [?    ?    ?] AP0
	//  rc=1    [CQ   ?    ?] AP27
	//  rc=2    [CQ   ?     ] AP44
	//  rc=3    [CALL ?    ?] AP29
	//  rc=4    [CALL ?     ] AP45
	//  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[QRA65_K], xdec[QRA65_K];
	int y[QRA65_N];
	float *rx;
	int rc,k;

	int ndecok[6] = { 0, 0, 0, 0, 0, 0};
	int ntx = 100,ndec=0;

	qra65codec *codec_iv3nwv = qra65_init(mode,CALL_IV3NWV);	// codec for IV3NWV
	qra65codec *codec_k1jt   = qra65_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);
	qra65_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);
	qra65_encode(codec_iv3nwv, y, x);

	printf("Simulating decodes by K1JT up to AP56 ...");

	for (k=0;k<ntx;k++) {
		printf(".");
		rx = mfskchannel(y,channel_type,EbNodB);
		rc = qra65_decode(codec_k1jt, xdec,rx);
		if (rc>=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("\nQRA65 Mode Tests\n");
	printf("2016, Nico Palermo - IV3NWV\n\n");
	printf("---------------------------\n\n");
	printf("Syntax: qra65 [-s<snrdb>] [-c<channel>] [-a<ap-type>] [-t<testtype>] [-h]\n");
	printf("Options: \n");
	printf("       -s<snrdb>   : set simulation SNR in 2500 Hz BW (default:-27.5 dB)\n");
	printf("       -c<channel> : set channel type 0=AWGN (default) 1=Rayleigh\n");
	printf("       -a<ap-type> : set decode type 0=NO_AP 1=AUTO_AP (default)\n");
	printf("       -t<testtype>: 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 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<nqso;k++) {
			printf("\n\n------------------------\n");
			rc = test_proc_1(channel, EbNodB, mode);
			if (rc==0)
				nok++;
			}
		printf("\n\n%d/%d QSOs to end without repetitions\n",nok,nqso);
		}
	else
		test_proc_2(channel, EbNodB, mode);

	printf("SNR = %.1fdB channel=%s ap-mode=%s\n\n",
			SNRdB,
			channel==CHANNEL_AWGN?"AWGN":"RAYLEIGH",
			mode==QRA_NOAP?"NO_AP":"AUTO_AP"
			);


	return 0;
}