code:
/*
qra64.c
Encoding/decoding functions for the QRA64 mode

(c) 2016 - Nico Palermo, IV3NWV

-------------------------------------------------------------------------------

   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/>.

-----------------------------------------------------------------------------

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

#include "qra64.h"
#include "../qracodes/qracodes.h"
#include "../qracodes/qra13_64_64_irr_e.h"
#include "../qracodes/pdmath.h"
#include "../qracodes/normrnd.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_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix);
static int qra64_do_decode(int *x, const float *pix, const int *ap_mask,
			    const int *ap_x);
static float qra64_fastfading_estim_noise_std(
				const float *rxen,
				const float esnometric,
				const int submode);

static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxen,
				const float *hptr,
				const int    hlen,
				const float sigma,
				const float EsNoMetric,
				const int submode);

static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxen,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode);

// a-priori information masks for fields in JT65-like msgs --------------------

// when defined limits the AP masks to reduce the false decode rate
#define LIMIT_AP_MASKS

#ifdef LIMIT_AP_MASKS
#define MASK_CQQRZ      0xFFFFFFC
#define MASK_CALL1      0xFFFFFFC
#define MASK_CALL2      0xFFFFFFC
#define MASK_GRIDFULL	0x3FFC
#define MASK_GRIDFULL12	0x3FFC
#define MASK_GRIDBIT	0x8000
#else
#define MASK_CQQRZ      0xFFFFFFC
#define MASK_CALL1      0xFFFFFFF
#define MASK_CALL2      0xFFFFFFF
#define MASK_GRIDFULL	0xFFFF
#define MASK_GRIDFULL12	0x3FFC
#define MASK_GRIDBIT	0x8000	  // b[15] is 1 for free text, 0 otherwise
#endif

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

qra64codec *qra64_init(int flags)
{

  // 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;

  memset(pcodec->apmsg_set,0,APTYPE_SIZE*sizeof(int));

  if (flags==QRA_NOAP)
    return pcodec;

  // for QRA_USERAP and QRA_AUTOAP modes we always enable [CQ/QRZ ? ?] mgs look-up.
  // encode CQ/QRZ AP messages
  // 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
  qra64_apset(pcodec, CALL_CQ, 0, GRID_BLANK, APTYPE_CQQRZ);

  // initialize masks for decoding with a-priori information
  encodemsg_jt65(pcodec->apmask_cqqrz,     MASK_CQQRZ, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cqqrz_ooo, MASK_CQQRZ, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1,     MASK_CALL1, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_ooo, MASK_CALL1, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call2,     0, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call2_ooo, 0, MASK_CALL2, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1_call2,     MASK_CALL1,MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_call2_grid,MASK_CALL1,MASK_CALL2, MASK_GRIDFULL12);
  encodemsg_jt65(pcodec->apmask_cq_call2,     MASK_CQQRZ, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cq_call2_ooo, MASK_CQQRZ, MASK_CALL2, MASK_GRIDFULL12);

  return pcodec;
}

void qra64_close(qra64codec *pcodec)
{
	free(pcodec);
}

int qra64_apset(qra64codec *pcodec, const int mycall, const int hiscall, const int grid, const int aptype)
{
// Set decoder a-priori knowledge accordingly to the type of the message to look up for
// arguments:
//		pcodec    = pointer to a qra64codec data structure as returned by qra64_init
//		mycall    = mycall to look for
//		hiscall   = hiscall to look for
//		grid      = grid to look for
//		aptype    = define and masks the type of AP to be set accordingly to the following:
//			APTYPE_CQQRZ     set [cq/qrz ?       ?/blank]
//			APTYPE_MYCALL    set [mycall ?       ?/blank]
//			APTYPE_HISCALL   set [?      hiscall ?/blank]
//			APTYPE_BOTHCALLS set [mycall hiscall ?]
//			APTYPE_FULL		 set [mycall hiscall grid]
//			APTYPE_CQHISCALL set [cq/qrz hiscall ?/blank] and [cq/qrz hiscall grid]
// returns:
//		0   on success
//      -1  when qra64_init was called with the QRA_NOAP flag
//		-2  invalid apytpe

	if (pcodec->apflags==QRA_NOAP)
		return -1;

	switch (aptype) {
		case APTYPE_CQQRZ:
			encodemsg_jt65(pcodec->apmsg_cqqrz,  CALL_CQ, 0, GRID_BLANK);
			break;
		case APTYPE_MYCALL:
			encodemsg_jt65(pcodec->apmsg_call1,  mycall,  0, GRID_BLANK);
			break;
		case APTYPE_HISCALL:
			encodemsg_jt65(pcodec->apmsg_call2,  0, hiscall, GRID_BLANK);
			break;
		case APTYPE_BOTHCALLS:
			encodemsg_jt65(pcodec->apmsg_call1_call2,  mycall, hiscall, GRID_BLANK);
			break;
		case APTYPE_FULL:
			encodemsg_jt65(pcodec->apmsg_call1_call2_grid,  mycall, hiscall, grid);
			break;
		case APTYPE_CQHISCALL:
			encodemsg_jt65(pcodec->apmsg_cq_call2,      CALL_CQ, hiscall, GRID_BLANK);
			encodemsg_jt65(pcodec->apmsg_cq_call2_grid, CALL_CQ, hiscall, grid);
			break;
		default:
			return -2;	// invalid ap type
		}

	  pcodec->apmsg_set[aptype]=1;	// signal the decoder to look-up for the specified type

	  return 0;
}
void qra64_apdisable(qra64codec *pcodec, const int aptype)
{
	if (pcodec->apflags==QRA_NOAP)
		return;

	if (aptype<APTYPE_CQQRZ || aptype>=APTYPE_SIZE)
		return;

	pcodec->apmsg_set[aptype] = 0;	//  signal the decoder not to look-up to the specified type
}

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 hiscall,mycall,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;

  // Here we handle the QRA_AUTOAP mode --------------------------------------------

  // When a [hiscall mycall ?] msg is detected we instruct the decoder
  // to look for [mycall hiscall ?] msgs
  // otherwise when a [cq mycall ?] msg is sent we reset the APTYPE_BOTHCALLS

  // 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, nothing to do

  // It's a [hiscall mycall grid] message

  // We assume that mycall is our call (but we don't check it)
  // hiscall the station we are calling or a general call (CQ/QRZ/etc..)
  decodemsg_jt65(&hiscall,&mycall,&grid,x);

  if ((hiscall>=CALL_CQ && hiscall<=CALL_CQ999) || hiscall==CALL_CQDX ||
      hiscall==CALL_DE) {
	// tell the decoder to look for msgs directed to us
	qra64_apset(pcodec,mycall,0,0,APTYPE_MYCALL);
    // We are making a general call and don't know who might reply
    // Reset APTYPE_BOTHCALLS so decoder won't look for [mycall hiscall ?] msgs
    qra64_apdisable(pcodec,APTYPE_BOTHCALLS);
  } else {
    // We are replying to someone named hiscall
    // Set APTYPE_BOTHCALLS so decoder will try for [mycall hiscall ?] msgs
    qra64_apset(pcodec,mycall, hiscall, GRID_BLANK, APTYPE_BOTHCALLS);
  }

}

#define EBNO_MIN -10.0f		// minimum Eb/No value returned by the decoder (in dB)
// AWGN metric decoder
int qra64_decode(qra64codec *pcodec, float *ebno, int *x, const float *rxen)
{
  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise variance
  float msge;						// estimated message energy
  float ebnoval;					// estimated Eb/No
  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
  noisestd = 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<QRA64_C;k++) {
    pd_init(dsttmp,srctmp,QRA64_M);
    dsttmp -=QRA64_M;
    srctmp -=QRA64_M;
  }
  // Initialize crc prob to a uniform distribution
  pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

  // Try to decode using all AP cases if required
  rc = qra64_decode_attempts(pcodec, xdec, ix);

  if (rc<0)
	  return rc;	// no success

  // successfull decode --------------------------------

  // copy decoded message (without crc) to output buffer
  memcpy(x,xdec,QRA64_K*sizeof(int));

  if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
	  return rc;

  // reencode message and estimate Eb/No
  qra_encode(&QRA64_CODE, ydec, xdec);
  // puncture crc
  memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));
  // compute total power of decoded message
  msge = 0;
  for (k=0;k<QRA64_N;k++) {
	  msge +=rxen[ydec[k]];	// add energy of current symbol
	  rxen+=QRA64_M;			// ptr to next symbol
	  }

  // NOTE:
  // To make a more accurate Eb/No estimation we should compute the noise variance
  // on all the rxen values but the transmitted symbols.
  // Noisestd is compute by qra_mfskbesselmetric assuming that
  // the signal power is much less than the total noise power in the QRA64_M tones
  // but this is true only if the Eb/No is low.
  // Here, in order to improve accuracy, we linearize the estimated Eb/No value empirically
  // (it gets compressed when it is very high as in this case the noise variance
  // is overestimated)

  // this would be the exact value if the noisestd were not overestimated at high Eb/No
  ebnoval = (0.5f/(QRA64_K*QRA64_m))*msge/(noisestd*noisestd)-1.0f;

  // Empirical linearization (to remove the noise variance overestimation)
  // the resulting SNR is accurate up to +20 dB (51 dB Eb/No)
  if (ebnoval>57.004f)
	  ebnoval=57.004f;
  ebnoval = ebnoval*57.03f/(57.03f-ebnoval);

  // compute value in dB
  if (ebnoval<=0) {
    ebnoval = EBNO_MIN; // assume a minimum, positive value
  }
  else {
    ebnoval = 10.0f*(float)log10(ebnoval);
    if (ebnoval<EBNO_MIN) {
      ebnoval = EBNO_MIN;
    }
  }

  *ebno = ebnoval;

  return rc;
}

//
// Fast-fading / Rayleigh channel metric decoder ----------------------------------------------
//
// Tables of fading energies coefficients for QRA64 (Ts=6912/12000)
// As the fading is assumed to be symmetric around the nominal frequency
// only the leftmost and the central coefficient are stored in the tables.
// (files have been generated with the Matlab code efgengaussenergy.m and efgenlorentzenergy.m)
#include "fadengauss.c"
#include "fadenlorentz.c"

int qra64_decode_fastfading(
				qra64codec *pcodec,		// ptr to the codec structure
				float *ebno,			// ptr to where the estimated Eb/No value will be saved
				int *x,					// ptr to decoded message
				const float *rxen,		// ptr to received symbol energies array
				const int submode,		// submode idx (0=QRA64A ... 4=QRA64E)
				const float B90,	    // spread bandwidth (90% fractional energy)
				const int fadingModel)  // 0=Gaussian 1=Lorentzian fade model

// Decode a QRA64 msg using a fast-fading metric
//
// rxen: The array of the received bin energies
//       Bins must be spaced by integer multiples of the symbol rate (1/Ts Hz)
//       The array must be an array of total length U = L x N where:
//			L: is the number of frequency bins per message symbol (see after)
//          N: is the number of symbols in a QRA64 msg (63)

//       The number of bins/symbol L depends on the selected submode accordingly to
//		 the following rule:
//			L = (64+64*2^submode+64) = 64*(2+2^submode)
//		 Tone 0 is always supposed to be at offset 64 in the array.
//		 The m-th tone nominal frequency is located at offset 64 + m*2^submode (m=0..63)
//
//		 Submode A: (2^submode = 1)
//          L = 64*3 = 196 bins/symbol
//          Total length of the energies array: U = 192*63 = 12096 floats
//
//		 Submode B: (2^submode = 2)
//          L = 64*4 = 256 bins/symbol
//          Total length of the energies array: U = 256*63 = 16128 floats
//
//		 Submode C: (2^submode = 4)
//          L = 64*6 = 384 bins/symbol
//          Total length of the energies array: U = 384*63 = 24192 floats
//
//		 Submode D: (2^submode = 8)
//          L = 64*10 = 640 bins/symbol
//          Total length of the energies array: U = 640*63 = 40320 floats
//
//		 Submode E: (2^submode = 16)
//          L = 64*18 = 1152 bins/symbol
//          Total length of the energies array: U = 1152*63 = 72576 floats
//
//		Note: The rxen array is modified and reused for internal calculations.
//
//
//	B90: spread fading bandwidth in Hz (90% fractional average energy)
//
//			B90 should be in the range 1 Hz ... 238 Hz
//			The value passed to the call is rounded to the closest value among the
//			64 available values:
//				B = 1.09^k Hz, with k=0,1,...,63
//
//			I.e. B90=27 Hz will be approximated in this way:
//				k = rnd(log(27)/log(1.09)) = 38
//              B90 = 1.09^k = 1.09^38 = 26.4 Hz
//
//          For any input value the maximum rounding error is not larger than +/- 5%
//

{

  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise std
  float esno,ebnoval;				// estimated Eb/No
  float tempf;
  float EsNoMetric;

  int rc;
  int hidx, hlen;
  const float *hptr;

	if (QRA64_NMSG!=QRA64_CODE.NMSG)
		return -16;					// QRA64_NMSG define is wrong

	if (submode<0 || submode>4)
		return -17;				// invalid submode

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate amplitude weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian energy weighting taps
		hlen = glen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian energy weighting taps
		hlen = glen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// compute (euristically) the optimal decoder metric accordingly the given spread amount
	// We assume that the decoder threshold is:
	//		Es/No(dB) = Es/No(AWGN)(dB) + 8*log(B90)/log(240)(dB)
	// that's to say, at the maximum Doppler spread bandwidth (240 Hz) there's a ~8 dB Es/No degradation
	// over the AWGN case
	tempf = 8.0f*(float)log((float)B90)/(float)log(240.0f);
	EsNoMetric = pcodec->decEsNoMetric*(float)pow(10.0f,tempf/10.0f);

	// Step 1 -----------------------------------------------------------------------------------
	// Evaluate the noise stdev from the received energies at nominal tone frequencies
	noisestd = qra64_fastfading_estim_noise_std(rxen, EsNoMetric, submode);

	// Step 2 -----------------------------------------------------------------------------------
	// Compute message symbols probability distributions
	qra64_fastfading_intrinsics(ix, rxen, hptr, hlen, noisestd, EsNoMetric, submode);

	// Step 3 ---------------------------------------------------------------------------
	// 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<QRA64_C;k++) {
		pd_init(dsttmp,srctmp,QRA64_M);
		dsttmp -=QRA64_M;
		srctmp -=QRA64_M;
		}
	// Initialize crc prob to a uniform distribution
	pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

	// Step 4 ---------------------------------------------------------------------------
	// Attempt to decode
	rc = qra64_decode_attempts(pcodec, xdec, ix);
	if (rc<0)
	  return rc;	// no success

	// copy decoded message (without crc) to output buffer
	memcpy(x,xdec,QRA64_K*sizeof(int));

	// Step 5 ----------------------------------------------------------------------------
	// Estimate the message Eb/No

	if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
		return rc;

	// reencode message to estimate Eb/No
	qra_encode(&QRA64_CODE, ydec, xdec);
	// puncture crc
	memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));

	// compute Es/N0 of decoded message
	esno = qra64_fastfading_msg_esno(ydec,rxen,noisestd, EsNoMetric, hlen,submode);

	// as the weigthing function include about 90% of the energy
	// we could compute the unbiased esno with:
	// esno = esno/0.9;

	// Es/N0 --> Eb/N0 conversion
	ebnoval = 1.0f/(1.0f*QRA64_K/QRA64_N*QRA64_m)*esno;

	// compute value in dB
	if (ebnoval<=0) {
	  ebnoval = EBNO_MIN; // assume a minimum, positive value
    }
	else {
	  ebnoval = 10.0f*(float)log10(ebnoval);
	  if (ebnoval<EBNO_MIN) {
		  ebnoval = EBNO_MIN;
      }
    }

	*ebno = ebnoval;

  return rc;

}

/*
int qra64_fastfading_channel(float **rxen, const int *xmsg, const int submode, const float EbN0dB, const float B90, const int fadingModel)
{
	// Simulate transmission over a fading channel and non coherent detection

	// Set rxen to point to an array of bin energies formatted as required
	// by the (fast-fading) decoding routine

	// returns 0 on success or negative values on error conditions

	static float *channel_out = NULL;
	static int    channel_submode = -1;

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int n,j,hidx, hlen;
	const float *hptr;
	float *cursym,*curtone;

	float iq[2];
	float *curi, *curq;
	float sigmasig[65];	// signal standard deviation taps

	float N0, EsN0, Es, sigmanoise;

	if (rxen==NULL)
		return -1;		// rxen must be a non-null ptr

	// allocate output buffer if not yet done or if submode changed
	if (channel_out==NULL || submode!=channel_submode) {

		// unallocate previous buffer
		if (channel_out)
			free(channel_out);

		// allocate new buffer
		// we allocate twice the mem so that we can store/compute complex amplitudes
		channel_out = (float*)malloc(bpm*sizeof(float)*2);
		if (channel_out==NULL)
			return -2;	// error allocating memory

		channel_submode = submode;
		}

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate energy weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian weighting taps
		hlen = glen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian weighting taps
		hlen = glen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// Compute the unfaded tone amplitudes from the Eb/No value passed to the call
	N0 = 1.0f;	// assume unitary noise PSD
	sigmanoise = (float)sqrt(N0/2);
	EsN0 = (float)pow(10.0f,EbN0dB/10.0f)*QRA64_m*QRA64_K/QRA64_N; // Es/No = m*R*Eb/No
	Es   = EsN0*N0;

	// compute signal bin sigmas
	for (n=0;n<hlen;n++)
		sigmasig[n] = (float)sqrt(hptr[n]*Es/2.0f);

	// Generate gaussian noise iq components
	normrnd_s(channel_out, bpm*2, 0 , sigmanoise);

	// Add symbols, bin by bin energies
	for (n=0;n<QRA64_N;n++) {

		cursym  = channel_out+n*bps + QRA64_M; // point to n-th symbol
		curtone = cursym+xmsg[n]*bpt;	 // point to encoded tone
		curi    = curtone-hlen+1;		 // point to real part of first bin
		curq    = curtone-hlen+1+bpm;	 // point to imag part of first bin

		// generate Rayleigh faded bins with given average energy and add to noise
		for (j=0;j<hlen;j++) {
			normrnd_s(iq, 2, 0 , sigmasig[j]);
			*curi++ += iq[0];
			*curq++ += iq[1];
			}
		for (j=hlen-2;j>=0;j--) {
			normrnd_s(iq, 2, 0 , sigmasig[j]);
			*curi++ += iq[0];
			*curq++ += iq[1];
			}

		}

	// compute total bin energies (S+N) and store in first half of buffer
	curi = channel_out;
	curq = channel_out+bpm;
	for (n=0;n<bpm;n++)
		channel_out[n] = curi[n]*curi[n] + curq[n]*curq[n];

	// set rxen to point to the channel output energies
	*rxen = channel_out;

	return 0;
}
*/

// Static functions definitions ----------------------------------------------

// fast-fading static functions --------------------------------------------------------------

static float qra64_fastfading_estim_noise_std(const float *rxen, const float esnometric, const int submode)
{
	// estimate the noise standard deviation from nominal frequency symbol bins
	// transform energies to amplitudes

	// rxen = message symbols energies (overwritten with symbols amplitudes)
	// esnometric = Es/No at nominal frequency bin for which we compute the decoder metric
	// submode = submode used (0=A...4=E)

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int k;
	float sigmaest;

	// estimate noise std
	sigmaest = 0;
	for (k=0;k<bpm;k++)
		sigmaest += rxen[k];

	sigmaest = sigmaest/bpm;
	sigmaest = (float)sqrt(sigmaest/(1.0f+esnometric/bps)/2.0f);

	// Note: sigma is overestimated by the (unknown) factor sqrt((1+esno(true)/bps)/(1+esnometric/bps))

	return sigmaest;
}

static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxen,
				const float *hptr,
				const int    hlen,
				const float sigmaest,
				const float EsNoMetric,
				const int submode)
{

	// For each symbol in a message:
	// a) Compute tones loglikelihoods as a sum of products between of the expected
	// energy fading coefficient and received energies.
	// b) Compute intrinsic symbols probability distributions from symbols loglikelihoods

	int n,k,j, bps, bpt;
	const float *cursym, *curbin;
	float *curix;
	float u, maxloglh, loglh, sumix,hh;
	float w[65];
	int hhsz  = hlen-1;
	int hlast = 2*hhsz;
	float npwrest = 2.0f*sigmaest*sigmaest;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	u = EsNoMetric;
	// compute weights from energy tables
	for (j=0;j<hlen;j++) {
		hh = hptr[j]*u;
		w[j] = hh/(1+hh)/npwrest;
		}

	for (n=0;n<QRA64_N;n++) {			// for each symbol in the message
		cursym = rxen+n*bps + QRA64_M;	// point to current symbol nominal bin
		maxloglh = 0;
		curix  = pix+n*QRA64_M;
		for (k=0;k<QRA64_M;k++) {   // for each tone in the current symbol
			curbin = cursym + k*bpt -hlen+1;
			// compute tone loglikelihood (symmetric fir with given weights)
			loglh = 0.f;
			for (j=0;j<hhsz;j++)
				loglh += w[j]*(curbin[j] + curbin[hlast-j]);
			loglh += w[hhsz]*curbin[hhsz];

			if (loglh>maxloglh)		// keep track of the max loglikelihood
				maxloglh = loglh;
			curix[k]=loglh;
			}

		// scale to likelihoods
		sumix = 0.f;
		for (k=0;k<QRA64_M;k++) {
			u = (float)exp(curix[k]-maxloglh);
			curix[k]=u;
			sumix +=u;
			}
		// scale to probabilities
		sumix = 1.0f/sumix;
		for (k=0;k<QRA64_M;k++)
			curix[k] = curix[k]*sumix;
		}
}

static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxen,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode)
{
	// Estimate msg Es/N0

	int n,j, bps, bpt;
	const float *cursym, *curtone, *curbin;
	float u, msgsn,esno;
	int tothlen = 2*hlen-1;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	msgsn = 0;
	for (n=0;n<QRA64_N;n++) {
		cursym  = rxen+n*bps + QRA64_M; // point to n-th symbol amplitudes
		curtone = cursym+ydec[n]*bpt;	 // point to decoded tone amplitudes
		curbin  = curtone-hlen+1;		 // point to first bin amplitude

		// sum bin energies
		for (j=0;j<tothlen;j++)
			msgsn += curbin[j];

		}

	msgsn =  msgsn/(QRA64_N*tothlen);	// avg msg energy per bin (noise included)

	// as sigma is overestimated (sigmatrue = sigma*sqrt((1+EsNoMetric/bps)/(1+EsNo/bps))
	// we have: msgsn = (1+x/hlen)/(1+x/bps)*2*sigma^2*(1+EsnoMetric/bps), where x = Es/N0(true)
	//
	// we can then write:
	// u = msgsn/2.0f/(sigma*sigma)/(1.0f+EsNoMetric/bps);
	// (1+x/hlen)/(1+x/bps) = u

	u = msgsn/(2.0f*sigma*sigma)/(1.0f+EsNoMetric/bps);

	// check u>1
	if (u<1)
		return 0.f;

	// check u<bps/tot hlen
	if (u>(bps/tothlen))
		return 10000.f;

	// solve for Es/No
	esno = (u-1.0f)/(1.0f/tothlen-u/bps);

	return esno;

}

#ifdef LIMIT_AP_MASKS

static int call1_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_CALL1   =   0xFFFFFFC
	int u = papmsg[4]^pdec[4];
	return (u&0x3C)==0;
}
static int call2_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_CALL2   =   0xFFFFFFC
	int u = papmsg[9]^pdec[9];
	return  (u&0x30)==0;
}
static int grid_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_GRIDFULL =	0x3FFC
	int u = papmsg[11]^pdec[11];
	int rc = (u&0x03)==0;

	u = papmsg[9]^pdec[9];

	return (u&0x0C)==0 && rc;
}

#else
#define call1_match(a,b) (1)
#define call2_match(a,b) (1)
#define grid_match(a,b) (1)
#endif

// Attempt to decode given intrisic information
static int qra64_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix)
{
  int rc;

  // Attempt to decode without a-priori info --------------------------------
  rc = qra64_do_decode(xdec, ix, NULL, NULL);
  if (rc>=0)
	  return 0; // successfull decode with AP0
  else
	  if (pcodec->apflags==QRA_NOAP)
		  // nothing more to do
		  return rc; // rc<0 = unsuccessful decode

  // Here we handle decoding with AP knowledge

  // Attempt to decode CQ calls
  rc = qra64_do_decode(xdec,ix,pcodec->apmask_cqqrz, pcodec->apmsg_cqqrz);
  if (rc>=0)
	  return 1;    // decoded [cq/qrz ? ?]

  rc = qra64_do_decode(xdec, ix, pcodec->apmask_cqqrz_ooo,
		       pcodec->apmsg_cqqrz);
  if (rc>=0)
	  // check that ooo really matches
	  if (grid_match(pcodec->apmsg_cqqrz,xdec))
		  return 2;    // decoded [cq/qrz ? ooo]

  // attempt to decode calls directed to us
  if (pcodec->apmsg_set[APTYPE_MYCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1,
		       pcodec->apmsg_call1);
	if (rc>=0)
		// check that mycall really matches
		if (call1_match(pcodec->apmsg_call1,xdec))
			return 3;    // decoded [mycall ? ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_ooo,
		       pcodec->apmsg_call1);
	if (rc>=0)
		// check that mycall and ooo really match
		if (call1_match(pcodec->apmsg_call1,xdec) &&
			grid_match(pcodec->apmsg_call1,xdec))
			return 4;    // decoded [mycall ? ooo]
	}

  // attempt to decode [mycall hiscall ?] msgs
  if (pcodec->apmsg_set[APTYPE_BOTHCALLS]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2,
		       pcodec->apmsg_call1_call2);
	if (rc>=0)
		// check that mycall and hiscall really match
		if (call1_match(pcodec->apmsg_call1_call2,xdec) &&
			call2_match(pcodec->apmsg_call1_call2,xdec))
			return 5;    // decoded [mycall srccall ?]
	}

  // attempt to decode [? hiscall ?/b] msgs
  if (pcodec->apmsg_set[APTYPE_HISCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2,
		       pcodec->apmsg_call2);
	if (rc>=0)
		// check that hiscall really match
		if (call2_match(pcodec->apmsg_call2,xdec))
			return 6;    // decoded [? hiscall ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2_ooo,
		       pcodec->apmsg_call2);
	if (rc>=0)
		// check that hiscall and ooo match
		if (call2_match(pcodec->apmsg_call2,xdec) &&
			grid_match(pcodec->apmsg_call2,xdec))
			return 7;    // decoded [? hiscall ooo]
	}

  // attempt to decode [cq/qrz hiscall ?/b/grid] msgs
  if (pcodec->apmsg_set[APTYPE_CQHISCALL]) {

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2,
				pcodec->apmsg_cq_call2);
	if (rc>=0)
		// check that hiscall matches
		if (call2_match(pcodec->apmsg_call2,xdec))
			return 9;	// decoded [cq/qrz hiscall ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2_grid,xdec, QRA64_K*sizeof(int))==0)
			return 11;		// decoded [cq/qrz hiscall grid]
		}

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2,xdec, QRA64_K*sizeof(int))==0)
			return 10;    // decoded [cq/qrz hiscall ]
		}
	}

  // attempt to decode [mycall hiscall grid]
  if (pcodec->apmsg_set[APTYPE_FULL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2_grid,
		       pcodec->apmsg_call1_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// All the three msg fields were given.
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_call1_call2_grid,xdec, QRA64_K*sizeof(int))==0)
			return 8;	   // decoded [mycall hiscall grid]
		}
	}

  // all decoding attempts failed
  return -1;
}

// Decode with given a-priori information
static int qra64_do_decode(int *xdec, 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];

  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

  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<sz;k++) {
    t = x[k];
    for (j=0;j<6;j++) {
      if ((t^sr)&0x01)
	sr = (sr>>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<QRA64_K;k++) {	// we can mask only information symbols distrib
    smask = mask[k];
    row = PD_ROWADDR(dst,QRA64_M,k);
    if (smask) {
      for (kk=0;kk<QRA64_M;kk++)
	if (((kk^x[k])&smask)!=0)
	  *(row+kk) = 0.f;

      pd_norm(row,QRA64_m);
    }
  }
}

// encode/decode msgs as done in JT65
void encodemsg_jt65(int *y, const int call1, const int call2, const int grid)
{
  y[0]= (call1>>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;
}


code:
// 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 <http://www.gnu.org/licenses/>.

#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
#define QRA_USERAP   2	// a-priori knowledge messages provided by the user

// 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

// Types of a-priori knowledge messages
#define APTYPE_CQQRZ     0  // [cq/qrz ?       ?/blank]
#define APTYPE_MYCALL    1  // [mycall ?       ?/blank]
#define APTYPE_HISCALL   2  // [?      hiscall ?/blank]
#define APTYPE_BOTHCALLS 3  // [mycall hiscall ?]
#define APTYPE_FULL	 4  // [mycall hiscall grid]
#define APTYPE_CQHISCALL 5  // [cq/qrz hiscall ?/blank]
#define APTYPE_SIZE	(APTYPE_CQHISCALL+1)

typedef struct {
  float decEsNoMetric;
  int apflags;
  int apmsg_set[APTYPE_SIZE];     // indicate which ap type knowledge has
                                  // been set by the user
// ap messages buffers
  int apmsg_cqqrz[12];		  // [cq/qrz ?       ?/blank]
  int apmsg_call1[12];		  // [mycall ?       ?/blank]
  int apmsg_call2[12];		  // [?      hiscall ?/blank]
  int apmsg_call1_call2[12];      // [mycall hiscall ?]
  int apmsg_call1_call2_grid[12]; // [mycall hiscall grid]
  int apmsg_cq_call2[12];         // [cq     hiscall ?/blank]
  int apmsg_cq_call2_grid[12];    // [cq     hiscall grid]

// ap messages masks
  int apmask_cqqrz[12];
  int apmask_cqqrz_ooo[12];
  int apmask_call1[12];
  int apmask_call1_ooo[12];
  int apmask_call2[12];
  int apmask_call2_ooo[12];
  int apmask_call1_call2[12];
  int apmask_call1_call2_grid[12];
  int apmask_cq_call2[12];
  int apmask_cq_call2_ooo[12];
} qra64codec;

#ifdef __cplusplus
extern "C" {
#endif

qra64codec *qra64_init(int flags);
// QRA64 mode initialization function
// arguments:
//    flags: set the decoder mode
//	     QRA_NOAP    use no a-priori information
//	     QRA_AUTOAP  use any relevant previous decodes
//	     QRA_USERAP  use a-priori information provided via qra64_apset(...)
// returns:
//    Pointer to initialized qra64codec data structure
//		this pointer 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 encoder
// arguments:
//    pcodec = pointer to a qra64codec data structure as returned by qra64_init
//    x      = pointer to the message to be encoded, int x[12]
//	       x must point to an array of integers (i.e. defined as int x[12])
//    y      = pointer to encoded message, int y[63]=
// ----------------------------------------------------------------------------

int  qra64_decode(qra64codec *pcodec, float *ebno, int *x, const float *r);
// QRA64 mode decoder
// arguments:
//    pcodec = pointer to a qra64codec data structure as returned by qra64_init
//    ebno   = pointer to a float where the avg Eb/No (in dB) will be stored
//             in case of successfull decoding
//             (pass a null pointer if not interested)
//    x      = pointer to decoded message, int x[12]
//    r      = pointer to received symbol energies (squared amplitudes)
//             r must point to an array of length QRA64_M*QRA64_N (=64*63=4032)
//	       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
//		rc=6    [?    CALL ?]    AP29
//		rc=7    [?    CALL  ]    AP45
//		rc=8    [CALL CALL GRID] AP72 (actually a AP68 mask to reduce false decodes)
//		rc=9    [CQ   CALL ?]    AP55
//		rc=10   [CQ   CALL  ]    AP70 (actaully a AP68 mask to reduce false decodes)

//  return codes in the range 1-10 indicate the amount and the type of a-priori
//  information was required to decode the received message.

// Decode a QRA64 msg using a fast-fading metric
int qra64_decode_fastfading(
				qra64codec *pcodec,		// ptr to the codec structure
				float *ebno,			// ptr to where the estimated Eb/No value will be saved
				int *x,					// ptr to decoded message
				const float *rxen,		// ptr to received symbol energies array
				const int submode,		// submode idx (0=QRA64A ... 4=QRA64E)
				const float B90,	    // spread bandwidth (90% fractional energy)
				const int fadingModel); // 0=Gaussian 1=Lorentzian fade model
//
// rxen: The array of the received bin energies
//       Bins must be spaced by integer multiples of the symbol rate (1/Ts Hz)
//       The array must be an array of total length U = L x N where:
//			L: is the number of frequency bins per message symbol (see after)
//          N: is the number of symbols in a QRA64 msg (63)
//
//       The number of bins/symbol L depends on the selected submode accordingly to
//		 the following rule:
//			L = (64+64*2^submode+64) = 64*(2+2^submode)
//		 Tone 0 is always supposed to be at offset 64 in the array.
//		 The m-th tone nominal frequency is located at offset 64 + m*2^submode (m=0..63)
//
//		 Submode A: (2^submode = 1)
//          L = 64*3 = 196 bins/symbol
//          Total length of the energies array: U = 192*63 = 12096 floats
//
//		 Submode B: (2^submode = 2)
//          L = 64*4 = 256 bins/symbol
//          Total length of the energies array: U = 256*63 = 16128 floats
//
//		 Submode C: (2^submode = 4)
//          L = 64*6 = 384 bins/symbol
//          Total length of the energies array: U = 384*63 = 24192 floats
//
//		 Submode D: (2^submode = 8)
//          L = 64*10 = 640 bins/symbol
//          Total length of the energies array: U = 640*63 = 40320 floats
//
//		 Submode E: (2^submode = 16)
//          L = 64*18 = 1152 bins/symbol
//          Total length of the energies array: U = 1152*63 = 72576 floats
//
//		Note: The rxen array is modified and reused for internal calculations.
//
//
//	B90: spread fading bandwidth in Hz (90% fractional average energy)
//
//			B90 should be in the range 1 Hz ... 238 Hz
//			The value passed to the call is rounded to the closest value among the
//			64 available values:
//				B = 1.09^k Hz, with k=0,1,...,63
//
//			I.e. B90=27 Hz will be approximated in this way:
//				k = rnd(log(27)/log(1.09)) = 38
//              B90 = 1.09^k = 1.09^38 = 26.4 Hz
//
//          For any input value the maximum rounding error is not larger than +/- 5%
//
// return codes: same return codes of qra64_decode (+some additional error codes)

// Simulate the fast-fading channel (to be used with qra64_decode_fastfading)
int qra64_fastfading_channel(
				float **rxen,
				const int *xmsg,
				const int submode,
				const float EbN0dB,
				const float B90,
				const int fadingModel);
// Simulate transmission over a fading channel with given B90, fading model and submode
// and non coherent detection.
// Sets rxen to point to an array of bin energies formatted as required
// by the (fast-fading) decoding routine.
// returns 0 on success or negative values on error conditions

int qra64_apset(qra64codec *pcodec, const int mycall, const int hiscall, const int grid, const int aptype);
// Set decoder a-priori knowledge accordingly to the type of the message to
// look up for
// arguments:
//    pcodec    = pointer to a qra64codec data structure as returned by qra64_init
//    mycall    = mycall to look for
//    hiscall   = hiscall to look for
//    grid      = grid to look for
//    aptype    = define the type of AP to be set:
//		APTYPE_CQQRZ     set [cq/qrz ?       ?/blank]
//		APTYPE_MYCALL    set [mycall ?       ?/blank]
//		APTYPE_HISCALL   set [?      hiscall ?/blank]
//		APTYPE_BOTHCALLS set [mycall hiscall ?]
//		APTYPE_FULL		 set [mycall hiscall grid]
//		APTYPE_CQHISCALL set [cq/qrz hiscall ?/blank]

// returns:
//	 0  on success
//  -1  when qra64_init was called with the QRA_NOAP flag
//	-2  invalid apytpe (valid range [APTYPE_CQQRZ..APTYPE_CQHISCALL]
//	    (APTYPE_CQQRZ [cq/qrz ? ?] is set by default )

void qra64_apdisable(qra64codec *pcodec, const int aptype);
// disable specific AP type
// arguments:
//    pcodec    = pointer to a qra64codec data structure as returned by qra64_init
//    aptype    = define the type of AP to be disabled
//		APTYPE_CQQRZ     disable [cq/qrz   ?      ?/blank]
//		APTYPE_MYCALL    disable [mycall   ?      ?/blank]
//		APTYPE_HISCALL   disable [   ?   hiscall  ?/blank]
//		APTYPE_BOTHCALLS disable [mycall hiscall     ?   ]
//		APTYPE_FULL	 disable [mycall hiscall     grid]
//		APTYPE_CQHISCALL set [cq/qrz hiscall ?/blank]

void qra64_close(qra64codec *pcodec);
// Free memory allocated by qra64_init
// arguments:
//    pcodec = pointer to a qra64codec data structure as returned by qra64_init

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

// 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_

license: mit


code:
/*
qra64.c
Encoding/decoding functions for the QRA64 mode

(c) 2016 - Nico Palermo, IV3NWV

-------------------------------------------------------------------------------

   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/>.

-----------------------------------------------------------------------------

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

#include "qra64.h"
#include "../qracodes/qracodes.h"
#include "../qracodes/qra13_64_64_irr_e.h"
#include "../qracodes/pdmath.h"
#include "../qracodes/normrnd.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_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix);
static int qra64_do_decode(int *x, const float *pix, const int *ap_mask,
			    const int *ap_x);
static float qra64_fastfading_estim_noise_std(
				float *rxen,
				const float esnometric,
				const int submode);
static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxamp,
				const float *hptr,
				const int    hlen,
				const float cmetric,
				const int submode);
static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxamp,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode);

// 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_GRIDFULL12	0x3FFC	  // less aggressive mask (to be used with full AP decoding)
#define MASK_GRIDBIT	0x8000	  // b[15] is 1 for free text, 0 otherwise
// ----------------------------------------------------------------------------

qra64codec *qra64_init(int flags)
{

  // 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;

  memset(pcodec->apmsg_set,0,APTYPE_SIZE*sizeof(int));

  if (flags==QRA_NOAP)
    return pcodec;

  // for QRA_USERAP and QRA_AUTOAP modes we always enable [CQ/QRZ ? ?] mgs look-up.
  // encode CQ/QRZ AP messages
  // 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
  qra64_apset(pcodec, CALL_CQ, 0, GRID_BLANK, APTYPE_CQQRZ);

  // initialize masks for decoding with a-priori information
  encodemsg_jt65(pcodec->apmask_cqqrz,     MASK_CQQRZ, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cqqrz_ooo, MASK_CQQRZ, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1,     MASK_CALL1, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_ooo, MASK_CALL1, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call2,     0, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call2_ooo, 0, MASK_CALL2, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1_call2,     MASK_CALL1,MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_call2_grid,MASK_CALL1,MASK_CALL2, MASK_GRIDFULL12);
  encodemsg_jt65(pcodec->apmask_cq_call2,     MASK_CQQRZ, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cq_call2_ooo, MASK_CQQRZ, MASK_CALL2, MASK_GRIDFULL12);

  return pcodec;
}

void qra64_close(qra64codec *pcodec)
{
	free(pcodec);
}

int qra64_apset(qra64codec *pcodec, const int mycall, const int hiscall, const int grid, const int aptype)
{
// Set decoder a-priori knowledge accordingly to the type of the message to look up for
// arguments:
//		pcodec    = pointer to a qra64codec data structure as returned by qra64_init
//		mycall    = mycall to look for
//		hiscall   = hiscall to look for
//		grid      = grid to look for
//		aptype    = define and masks the type of AP to be set accordingly to the following:
//			APTYPE_CQQRZ     set [cq/qrz ?       ?/blank]
//			APTYPE_MYCALL    set [mycall ?       ?/blank]
//			APTYPE_HISCALL   set [?      hiscall ?/blank]
//			APTYPE_BOTHCALLS set [mycall hiscall ?]
//			APTYPE_FULL		 set [mycall hiscall grid]
//			APTYPE_CQHISCALL set [cq/qrz hiscall ?/blank] and [cq/qrz hiscall grid]
// returns:
//		0   on success
//      -1  when qra64_init was called with the QRA_NOAP flag
//		-2  invalid apytpe

	if (pcodec->apflags==QRA_NOAP)
		return -1;

	switch (aptype) {
		case APTYPE_CQQRZ:
			encodemsg_jt65(pcodec->apmsg_cqqrz,  CALL_CQ, 0, GRID_BLANK);
			break;
		case APTYPE_MYCALL:
			encodemsg_jt65(pcodec->apmsg_call1,  mycall,  0, GRID_BLANK);
			break;
		case APTYPE_HISCALL:
			encodemsg_jt65(pcodec->apmsg_call2,  0, hiscall, GRID_BLANK);
			break;
		case APTYPE_BOTHCALLS:
			encodemsg_jt65(pcodec->apmsg_call1_call2,  mycall, hiscall, GRID_BLANK);
			break;
		case APTYPE_FULL:
			encodemsg_jt65(pcodec->apmsg_call1_call2_grid,  mycall, hiscall, grid);
			break;
		case APTYPE_CQHISCALL:
			encodemsg_jt65(pcodec->apmsg_cq_call2,      CALL_CQ, hiscall, GRID_BLANK);
			encodemsg_jt65(pcodec->apmsg_cq_call2_grid, CALL_CQ, hiscall, grid);
			break;
		default:
			return -2;	// invalid ap type
		}

	  pcodec->apmsg_set[aptype]=1;	// signal the decoder to look-up for the specified type

	  return 0;
}
void qra64_apdisable(qra64codec *pcodec, const int aptype)
{
	if (pcodec->apflags==QRA_NOAP)
		return;

	if (aptype<APTYPE_CQQRZ || aptype>=APTYPE_SIZE)
		return;

	pcodec->apmsg_set[aptype] = 0;	//  signal the decoder not to look-up to the specified type
}

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 hiscall,mycall,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;

  // Here we handle the QRA_AUTOAP mode --------------------------------------------

  // When a [hiscall mycall ?] msg is detected we instruct the decoder
  // to look for [mycall hiscall ?] msgs
  // otherwise when a [cq mycall ?] msg is sent we reset the APTYPE_BOTHCALLS

  // 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, nothing to do

  // It's a [hiscall mycall grid] message

  // We assume that mycall is our call (but we don't check it)
  // hiscall the station we are calling or a general call (CQ/QRZ/etc..)
  decodemsg_jt65(&hiscall,&mycall,&grid,x);

  if ((hiscall>=CALL_CQ && hiscall<=CALL_CQ999) || hiscall==CALL_CQDX ||
      hiscall==CALL_DE) {
	// tell the decoder to look for msgs directed to us
	qra64_apset(pcodec,mycall,0,0,APTYPE_MYCALL);
    // We are making a general call and don't know who might reply
    // Reset APTYPE_BOTHCALLS so decoder won't look for [mycall hiscall ?] msgs
    qra64_apdisable(pcodec,APTYPE_BOTHCALLS);
  } else {
    // We are replying to someone named hiscall
    // Set APTYPE_BOTHCALLS so decoder will try for [mycall hiscall ?] msgs
    qra64_apset(pcodec,mycall, hiscall, GRID_BLANK, APTYPE_BOTHCALLS);
  }

}

#define EBNO_MIN -10.0f		// minimum Eb/No value returned by the decoder (in dB)
int qra64_decode(qra64codec *pcodec, float *ebno, int *x, const float *rxen)
{
  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise variance
  float msge;						// estimated message energy
  float ebnoval;					// estimated Eb/No
  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
  noisestd = 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<QRA64_C;k++) {
    pd_init(dsttmp,srctmp,QRA64_M);
    dsttmp -=QRA64_M;
    srctmp -=QRA64_M;
  }
  // Initialize crc prob to a uniform distribution
  pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

  // Try to decode using all AP cases if required
  rc = qra64_decode_attempts(pcodec, xdec, ix);

  if (rc<0)
	  return rc;	// no success

  // successfull decode --------------------------------

  // copy decoded message (without crc) to output buffer
  memcpy(x,xdec,QRA64_K*sizeof(int));

  if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
	  return rc;

  // reencode message and estimate Eb/No
  qra_encode(&QRA64_CODE, ydec, xdec);
  // puncture crc
  memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));
  // compute total power of decoded message
  msge = 0;
  for (k=0;k<QRA64_N;k++) {
	  msge +=rxen[ydec[k]];	// add energy of current symbol
	  rxen+=QRA64_M;			// ptr to next symbol
	  }

  // NOTE:
  // To make a more accurate Eb/No estimation we should compute the noise variance
  // on all the rxen values but the transmitted symbols.
  // Noisestd is compute by qra_mfskbesselmetric assuming that
  // the signal power is much less than the total noise power in the QRA64_M tones
  // but this is true only if the Eb/No is low.
  // Here, in order to improve accuracy, we linearize the estimated Eb/No value empirically
  // (it gets compressed when it is very high as in this case the noise variance
  // is overestimated)

  // this would be the exact value if the noisestd were not overestimated at high Eb/No
  ebnoval = (0.5f/(QRA64_K*QRA64_m))*msge/(noisestd*noisestd)-1.0f;

  // Empirical linearization (to remove the noise variance overestimation)
  // the resulting SNR is accurate up to +20 dB (51 dB Eb/No)
  if (ebnoval>57.004f)
	  ebnoval=57.004f;
  ebnoval = ebnoval*57.03f/(57.03f-ebnoval);

  // compute value in dB
  if (ebnoval<=0)
	  ebnoval = EBNO_MIN; // assume a minimum, positive value
  else
	  ebnoval = 10.0f*(float)log10(ebnoval);
	  if (ebnoval<EBNO_MIN)
		  ebnoval = EBNO_MIN;

  *ebno = ebnoval;

  return rc;
}

// Tables of fading amplitudes coefficients for QRA64 (Ts=6912/12000)
// As the fading is assumed to be symmetric around the nominal frequency
// only the leftmost and the central coefficient are stored in the tables.
// (files have been generated with the Matlab code efgengauss.m and efgenlorentz.m)
#include "fadampgauss.c"
#include "fadamplorentz.c"

int qra64_decode_fastfading(
				qra64codec *pcodec,		// ptr to the codec structure
				float *ebno,			// ptr to where the estimated Eb/No value will be saved
				int *x,					// ptr to decoded message
				float *rxen,		    // ptr to received symbol energies array
				const int submode,		// submode idx (0=QRA64A ... 4=QRA64E)
				const float B90,	    // spread bandwidth (90% fractional energy)
				const int fadingModel)  // 0=Gaussian 1=Lorentzian fade model

// Decode a QRA64 msg using a fast-fading metric
//
// rxen: The array of the received bin energies
//       Bins must be spaced by integer multiples of the symbol rate (1/Ts Hz)
//       The array must be an array of total length U = L x N where:
//			L: is the number of frequency bins per message symbol (see after)
//          N: is the number of symbols in a QRA64 msg (63)

//       The number of bins/symbol L depends on the selected submode accordingly to
//		 the following rule:
//			L = (64+64*2^submode+64) = 64*(2+2^submode)
//		 Tone 0 is always supposed to be at offset 64 in the array.
//		 The m-th tone nominal frequency is located at offset 64 + m*2^submode (m=0..63)
//
//		 Submode A: (2^submode = 1)
//          L = 64*3 = 196 bins/symbol
//          Total length of the energies array: U = 192*63 = 12096 floats
//
//		 Submode B: (2^submode = 2)
//          L = 64*4 = 256 bins/symbol
//          Total length of the energies array: U = 256*63 = 16128 floats
//
//		 Submode C: (2^submode = 4)
//          L = 64*6 = 384 bins/symbol
//          Total length of the energies array: U = 384*63 = 24192 floats
//
//		 Submode D: (2^submode = 8)
//          L = 64*10 = 640 bins/symbol
//          Total length of the energies array: U = 640*63 = 40320 floats
//
//		 Submode E: (2^submode = 16)
//          L = 64*18 = 1152 bins/symbol
//          Total length of the energies array: U = 1152*63 = 72576 floats
//
//		Note: The rxen array is modified and reused for internal calculations.
//
//
//	B90: spread fading bandwidth in Hz (90% fractional average energy)
//
//			B90 should be in the range 1 Hz ... 238 Hz
//			The value passed to the call is rounded to the closest value among the
//			64 available values:
//				B = 1.09^k Hz, with k=0,1,...,63
//
//			I.e. B90=27 Hz will be approximated in this way:
//				k = rnd(log(27)/log(1.09)) = 38
//              B90 = 1.09^k = 1.09^38 = 26.4 Hz
//
//          For any input value the maximum rounding error is not larger than +/- 5%
//

{

  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise std
  float esno,ebnoval;				// estimated Eb/No
  float tempf;
  float EsNoMetric, cmetric;
  int rc;
  int hidx, hlen;
  const float *hptr;

	if (QRA64_NMSG!=QRA64_CODE.NMSG)
		return -16;					// QRA64_NMSG define is wrong

	if (submode<0 || submode>4)
		return -17;				// invalid submode

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate amplitude weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian weighting taps
		hlen = hlen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = hptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian weighting taps
		hlen = hlen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = hptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// compute (euristically) the optimal decoder metric accordingly the given spread amount
	// We assume that the decoder threshold is:
	//		Es/No(dB) = Es/No(AWGN)(dB) + 8*log(B90)/log(240)(dB)
	// that's to say, at the maximum Doppler spread bandwidth (240 Hz) there's a ~8 dB Es/No degradation
	// over the AWGN case
	tempf = 8.0f*(float)log((float)B90)/(float)log(240.0f);
	EsNoMetric = pcodec->decEsNoMetric*(float)pow(10.0f,tempf/10.0f);

	// Step 1 -----------------------------------------------------------------------------------
	// Evaluate the noise stdev from the received energies at nominal tone frequencies
    // and transform energies to amplitudes
	tempf = hptr[hlen-1];				// amplitude weigth at nominal freq;
	tempf = tempf*tempf;				// fractional energy at nominal freq. bin

	noisestd = qra64_fastfading_estim_noise_std(rxen, EsNoMetric, submode);
	cmetric = (float)sqrt(M_PI_2*EsNoMetric)/noisestd;

	// Step 2 -----------------------------------------------------------------------------------
	// Compute message symbols probability distributions
	qra64_fastfading_intrinsics(ix, rxen, hptr, hlen, cmetric, submode);

	// Step 3 ---------------------------------------------------------------------------
	// 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<QRA64_C;k++) {
		pd_init(dsttmp,srctmp,QRA64_M);
		dsttmp -=QRA64_M;
		srctmp -=QRA64_M;
		}
	// Initialize crc prob to a uniform distribution
	pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

	// Step 4 ---------------------------------------------------------------------------
	// Attempt to decode
	rc = qra64_decode_attempts(pcodec, xdec, ix);
	if (rc<0)
	  return rc;	// no success

	// copy decoded message (without crc) to output buffer
	memcpy(x,xdec,QRA64_K*sizeof(int));

	// Step 5 ----------------------------------------------------------------------------
	// Estimate the message Eb/No

	if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
		return rc;

	// reencode message to estimate Eb/No
	qra_encode(&QRA64_CODE, ydec, xdec);
	// puncture crc
	memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));

	// compute Es/N0 of decoded message
	esno = qra64_fastfading_msg_esno(ydec,rxen,noisestd, EsNoMetric, hlen,submode);

	// as the weigthing function include about 90% of the energy
	// we could compute the unbiased esno with:
	// esno = esno/0.9;

	// this would be the exact value if the noisestd were not overestimated at high Eb/No
	ebnoval = 1.0f/(1.0f*QRA64_K/QRA64_N*QRA64_m)*esno;

	// compute value in dB
	if (ebnoval<=0)
	  ebnoval = EBNO_MIN; // assume a minimum, positive value
	else
	  ebnoval = 10.0f*(float)log10(ebnoval);
	  if (ebnoval<EBNO_MIN)
		  ebnoval = EBNO_MIN;

	*ebno = ebnoval;

  return rc;

}

int qra64_fastfading_channel(float **rxen, const int *xmsg, const int submode, const float EbN0dB, const float B90, const int fadingModel)
{
	// Simulate transmission over a fading channel and non coherent detection

	// Set rxen to point to an array of bin energies formatted as required
	// by the (fast-fading) decoding routine

	// returns 0 on success or negative values on error conditions

	static float *channel_out = NULL;
	static int    channel_submode = -1;

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int n,j,hidx, hlen;
	const float *hptr;
	float *cursym,*curtone;

	float iq[2];
	float *curi, *curq;

//	float tote=0;	// debug

	float N0, EsN0, Es, A, sigmanoise, sigmasig;

	if (rxen==NULL)
		return -1;		// rxen must be a non-null ptr

	// allocate output buffer if not yet done or if submode changed
	if (channel_out==NULL || submode!=channel_submode) {

		// unallocate previous buffer
		if (channel_out)
			free(channel_out);

		// allocate new buffer
		// we allocate twice the mem so that we can store/compute complex amplitudes
		channel_out = (float*)malloc(bpm*sizeof(float)*2);
		if (channel_out==NULL)
			return -2;	// error allocating memory

		channel_submode = submode;
		}

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate amplitude weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian weighting taps
		hlen = hlen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = hptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian weighting taps
		hlen = hlen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = hptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// Compute the unfaded tone amplitudes from the Eb/No value passed to the call
	N0 = 1.0f;	// assume unitary noise PSD
	sigmanoise = (float)sqrt(N0/2);
	EsN0 = (float)pow(10.0f,EbN0dB/10.0f)*QRA64_m*QRA64_K/QRA64_N; // Es/No = m*R*Eb/No
	Es   = EsN0*N0;
	A    = (float)sqrt(Es/2.0f);	// unfaded tone amplitude (i^2+q^2 = Es/2+Es/2 = Es)

	// Generate gaussian noise iq components
	normrnd_s(channel_out, bpm*2, 0 , sigmanoise);

	// Add message symbols energies
	for (n=0;n<QRA64_N;n++) {

		cursym  = channel_out+n*bps + QRA64_M; // point to n-th symbol
		curtone = cursym+xmsg[n]*bpt;	 // point to encoded tone
		curi    = curtone-hlen+1;		 // point to real part of first bin
		curq    = curtone-hlen+1+bpm;	 // point to imag part of first bin

		// generate Rayleigh faded bins with given average energy and add to noise
		for (j=0;j<hlen;j++) {
			sigmasig = A*hptr[j];
			normrnd_s(iq, 2, 0 , sigmasig);
//			iq[0]=sigmasig*sqrt(2); iq[1]=0;	debug: used to verify Eb/No
			*curi++ += iq[0];
			*curq++ += iq[1];
//			tote +=iq[0]*iq[0]+iq[1]*iq[1];		// debug
			}
		for (j=hlen-2;j>=0;j--) {
			sigmasig = A*hptr[j];
			normrnd_s(iq, 2, 0 , sigmasig);
//			iq[0]=sigmasig*sqrt(2); iq[1]=0;	debug: used to verify Eb/No
			*curi++ += iq[0];
			*curq++ += iq[1];
//			tote +=iq[0]*iq[0]+iq[1]*iq[1];		// debug
			}

		}

//	tote = tote/QRA64_N;	// debug

	// compute total bin energies (S+N) and store in first half of buffer
	curi = channel_out;
	curq = channel_out+bpm;
	for (n=0;n<bpm;n++)
		channel_out[n] = curi[n]*curi[n] + curq[n]*curq[n];

	// set rxen to point to the channel output energies
	*rxen = channel_out;

	return 0;
}

// Static functions definitions ----------------------------------------------

// fast-fading static functions --------------------------------------------------------------

static float qra64_fastfading_estim_noise_std(float *rxen, const float esnometric, const int submode)
{
	// estimate the noise standard deviation from nominal frequency symbol bins
	// transform energies to amplitudes

	// rxen = message symbols energies (overwritten with symbols amplitudes)
	// esnometric = Es/No at nominal frequency bin for which we compute the decoder metric
	// submode = submode used (0=A...4=E)

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int k;
	float sigmaest;

	// estimate noise std
	sigmaest = 0;
	for (k=0;k<bpm;k++) {
		sigmaest += rxen[k];
		// convert energies to amplitudes for later use
		rxen[k] = (float)sqrt(rxen[k]);	// we do it in place to avoid memory allocations
		}
	sigmaest = sigmaest/bpm;
	sigmaest = (float)sqrt(sigmaest/(1.0f+esnometric/bps)/2.0f);

	// Note: sigma is overestimated by the (unknown) factor sqrt((1+esno(true)/bps)/(1+esnometric/bps))

	return sigmaest;
}

static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxamp,
				const float *hptr,
				const int    hlen,
				const float cmetric,
				const int submode)
{

	// For each symbol in a message:
	// a) Compute tones loglikelihoods as a sum of products between of the expected
	// amplitude fading coefficient and received amplitudes.
	// Each product is computed as log(I0(hk*xk*cmetric)) where hk is the average fading amplitude,
	// xk is the received amplitude at bin offset k, and cmetric is a constant dependend on the
	// Eb/N0 value for which the metric is optimized
	// The function y = log(I0(x)) is approximated as y = x^2/(x+e)
	// b) Compute intrinsic symbols probability distributions from symbols loglikelihoods

	int n,k,j, bps, bpt;
	const float *cursym, *curbin;
	float *curix;
	float u, maxloglh, loglh, sumix;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	for (n=0;n<QRA64_N;n++) {			// for each symbol in the message
		cursym = rxamp+n*bps + QRA64_M;	// point to current symbol nominal bin
		maxloglh = 0;
		curix  = pix+n*QRA64_M;
		for (k=0;k<QRA64_M;k++) {   // for each tone in the current symbol
			curbin = cursym + k*bpt -hlen+1;	// ptr to lowest bin of the current tone
			// compute tone loglikelihood as a weighted sum of bins loglikelihoods
			loglh = 0.f;
			for (j=0;j<hlen;j++) {
				u = *curbin++ * hptr[j]*cmetric;
				u = u*u/(u+(float)M_E);	// log(I0(u)) approx.
				loglh = loglh + u;
				}
			for (j=hlen-2;j>=0;j--) {
				u = *curbin++ * hptr[j]*cmetric;
				u = u*u/(u+(float)M_E);	// log(I0(u)) approx.
				loglh = loglh + u;
				}
			if (loglh>maxloglh)		// keep track of the max loglikelihood
				maxloglh = loglh;
			curix[k]=loglh;
			}
		// scale to likelihoods
		sumix = 0.f;
		for (k=0;k<QRA64_M;k++) {
			u = (float)exp(curix[k]-maxloglh);
			curix[k]=u;
			sumix +=u;
			}
		// scale to probabilities
		sumix = 1.0f/sumix;
		for (k=0;k<QRA64_M;k++)
			curix[k] = curix[k]*sumix;
		}
}

static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxamp,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode)
{
	// Estimate msg Es/N0

	int n,j, bps, bpt;
	const float *cursym, *curtone, *curbin;
	float u, msgsn,esno;
	int tothlen = 2*hlen-1;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	msgsn = 0;
	for (n=0;n<QRA64_N;n++) {
		cursym  = rxamp+n*bps + QRA64_M; // point to n-th symbol amplitudes
		curtone = cursym+ydec[n]*bpt;	 // point to decoded tone amplitudes
		curbin  = curtone-hlen+1;		 // point to first bin amplitude

		// sum bin energies
		for (j=0;j<hlen;j++) {
			u = *curbin++;
			msgsn += u*u;
			}
		for (j=hlen-2;j>=0;j--) {
			u = *curbin++;
			msgsn += u*u;
			}

		}

	msgsn =  msgsn/(QRA64_N*tothlen);	// avg msg energy per bin (noise included)

	// as sigma is overestimated (sigmatrue = sigma*sqrt((1+EsNoMetric/bps)/(1+EsNo/bps))
	// we have: msgsn = (1+x/hlen)/(1+x/bps)*2*sigma^2*(1+EsnoMetric/bps), where x = Es/N0(true)
	//
	// we can then write:
	// u = msgsn/2.0f/(sigma*sigma)/(1.0f+EsNoMetric/bps);
	// (1+x/hlen)/(1+x/bps) = u

	u = msgsn/(2.0f*sigma*sigma)/(1.0f+EsNoMetric/bps);

	// check u>1
	if (u<1)
		return 0.f;

	// check u<bps/tot hlen
	if (u>(bps/tothlen))
		return 10000.f;

	// solve for Es/No
	esno = (u-1.0f)/(1.0f/tothlen-u/bps);

	return esno;

}

// Attempt to decode given intrisic information
static int qra64_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix)
{
  int rc;

  // Attempt to decode without a-priori info --------------------------------
  rc = qra64_do_decode(xdec, ix, NULL, NULL);
  if (rc>=0)
	  return 0; // successfull decode with AP0
  else
	  if (pcodec->apflags==QRA_NOAP)
		  // nothing more to do
		  return rc; // rc<0 = unsuccessful decode

  // Here we handle decoding with AP knowledge

  // Attempt to decode CQ calls
  rc = qra64_do_decode(xdec,ix,pcodec->apmask_cqqrz, pcodec->apmsg_cqqrz);
  if (rc>=0) return 1;    // decoded [cq/qrz ? ?]

  rc = qra64_do_decode(xdec, ix, pcodec->apmask_cqqrz_ooo,
		       pcodec->apmsg_cqqrz);
  if (rc>=0) return 2;    // decoded [cq ? ooo]

  // attempt to decode calls directed to us
  if (pcodec->apmsg_set[APTYPE_MYCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1,
		       pcodec->apmsg_call1);
	if (rc>=0) return 3;    // decoded [mycall ? ?]
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_ooo,
		       pcodec->apmsg_call1);
	if (rc>=0) return 4;    // decoded [mycall ? ooo]
	}

  // attempt to decode [mycall srccall ?] msgs
  if (pcodec->apmsg_set[APTYPE_BOTHCALLS]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2,
		       pcodec->apmsg_call1_call2);
	if (rc>=0) return 5;    // decoded [mycall srccall ?]
	}

  // attempt to decode [? hiscall ?/b] msgs
  if (pcodec->apmsg_set[APTYPE_HISCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2,
		       pcodec->apmsg_call2);
	if (rc>=0) return 6;    // decoded [? hiscall ?]
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2_ooo,
		       pcodec->apmsg_call2);
	if (rc>=0) return 7;    // decoded [? hiscall ooo]
	}

  // attempt to decode [cq/qrz hiscall ?/b/grid] msgs
  if (pcodec->apmsg_set[APTYPE_CQHISCALL]) {

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2,
				pcodec->apmsg_cq_call2);
	if (rc>=0) return 9;	// decoded [cq/qrz hiscall ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2_grid,xdec, QRA64_K*sizeof(int))!=0)
			return -1;
		else
			return 11;		// decoded [cq/qrz hiscall grid]
		};

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2,xdec, QRA64_K*sizeof(int))!=0)
			return -1;
		else
			return 10;    // decoded [cq/qrz hiscall ]
		}
	}

  // attempt to decode [mycall hiscall grid]
  if (pcodec->apmsg_set[APTYPE_FULL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2_grid,
		       pcodec->apmsg_call1_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// All the three msg fields were given.
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_call1_call2_grid,xdec, QRA64_K*sizeof(int))!=0)
			return -1;
		else
			return 8;	   // decoded [mycall hiscall grid]
		}
	}

  // all decoding attempts failed
  return rc;
}

// Decode with given a-priori information
static int qra64_do_decode(int *xdec, 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];

  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

  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<sz;k++) {
    t = x[k];
    for (j=0;j<6;j++) {
      if ((t^sr)&0x01)
	sr = (sr>>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<QRA64_K;k++) {	// we can mask only information symbols distrib
    smask = mask[k];
    row = PD_ROWADDR(dst,QRA64_M,k);
    if (smask) {
      for (kk=0;kk<QRA64_M;kk++)
	if (((kk^x[k])&smask)!=0)
	  *(row+kk) = 0.f;

      pd_norm(row,QRA64_m);
    }
  }
}

// encode/decode msgs as done in JT65
void encodemsg_jt65(int *y, const int call1, const int call2, const int grid)
{
  y[0]= (call1>>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;
}


code:
// qra64_subs.c
// Fortran interface routines for QRA64

#include "qra64.h"
#include <stdio.h>

static qra64codec *pqra64codec = NULL;

void qra64_enc_(int x[], int y[])
{
  if (pqra64codec==NULL) pqra64codec = qra64_init(QRA_USERAP);
  qra64_encode(pqra64codec, y, x);
}

void qra64_dec_(float r[], int* nc1, int* nc2, int* ng2, int* APtype,
		int* iset, int* ns0, float* b0, int* nf0,
		int xdec[], float* snr, int* rc)
{
/*
  APtype:                         AP
-----------------------------------------------------------------------
     -1                            0     (no AP information)
      0    [CQ/QRZ    ?      ? ]  25/37
      1    [MyCall    ?      ? ]  25/37
      2    [  ?    HisCall   ? ]  25/37
      3    [MyCall HisCall   ? ]  49/68
      4    [MyCall HisCall grid]  68
      5    [CQ/QRZ HisCall   ? ]  49/68

     rc    Message format         AP APTYPE Comments
------------------------------------------------------------------------
    -16                                     Failed sanity check
     -2                                     Decoded but CRC failed
     -1                                     No decode
      0    [   ?      ?      ? ]   0   -1   Decode with no AP info
      1    [CQ/QRZ    ?      ? ]  25    0
      2    [CQ/QRZ    ?      _ ]  37    0
      3    [MyCall    ?      ? ]  25    1
      4    [MyCall    ?      _ ]  37    1
      5    [MyCall HisCall   ? ]  49    3
      6    [   ?   HisCall   ? ]  25    2   Optional
      7    [   ?   HisCall   _ ]  37    2   Optional
      8    [MyCall HisCall Grid]  68    4
      9    [CQ/QRZ HisCall   ? ]  49    5   Optional (not needed?)
     10    [CQ/QRZ HisCall   _ ]  68    5   Optional
     11    [CQ/QRZ HisCall Grid]  68    ?   Optional
*/

  float EbNodBEstimated;
  int err=0;
  int nSubmode=*ns0;
  float b90=*b0;
  int nFadingModel=*nf0;

  if(pqra64codec==NULL) pqra64codec = qra64_init(QRA_USERAP);
  err=qra64_apset(pqra64codec,*nc1,*nc2,*ng2,*APtype);
  if(err<0) printf("ERROR: qra64_apset returned %d\n",err);

  if(*iset==0) {
    *rc = qra64_decode_fastfading(pqra64codec,&EbNodBEstimated,xdec,r,
    				  nSubmode,b90,nFadingModel);
    *snr = EbNodBEstimated - 31.0;
  }
}


code:
$ qra64_nico -h

QRA64 Mode Tests
2016, Nico Palermo - IV3NWV

---------------------------

Syntax: qra64 [-s<snrdb>] [-c<channel>] [-a<ap-type>] [-t<testtype>] [-h]
Options:
       -s<snrdb>   : set simulation SNR in 2500 Hz BW (default:-27.5 dB)
       -c<channel> : set channel type 0=AWGN (default) 1=Rayleigh
       -a<ap-type> : set decode type 0=NOAP 1=AUTOAP (default) 2=USERAP
       -t<testtype>: 0=simulate seq of msgs between IV3NWV and K1JT (default)
                     1=simulate K1JT receiving K1JT IV3NWV JN66
                     2=simulate fast-fading routines (option -c ignored)
Options used only for fast-fading simulations:
       -b          : 90% fading bandwidth in Hz [1..230 Hz] (default = 2.5 Hz)
       -m          : fading model. 0=Gauss, 1=Lorentz (default = Lorentz)
       -q          : qra64 submode. 0=QRA64A,... 4=QRA64E (default = QRA64A)
       -h: this help

#############################################################################
Usage example:

qra64 -c2 -t2 -a2 -b50 -m1 -q2 -s-26.0 -n50000

Simulate a fast-fading channel (-c2)
Evaluate decoder performance (-t2) with
USER_AP (-a2)
B90 = 50 Hz (-b50)
Lorentz model (-m1)
submode QRA64C (-q2)
Input SNR = -26.0 dB (-s-26.0)
Simulate 50000 transmissions (-n50000)

(type qra64 -h for command syntax)

Command Output:

Input SNR = -26.0dB ap-mode=USER AP

Simulating the fast-fading channel
B90=50.00 Hz - Fading Model=Lorentz - Submode=QRA64C
Decoder metric = Matched to fast-fading signal

Encoding msg [K1JT IV3NWV JN66]
50000 transmissions will be simulated

K1JT decoder enabled for [CQ    ?     ?/blank]
K1JT decoder enabled for [K1JT  ?     ?/blank]
K1JT decoder enabled for [?    IV3NWV ?/blank]
K1JT decoder enabled for [K1JT IV3NWV ?]
K1JT decoder enabled for [K1JT IV3NWV JN66]
K1JT decoder enabled for [CQ   IV3NWV ?/b/JN66]

Now decoding with K1JT's decoder...
    5.5 %
Undetected error with rc=6
   13.1 %
Undetected error with rc=7
   70.8 %
Undetected error with rc=1
   75.8 %
Undetected error with rc=9
   88.9 %
Undetected error with rc=6
  100.0 %

Msgs transmitted:50000
Msg decoded:

rc= 0     0 with [?    ?    ?] AP0
rc= 1     0 with [CQ   ?    ?] AP27
rc= 2     0 with [CQ   ?     ] AP42
rc= 3   145 with [CALL ?    ?] AP29
rc= 4     0 with [CALL ?     ] AP44
rc= 5   12085 with [CALL CALL ?] AP57
rc= 6     0 with [?    CALL ?] AP29
rc= 7     0 with [?    CALL  ] AP44
rc= 8   24307 with [CALL CALL G] AP72
rc= 9     0 with [CQ   CALL ?] AP55
rc=10     0 with [CQ   CALL  ] AP70
rc=11     0 with [CQ   CALL G] AP70

Total: 36537/50000 (5 undetected errors)

Estimated SNR (average in dB) = -26.26 dB


code:
program qra64sim

! Generate simulated QRA64 data for testing the decoder.

  use wavhdr
  use packjt
  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,arg*12
  character msgsent*22

  nargs=iargc()
  if(nargs.ne. 7) then
     print *, 'Usage:   qra64sim         "msg"     A-E Nsigs fDop DT Nfiles SNR'
     print *, 'Example  qra64sim "K1ABC W9XYZ EN37" A   10   0.2 0.0   1     0'
     go to 999
  endif
  call getarg(1,msg)
  call getarg(2,csubmode)
  mode64=2**(ichar(csubmode)-ichar('A'))
  call getarg(3,arg)
  read(arg,*) nsigs
  call getarg(4,arg)
  read(arg,*) fspread
  call getarg(5,arg)
  read(arg,*) xdt
  call getarg(6,arg)
  read(arg,*) nfiles
  call getarg(7,arg)
  read(arg,*) snrdb

  if(mode64.ge.8) nsigs=1
  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

  write(*,1000)
1000 format('File  Sig    Freq  A-E   S/N   DT   Dop    Message'/60('-'))

  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

     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)
        if(nsigs.eq.1) f0=1000.0
        xsnr=snrdb
        if(snrdb.eq.0.0) xsnr=-20 - 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,f6.1,1x,a22)

        phi=0.d0
        dphi=0.d0
        k=(xdt+1.0)*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*mode64
              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.
        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(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


code:
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:	 libqra64.a qracodes.exe

OBJS1 = normrnd.o npfwht.o pdmath.o qra12_63_64_irr_b.o \
	qra13_64_64_irr_e.o qracodes.o

libqra64.a: $(OBJS1)
	ar cr libqra64.a $(OBJS1)
	ranlib libqra64.a

OBJS2 = main.o

qracodes.exe: $(OBJS2)
	${CC} -o qracodes.exe $(OBJS2) libqra64.a -lm

.PHONY : clean

clean:
	$(RM) *.o libqra64.a qracodes.exe


code:
# 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
0.6 10000
1.1 10000
1.6 10000
2.1 10000


code:
# 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


code:
# 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