From a8200fdfdd7d92d31a45a52affb5d610812cbee5 Mon Sep 17 00:00:00 2001 From: Joe Taylor <joe@princeton.edu> Date: Tue, 23 Feb 2021 09:20:46 -0500 Subject: [PATCH] Rempove remaining QRA64 code and direcories. --- lib/qra/qra64/Makefile.Win | 30 - lib/qra/qra64/fadengauss.c | 302 --------- lib/qra/qra64/fadenlorentz.c | 304 --------- lib/qra/qra64/main.c | 746 --------------------- lib/qra/qra64/qra64.c | 1121 -------------------------------- lib/qra/qra64/qra64.h | 269 -------- lib/qra/qra64/qra64_all.c | 1050 ------------------------------ lib/qra/qra64/qra64_subs.c | 65 -- lib/qra/qra64/qra64example.txt | 88 --- lib/qra/qra64/qra64sim.f90 | 201 ------ lib/qra_loops.f90 | 137 ---- 11 files changed, 4313 deletions(-) delete mode 100644 lib/qra/qra64/Makefile.Win delete mode 100644 lib/qra/qra64/fadengauss.c delete mode 100644 lib/qra/qra64/fadenlorentz.c delete mode 100644 lib/qra/qra64/main.c delete mode 100644 lib/qra/qra64/qra64.c delete mode 100644 lib/qra/qra64/qra64.h delete mode 100644 lib/qra/qra64/qra64_all.c delete mode 100644 lib/qra/qra64/qra64_subs.c delete mode 100644 lib/qra/qra64/qra64example.txt delete mode 100644 lib/qra/qra64/qra64sim.f90 delete mode 100644 lib/qra_loops.f90 diff --git a/lib/qra/qra64/Makefile.Win b/lib/qra/qra64/Makefile.Win deleted file mode 100644 index 7bd10c2d1..000000000 --- a/lib/qra/qra64/Makefile.Win +++ /dev/null @@ -1,30 +0,0 @@ -FC = gfortran -CC = gcc -CFLAGS = -O2 -Wall -I. -D_WIN32 - -# Default rules -%.o: %.c - ${CC} ${CFLAGS} -c $< -%.o: %.f - ${FC} ${FFLAGS} -c $< -%.o: %.F - ${FC} ${FFLAGS} -c $< -%.o: %.f90 - ${FC} ${FFLAGS} -c $< -%.o: %.F90 - ${FC} ${FFLAGS} -c $< - -all: qra64.exe - -OBJS1 = main.o qra64.o -qra64.exe: $(OBJS1) - ${CC} -o qra64.exe $(OBJS1) ../qracodes/libqra64.a -lm - -OBJS2 = qra64sim.o options.o wavhdr.o -qra64sim.exe: $(OBJS2) - ${FC} -o qra64sim.exe $(OBJS2) ../qracodes/libqra64.a -lm - -.PHONY : clean - -clean: - $(RM) *.o qra64.exe qra64sim.exe diff --git a/lib/qra/qra64/fadengauss.c b/lib/qra/qra64/fadengauss.c deleted file mode 100644 index c229626f4..000000000 --- a/lib/qra/qra64/fadengauss.c +++ /dev/null @@ -1,302 +0,0 @@ -// Gaussian energy fading tables for QRA64 -static const int glen_tab_gauss[64] = { - 2, 2, 2, 2, 2, 2, 2, 2, - 2, 2, 2, 2, 2, 2, 2, 2, - 3, 3, 3, 3, 3, 3, 3, 3, - 4, 4, 4, 4, 5, 5, 5, 6, - 6, 6, 7, 7, 8, 8, 9, 10, - 10, 11, 12, 13, 14, 15, 17, 18, - 19, 21, 23, 25, 27, 29, 32, 34, - 37, 41, 44, 48, 52, 57, 62, 65 -}; -static const float ggauss1[2] = { -0.0296f, 0.9101f -}; -static const float ggauss2[2] = { -0.0350f, 0.8954f -}; -static const float ggauss3[2] = { -0.0411f, 0.8787f -}; -static const float ggauss4[2] = { -0.0483f, 0.8598f -}; -static const float ggauss5[2] = { -0.0566f, 0.8387f -}; -static const float ggauss6[2] = { -0.0660f, 0.8154f -}; -static const float ggauss7[2] = { -0.0767f, 0.7898f -}; -static const float ggauss8[2] = { -0.0886f, 0.7621f -}; -static const float ggauss9[2] = { -0.1017f, 0.7325f -}; -static const float ggauss10[2] = { -0.1159f, 0.7012f -}; -static const float ggauss11[2] = { -0.1310f, 0.6687f -}; -static const float ggauss12[2] = { -0.1465f, 0.6352f -}; -static const float ggauss13[2] = { -0.1621f, 0.6013f -}; -static const float ggauss14[2] = { -0.1771f, 0.5674f -}; -static const float ggauss15[2] = { -0.1911f, 0.5339f -}; -static const float ggauss16[2] = { -0.2034f, 0.5010f -}; -static const float ggauss17[3] = { -0.0299f, 0.2135f, 0.4690f -}; -static const float ggauss18[3] = { -0.0369f, 0.2212f, 0.4383f -}; -static const float ggauss19[3] = { -0.0454f, 0.2263f, 0.4088f -}; -static const float ggauss20[3] = { -0.0552f, 0.2286f, 0.3806f -}; -static const float ggauss21[3] = { -0.0658f, 0.2284f, 0.3539f -}; -static const float ggauss22[3] = { -0.0766f, 0.2258f, 0.3287f -}; -static const float ggauss23[3] = { -0.0869f, 0.2212f, 0.3049f -}; -static const float ggauss24[3] = { -0.0962f, 0.2148f, 0.2826f -}; -static const float ggauss25[4] = { -0.0351f, 0.1041f, 0.2071f, 0.2616f -}; -static const float ggauss26[4] = { -0.0429f, 0.1102f, 0.1984f, 0.2420f -}; -static const float ggauss27[4] = { -0.0508f, 0.1145f, 0.1890f, 0.2237f -}; -static const float ggauss28[4] = { -0.0582f, 0.1169f, 0.1791f, 0.2067f -}; -static const float ggauss29[5] = { -0.0289f, 0.0648f, 0.1176f, 0.1689f, 0.1908f -}; -static const float ggauss30[5] = { -0.0351f, 0.0703f, 0.1168f, 0.1588f, 0.1760f -}; -static const float ggauss31[5] = { -0.0411f, 0.0745f, 0.1146f, 0.1488f, 0.1623f -}; -static const float ggauss32[6] = { -0.0246f, 0.0466f, 0.0773f, 0.1115f, 0.1390f, 0.1497f -}; -static const float ggauss33[6] = { -0.0297f, 0.0512f, 0.0788f, 0.1075f, 0.1295f, 0.1379f -}; -static const float ggauss34[6] = { -0.0345f, 0.0549f, 0.0791f, 0.1029f, 0.1205f, 0.1270f -}; -static const float ggauss35[7] = { -0.0240f, 0.0387f, 0.0575f, 0.0784f, 0.0979f, 0.1118f, 0.1169f -}; -static const float ggauss36[7] = { -0.0281f, 0.0422f, 0.0590f, 0.0767f, 0.0926f, 0.1037f, 0.1076f -}; -static const float ggauss37[8] = { -0.0212f, 0.0318f, 0.0449f, 0.0596f, 0.0744f, 0.0872f, 0.0960f, 0.0991f -}; -static const float ggauss38[8] = { -0.0247f, 0.0348f, 0.0467f, 0.0593f, 0.0716f, 0.0819f, 0.0887f, 0.0911f -}; -static const float ggauss39[9] = { -0.0199f, 0.0278f, 0.0372f, 0.0476f, 0.0584f, 0.0684f, 0.0766f, 0.0819f, -0.0838f -}; -static const float ggauss40[10] = { -0.0166f, 0.0228f, 0.0303f, 0.0388f, 0.0478f, 0.0568f, 0.0649f, 0.0714f, -0.0756f, 0.0771f -}; -static const float ggauss41[10] = { -0.0193f, 0.0254f, 0.0322f, 0.0397f, 0.0474f, 0.0548f, 0.0613f, 0.0664f, -0.0697f, 0.0709f -}; -static const float ggauss42[11] = { -0.0168f, 0.0217f, 0.0273f, 0.0335f, 0.0399f, 0.0464f, 0.0524f, 0.0576f, -0.0617f, 0.0643f, 0.0651f -}; -static const float ggauss43[12] = { -0.0151f, 0.0191f, 0.0237f, 0.0288f, 0.0342f, 0.0396f, 0.0449f, 0.0498f, -0.0540f, 0.0572f, 0.0592f, 0.0599f -}; -static const float ggauss44[13] = { -0.0138f, 0.0171f, 0.0210f, 0.0252f, 0.0297f, 0.0343f, 0.0388f, 0.0432f, -0.0471f, 0.0504f, 0.0529f, 0.0545f, 0.0550f -}; -static const float ggauss45[14] = { -0.0128f, 0.0157f, 0.0189f, 0.0224f, 0.0261f, 0.0300f, 0.0339f, 0.0377f, -0.0412f, 0.0444f, 0.0470f, 0.0489f, 0.0501f, 0.0505f -}; -static const float ggauss46[15] = { -0.0121f, 0.0146f, 0.0173f, 0.0202f, 0.0234f, 0.0266f, 0.0299f, 0.0332f, -0.0363f, 0.0391f, 0.0416f, 0.0437f, 0.0452f, 0.0461f, 0.0464f -}; -static const float ggauss47[17] = { -0.0097f, 0.0116f, 0.0138f, 0.0161f, 0.0186f, 0.0212f, 0.0239f, 0.0267f, -0.0294f, 0.0321f, 0.0346f, 0.0369f, 0.0389f, 0.0405f, 0.0417f, 0.0424f, -0.0427f -}; -static const float ggauss48[18] = { -0.0096f, 0.0113f, 0.0131f, 0.0151f, 0.0172f, 0.0194f, 0.0217f, 0.0241f, -0.0264f, 0.0287f, 0.0308f, 0.0329f, 0.0347f, 0.0362f, 0.0375f, 0.0384f, -0.0390f, 0.0392f -}; -static const float ggauss49[19] = { -0.0095f, 0.0110f, 0.0126f, 0.0143f, 0.0161f, 0.0180f, 0.0199f, 0.0219f, -0.0239f, 0.0258f, 0.0277f, 0.0294f, 0.0310f, 0.0325f, 0.0337f, 0.0347f, -0.0354f, 0.0358f, 0.0360f -}; -static const float ggauss50[21] = { -0.0083f, 0.0095f, 0.0108f, 0.0122f, 0.0136f, 0.0152f, 0.0168f, 0.0184f, -0.0201f, 0.0217f, 0.0234f, 0.0250f, 0.0265f, 0.0279f, 0.0292f, 0.0303f, -0.0313f, 0.0320f, 0.0326f, 0.0329f, 0.0330f -}; -static const float ggauss51[23] = { -0.0074f, 0.0084f, 0.0095f, 0.0106f, 0.0118f, 0.0131f, 0.0144f, 0.0157f, -0.0171f, 0.0185f, 0.0199f, 0.0213f, 0.0227f, 0.0240f, 0.0252f, 0.0263f, -0.0273f, 0.0282f, 0.0290f, 0.0296f, 0.0300f, 0.0303f, 0.0303f -}; -static const float ggauss52[25] = { -0.0068f, 0.0076f, 0.0085f, 0.0094f, 0.0104f, 0.0115f, 0.0126f, 0.0137f, -0.0149f, 0.0160f, 0.0172f, 0.0184f, 0.0196f, 0.0207f, 0.0218f, 0.0228f, -0.0238f, 0.0247f, 0.0255f, 0.0262f, 0.0268f, 0.0273f, 0.0276f, 0.0278f, -0.0279f -}; -static const float ggauss53[27] = { -0.0063f, 0.0070f, 0.0078f, 0.0086f, 0.0094f, 0.0103f, 0.0112f, 0.0121f, -0.0131f, 0.0141f, 0.0151f, 0.0161f, 0.0170f, 0.0180f, 0.0190f, 0.0199f, -0.0208f, 0.0216f, 0.0224f, 0.0231f, 0.0237f, 0.0243f, 0.0247f, 0.0251f, -0.0254f, 0.0255f, 0.0256f -}; -static const float ggauss54[29] = { -0.0060f, 0.0066f, 0.0072f, 0.0079f, 0.0086f, 0.0093f, 0.0101f, 0.0109f, -0.0117f, 0.0125f, 0.0133f, 0.0142f, 0.0150f, 0.0159f, 0.0167f, 0.0175f, -0.0183f, 0.0190f, 0.0197f, 0.0204f, 0.0210f, 0.0216f, 0.0221f, 0.0225f, -0.0228f, 0.0231f, 0.0233f, 0.0234f, 0.0235f -}; -static const float ggauss55[32] = { -0.0053f, 0.0058f, 0.0063f, 0.0068f, 0.0074f, 0.0080f, 0.0086f, 0.0093f, -0.0099f, 0.0106f, 0.0113f, 0.0120f, 0.0127f, 0.0134f, 0.0141f, 0.0148f, -0.0155f, 0.0162f, 0.0168f, 0.0174f, 0.0180f, 0.0186f, 0.0191f, 0.0196f, -0.0201f, 0.0204f, 0.0208f, 0.0211f, 0.0213f, 0.0214f, 0.0215f, 0.0216f -}; -static const float ggauss56[34] = { -0.0052f, 0.0056f, 0.0060f, 0.0065f, 0.0070f, 0.0075f, 0.0080f, 0.0086f, -0.0091f, 0.0097f, 0.0103f, 0.0109f, 0.0115f, 0.0121f, 0.0127f, 0.0133f, -0.0138f, 0.0144f, 0.0150f, 0.0155f, 0.0161f, 0.0166f, 0.0170f, 0.0175f, -0.0179f, 0.0183f, 0.0186f, 0.0189f, 0.0192f, 0.0194f, 0.0196f, 0.0197f, -0.0198f, 0.0198f -}; -static const float ggauss57[37] = { -0.0047f, 0.0051f, 0.0055f, 0.0058f, 0.0063f, 0.0067f, 0.0071f, 0.0076f, -0.0080f, 0.0085f, 0.0090f, 0.0095f, 0.0100f, 0.0105f, 0.0110f, 0.0115f, -0.0120f, 0.0125f, 0.0130f, 0.0134f, 0.0139f, 0.0144f, 0.0148f, 0.0152f, -0.0156f, 0.0160f, 0.0164f, 0.0167f, 0.0170f, 0.0173f, 0.0175f, 0.0177f, -0.0179f, 0.0180f, 0.0181f, 0.0181f, 0.0182f -}; -static const float ggauss58[41] = { -0.0041f, 0.0044f, 0.0047f, 0.0050f, 0.0054f, 0.0057f, 0.0060f, 0.0064f, -0.0068f, 0.0072f, 0.0076f, 0.0080f, 0.0084f, 0.0088f, 0.0092f, 0.0096f, -0.0101f, 0.0105f, 0.0109f, 0.0113f, 0.0117f, 0.0121f, 0.0125f, 0.0129f, -0.0133f, 0.0137f, 0.0140f, 0.0144f, 0.0147f, 0.0150f, 0.0153f, 0.0155f, -0.0158f, 0.0160f, 0.0162f, 0.0163f, 0.0164f, 0.0165f, 0.0166f, 0.0167f, -0.0167f -}; -static const float ggauss59[44] = { -0.0039f, 0.0042f, 0.0044f, 0.0047f, 0.0050f, 0.0053f, 0.0056f, 0.0059f, -0.0062f, 0.0065f, 0.0068f, 0.0072f, 0.0075f, 0.0079f, 0.0082f, 0.0086f, -0.0089f, 0.0093f, 0.0096f, 0.0100f, 0.0104f, 0.0107f, 0.0110f, 0.0114f, -0.0117f, 0.0120f, 0.0124f, 0.0127f, 0.0130f, 0.0132f, 0.0135f, 0.0138f, -0.0140f, 0.0142f, 0.0144f, 0.0146f, 0.0148f, 0.0149f, 0.0150f, 0.0151f, -0.0152f, 0.0153f, 0.0153f, 0.0153f -}; -static const float ggauss60[48] = { -0.0036f, 0.0038f, 0.0040f, 0.0042f, 0.0044f, 0.0047f, 0.0049f, 0.0052f, -0.0055f, 0.0057f, 0.0060f, 0.0063f, 0.0066f, 0.0068f, 0.0071f, 0.0074f, -0.0077f, 0.0080f, 0.0083f, 0.0086f, 0.0089f, 0.0092f, 0.0095f, 0.0098f, -0.0101f, 0.0104f, 0.0107f, 0.0109f, 0.0112f, 0.0115f, 0.0117f, 0.0120f, -0.0122f, 0.0124f, 0.0126f, 0.0128f, 0.0130f, 0.0132f, 0.0134f, 0.0135f, -0.0136f, 0.0137f, 0.0138f, 0.0139f, 0.0140f, 0.0140f, 0.0140f, 0.0140f -}; -static const float ggauss61[52] = { -0.0033f, 0.0035f, 0.0037f, 0.0039f, 0.0041f, 0.0043f, 0.0045f, 0.0047f, -0.0049f, 0.0051f, 0.0053f, 0.0056f, 0.0058f, 0.0060f, 0.0063f, 0.0065f, -0.0068f, 0.0070f, 0.0073f, 0.0075f, 0.0078f, 0.0080f, 0.0083f, 0.0085f, -0.0088f, 0.0090f, 0.0093f, 0.0095f, 0.0098f, 0.0100f, 0.0102f, 0.0105f, -0.0107f, 0.0109f, 0.0111f, 0.0113f, 0.0115f, 0.0116f, 0.0118f, 0.0120f, -0.0121f, 0.0122f, 0.0124f, 0.0125f, 0.0126f, 0.0126f, 0.0127f, 0.0128f, -0.0128f, 0.0129f, 0.0129f, 0.0129f -}; -static const float ggauss62[57] = { -0.0030f, 0.0031f, 0.0033f, 0.0034f, 0.0036f, 0.0038f, 0.0039f, 0.0041f, -0.0043f, 0.0045f, 0.0047f, 0.0048f, 0.0050f, 0.0052f, 0.0054f, 0.0056f, -0.0058f, 0.0060f, 0.0063f, 0.0065f, 0.0067f, 0.0069f, 0.0071f, 0.0073f, -0.0075f, 0.0077f, 0.0080f, 0.0082f, 0.0084f, 0.0086f, 0.0088f, 0.0090f, -0.0092f, 0.0094f, 0.0096f, 0.0097f, 0.0099f, 0.0101f, 0.0103f, 0.0104f, -0.0106f, 0.0107f, 0.0108f, 0.0110f, 0.0111f, 0.0112f, 0.0113f, 0.0114f, -0.0115f, 0.0116f, 0.0116f, 0.0117f, 0.0117f, 0.0118f, 0.0118f, 0.0118f, -0.0118f -}; -static const float ggauss63[62] = { -0.0027f, 0.0029f, 0.0030f, 0.0031f, 0.0032f, 0.0034f, 0.0035f, 0.0037f, -0.0038f, 0.0040f, 0.0041f, 0.0043f, 0.0045f, 0.0046f, 0.0048f, 0.0049f, -0.0051f, 0.0053f, 0.0055f, 0.0056f, 0.0058f, 0.0060f, 0.0062f, 0.0063f, -0.0065f, 0.0067f, 0.0069f, 0.0071f, 0.0072f, 0.0074f, 0.0076f, 0.0078f, -0.0079f, 0.0081f, 0.0083f, 0.0084f, 0.0086f, 0.0088f, 0.0089f, 0.0091f, -0.0092f, 0.0094f, 0.0095f, 0.0096f, 0.0098f, 0.0099f, 0.0100f, 0.0101f, -0.0102f, 0.0103f, 0.0104f, 0.0105f, 0.0105f, 0.0106f, 0.0107f, 0.0107f, -0.0108f, 0.0108f, 0.0108f, 0.0108f, 0.0109f, 0.0109f -}; -static const float ggauss64[65] = { -0.0028f, 0.0029f, 0.0030f, 0.0031f, 0.0032f, 0.0034f, 0.0035f, 0.0036f, -0.0037f, 0.0039f, 0.0040f, 0.0041f, 0.0043f, 0.0044f, 0.0046f, 0.0047f, -0.0048f, 0.0050f, 0.0051f, 0.0053f, 0.0054f, 0.0056f, 0.0057f, 0.0059f, -0.0060f, 0.0062f, 0.0063f, 0.0065f, 0.0066f, 0.0068f, 0.0069f, 0.0071f, -0.0072f, 0.0074f, 0.0075f, 0.0077f, 0.0078f, 0.0079f, 0.0081f, 0.0082f, -0.0083f, 0.0084f, 0.0086f, 0.0087f, 0.0088f, 0.0089f, 0.0090f, 0.0091f, -0.0092f, 0.0093f, 0.0094f, 0.0094f, 0.0095f, 0.0096f, 0.0097f, 0.0097f, -0.0098f, 0.0098f, 0.0099f, 0.0099f, 0.0099f, 0.0099f, 0.0100f, 0.0100f, -0.0100f -}; -static const float *gptr_tab_gauss[64] = { -ggauss1, ggauss2, ggauss3, ggauss4, -ggauss5, ggauss6, ggauss7, ggauss8, -ggauss9, ggauss10, ggauss11, ggauss12, -ggauss13, ggauss14, ggauss15, ggauss16, -ggauss17, ggauss18, ggauss19, ggauss20, -ggauss21, ggauss22, ggauss23, ggauss24, -ggauss25, ggauss26, ggauss27, ggauss28, -ggauss29, ggauss30, ggauss31, ggauss32, -ggauss33, ggauss34, ggauss35, ggauss36, -ggauss37, ggauss38, ggauss39, ggauss40, -ggauss41, ggauss42, ggauss43, ggauss44, -ggauss45, ggauss46, ggauss47, ggauss48, -ggauss49, ggauss50, ggauss51, ggauss52, -ggauss53, ggauss54, ggauss55, ggauss56, -ggauss57, ggauss58, ggauss59, ggauss60, -ggauss61, ggauss62, ggauss63, ggauss64 -}; diff --git a/lib/qra/qra64/fadenlorentz.c b/lib/qra/qra64/fadenlorentz.c deleted file mode 100644 index 673a89033..000000000 --- a/lib/qra/qra64/fadenlorentz.c +++ /dev/null @@ -1,304 +0,0 @@ -// Lorentz energy fading tables for QRA64 -static const int glen_tab_lorentz[64] = { - 2, 2, 2, 2, 2, 2, 2, 2, - 2, 2, 2, 2, 2, 2, 3, 3, - 3, 3, 3, 3, 3, 4, 4, 4, - 4, 4, 5, 5, 5, 5, 6, 6, - 7, 7, 7, 8, 8, 9, 10, 10, - 11, 12, 13, 14, 15, 16, 17, 19, - 20, 22, 23, 25, 27, 30, 32, 35, - 38, 41, 45, 49, 53, 57, 62, 65 -}; -static const float glorentz1[2] = { -0.0214f, 0.9107f -}; -static const float glorentz2[2] = { -0.0244f, 0.9030f -}; -static const float glorentz3[2] = { -0.0280f, 0.8950f -}; -static const float glorentz4[2] = { -0.0314f, 0.8865f -}; -static const float glorentz5[2] = { -0.0349f, 0.8773f -}; -static const float glorentz6[2] = { -0.0388f, 0.8675f -}; -static const float glorentz7[2] = { -0.0426f, 0.8571f -}; -static const float glorentz8[2] = { -0.0463f, 0.8459f -}; -static const float glorentz9[2] = { -0.0500f, 0.8339f -}; -static const float glorentz10[2] = { -0.0538f, 0.8210f -}; -static const float glorentz11[2] = { -0.0579f, 0.8074f -}; -static const float glorentz12[2] = { -0.0622f, 0.7930f -}; -static const float glorentz13[2] = { -0.0668f, 0.7777f -}; -static const float glorentz14[2] = { -0.0715f, 0.7616f -}; -static const float glorentz15[3] = { -0.0196f, 0.0765f, 0.7445f -}; -static const float glorentz16[3] = { -0.0210f, 0.0816f, 0.7267f -}; -static const float glorentz17[3] = { -0.0226f, 0.0870f, 0.7080f -}; -static const float glorentz18[3] = { -0.0242f, 0.0925f, 0.6885f -}; -static const float glorentz19[3] = { -0.0259f, 0.0981f, 0.6682f -}; -static const float glorentz20[3] = { -0.0277f, 0.1039f, 0.6472f -}; -static const float glorentz21[3] = { -0.0296f, 0.1097f, 0.6255f -}; -static const float glorentz22[4] = { -0.0143f, 0.0316f, 0.1155f, 0.6031f -}; -static const float glorentz23[4] = { -0.0153f, 0.0337f, 0.1213f, 0.5803f -}; -static const float glorentz24[4] = { -0.0163f, 0.0358f, 0.1270f, 0.5570f -}; -static const float glorentz25[4] = { -0.0174f, 0.0381f, 0.1325f, 0.5333f -}; -static const float glorentz26[4] = { -0.0186f, 0.0405f, 0.1378f, 0.5095f -}; -static const float glorentz27[5] = { -0.0113f, 0.0198f, 0.0429f, 0.1428f, 0.4855f -}; -static const float glorentz28[5] = { -0.0120f, 0.0211f, 0.0455f, 0.1473f, 0.4615f -}; -static const float glorentz29[5] = { -0.0129f, 0.0225f, 0.0481f, 0.1514f, 0.4376f -}; -static const float glorentz30[5] = { -0.0137f, 0.0239f, 0.0508f, 0.1549f, 0.4140f -}; -static const float glorentz31[6] = { -0.0095f, 0.0147f, 0.0254f, 0.0536f, 0.1578f, 0.3907f -}; -static const float glorentz32[6] = { -0.0101f, 0.0156f, 0.0270f, 0.0564f, 0.1600f, 0.3680f -}; -static const float glorentz33[7] = { -0.0076f, 0.0109f, 0.0167f, 0.0287f, 0.0592f, 0.1614f, 0.3458f -}; -static const float glorentz34[7] = { -0.0081f, 0.0116f, 0.0178f, 0.0305f, 0.0621f, 0.1620f, 0.3243f -}; -static const float glorentz35[7] = { -0.0087f, 0.0124f, 0.0190f, 0.0324f, 0.0649f, 0.1618f, 0.3035f -}; -static const float glorentz36[8] = { -0.0069f, 0.0093f, 0.0133f, 0.0203f, 0.0343f, 0.0676f, 0.1607f, 0.2836f -}; -static const float glorentz37[8] = { -0.0074f, 0.0100f, 0.0142f, 0.0216f, 0.0362f, 0.0702f, 0.1588f, 0.2645f -}; -static const float glorentz38[9] = { -0.0061f, 0.0080f, 0.0107f, 0.0152f, 0.0230f, 0.0382f, 0.0726f, 0.1561f, -0.2464f -}; -static const float glorentz39[10] = { -0.0052f, 0.0066f, 0.0086f, 0.0115f, 0.0162f, 0.0244f, 0.0402f, 0.0747f, -0.1526f, 0.2291f -}; -static const float glorentz40[10] = { -0.0056f, 0.0071f, 0.0092f, 0.0123f, 0.0173f, 0.0259f, 0.0422f, 0.0766f, -0.1484f, 0.2128f -}; -static const float glorentz41[11] = { -0.0049f, 0.0061f, 0.0076f, 0.0098f, 0.0132f, 0.0184f, 0.0274f, 0.0441f, -0.0780f, 0.1437f, 0.1975f -}; -static const float glorentz42[12] = { -0.0044f, 0.0053f, 0.0065f, 0.0082f, 0.0106f, 0.0141f, 0.0196f, 0.0290f, -0.0460f, 0.0791f, 0.1384f, 0.1831f -}; -static const float glorentz43[13] = { -0.0040f, 0.0048f, 0.0057f, 0.0070f, 0.0088f, 0.0113f, 0.0150f, 0.0209f, -0.0305f, 0.0477f, 0.0797f, 0.1327f, 0.1695f -}; -static const float glorentz44[14] = { -0.0037f, 0.0043f, 0.0051f, 0.0062f, 0.0075f, 0.0094f, 0.0121f, 0.0160f, -0.0221f, 0.0321f, 0.0493f, 0.0799f, 0.1267f, 0.1568f -}; -static const float glorentz45[15] = { -0.0035f, 0.0040f, 0.0047f, 0.0055f, 0.0066f, 0.0081f, 0.0101f, 0.0129f, -0.0171f, 0.0234f, 0.0335f, 0.0506f, 0.0795f, 0.1204f, 0.1450f -}; -static const float glorentz46[16] = { -0.0033f, 0.0037f, 0.0043f, 0.0050f, 0.0059f, 0.0071f, 0.0087f, 0.0108f, -0.0138f, 0.0181f, 0.0246f, 0.0349f, 0.0517f, 0.0786f, 0.1141f, 0.1340f -}; -static const float glorentz47[17] = { -0.0031f, 0.0035f, 0.0040f, 0.0046f, 0.0054f, 0.0064f, 0.0077f, 0.0093f, -0.0116f, 0.0147f, 0.0192f, 0.0259f, 0.0362f, 0.0525f, 0.0773f, 0.1076f, -0.1237f -}; -static const float glorentz48[19] = { -0.0027f, 0.0030f, 0.0034f, 0.0038f, 0.0043f, 0.0050f, 0.0058f, 0.0069f, -0.0082f, 0.0100f, 0.0123f, 0.0156f, 0.0203f, 0.0271f, 0.0374f, 0.0530f, -0.0755f, 0.1013f, 0.1141f -}; -static const float glorentz49[20] = { -0.0026f, 0.0029f, 0.0032f, 0.0036f, 0.0041f, 0.0047f, 0.0054f, 0.0063f, -0.0074f, 0.0088f, 0.0107f, 0.0131f, 0.0165f, 0.0213f, 0.0282f, 0.0383f, -0.0531f, 0.0734f, 0.0950f, 0.1053f -}; -static const float glorentz50[22] = { -0.0023f, 0.0025f, 0.0028f, 0.0031f, 0.0035f, 0.0039f, 0.0044f, 0.0050f, -0.0058f, 0.0067f, 0.0079f, 0.0094f, 0.0114f, 0.0139f, 0.0175f, 0.0223f, -0.0292f, 0.0391f, 0.0529f, 0.0709f, 0.0889f, 0.0971f -}; -static const float glorentz51[23] = { -0.0023f, 0.0025f, 0.0027f, 0.0030f, 0.0034f, 0.0037f, 0.0042f, 0.0048f, -0.0054f, 0.0062f, 0.0072f, 0.0085f, 0.0100f, 0.0121f, 0.0148f, 0.0184f, -0.0233f, 0.0301f, 0.0396f, 0.0524f, 0.0681f, 0.0829f, 0.0894f -}; -static const float glorentz52[25] = { -0.0021f, 0.0023f, 0.0025f, 0.0027f, 0.0030f, 0.0033f, 0.0036f, 0.0040f, -0.0045f, 0.0051f, 0.0058f, 0.0067f, 0.0077f, 0.0090f, 0.0107f, 0.0128f, -0.0156f, 0.0192f, 0.0242f, 0.0308f, 0.0398f, 0.0515f, 0.0650f, 0.0772f, -0.0824f -}; -static const float glorentz53[27] = { -0.0019f, 0.0021f, 0.0022f, 0.0024f, 0.0027f, 0.0029f, 0.0032f, 0.0035f, -0.0039f, 0.0044f, 0.0049f, 0.0055f, 0.0062f, 0.0072f, 0.0083f, 0.0096f, -0.0113f, 0.0135f, 0.0164f, 0.0201f, 0.0249f, 0.0314f, 0.0398f, 0.0502f, -0.0619f, 0.0718f, 0.0759f -}; -static const float glorentz54[30] = { -0.0017f, 0.0018f, 0.0019f, 0.0021f, 0.0022f, 0.0024f, 0.0026f, 0.0029f, -0.0031f, 0.0034f, 0.0038f, 0.0042f, 0.0047f, 0.0052f, 0.0059f, 0.0067f, -0.0076f, 0.0088f, 0.0102f, 0.0120f, 0.0143f, 0.0171f, 0.0208f, 0.0256f, -0.0317f, 0.0395f, 0.0488f, 0.0586f, 0.0666f, 0.0698f -}; -static const float glorentz55[32] = { -0.0016f, 0.0017f, 0.0018f, 0.0019f, 0.0021f, 0.0022f, 0.0024f, 0.0026f, -0.0028f, 0.0031f, 0.0034f, 0.0037f, 0.0041f, 0.0045f, 0.0050f, 0.0056f, -0.0063f, 0.0071f, 0.0081f, 0.0094f, 0.0108f, 0.0127f, 0.0149f, 0.0178f, -0.0214f, 0.0261f, 0.0318f, 0.0389f, 0.0470f, 0.0553f, 0.0618f, 0.0643f -}; -static const float glorentz56[35] = { -0.0014f, 0.0015f, 0.0016f, 0.0017f, 0.0018f, 0.0020f, 0.0021f, 0.0023f, -0.0024f, 0.0026f, 0.0028f, 0.0031f, 0.0033f, 0.0036f, 0.0040f, 0.0044f, -0.0049f, 0.0054f, 0.0060f, 0.0067f, 0.0076f, 0.0087f, 0.0099f, 0.0114f, -0.0133f, 0.0156f, 0.0184f, 0.0220f, 0.0264f, 0.0318f, 0.0381f, 0.0451f, -0.0520f, 0.0572f, 0.0591f -}; -static const float glorentz57[38] = { -0.0013f, 0.0014f, 0.0015f, 0.0016f, 0.0017f, 0.0018f, 0.0019f, 0.0020f, -0.0021f, 0.0023f, 0.0024f, 0.0026f, 0.0028f, 0.0031f, 0.0033f, 0.0036f, -0.0039f, 0.0043f, 0.0047f, 0.0052f, 0.0058f, 0.0064f, 0.0072f, 0.0081f, -0.0092f, 0.0104f, 0.0120f, 0.0139f, 0.0162f, 0.0190f, 0.0224f, 0.0265f, -0.0315f, 0.0371f, 0.0431f, 0.0487f, 0.0529f, 0.0544f -}; -static const float glorentz58[41] = { -0.0012f, 0.0013f, 0.0014f, 0.0014f, 0.0015f, 0.0016f, 0.0017f, 0.0018f, -0.0019f, 0.0020f, 0.0022f, 0.0023f, 0.0025f, 0.0026f, 0.0028f, 0.0030f, -0.0033f, 0.0036f, 0.0039f, 0.0042f, 0.0046f, 0.0050f, 0.0056f, 0.0061f, -0.0068f, 0.0076f, 0.0086f, 0.0097f, 0.0110f, 0.0125f, 0.0144f, 0.0167f, -0.0194f, 0.0226f, 0.0265f, 0.0309f, 0.0359f, 0.0409f, 0.0455f, 0.0488f, -0.0500f -}; -static const float glorentz59[45] = { -0.0011f, 0.0012f, 0.0012f, 0.0013f, 0.0013f, 0.0014f, 0.0015f, 0.0016f, -0.0016f, 0.0017f, 0.0018f, 0.0019f, 0.0021f, 0.0022f, 0.0023f, 0.0025f, -0.0026f, 0.0028f, 0.0030f, 0.0033f, 0.0035f, 0.0038f, 0.0041f, 0.0045f, -0.0049f, 0.0054f, 0.0059f, 0.0065f, 0.0072f, 0.0081f, 0.0090f, 0.0102f, -0.0115f, 0.0130f, 0.0149f, 0.0171f, 0.0197f, 0.0227f, 0.0263f, 0.0302f, -0.0345f, 0.0387f, 0.0425f, 0.0451f, 0.0460f -}; -static const float glorentz60[49] = { -0.0010f, 0.0011f, 0.0011f, 0.0012f, 0.0012f, 0.0013f, 0.0013f, 0.0014f, -0.0014f, 0.0015f, 0.0016f, 0.0017f, 0.0018f, 0.0019f, 0.0020f, 0.0021f, -0.0022f, 0.0024f, 0.0025f, 0.0027f, 0.0028f, 0.0030f, 0.0033f, 0.0035f, -0.0038f, 0.0041f, 0.0044f, 0.0048f, 0.0052f, 0.0057f, 0.0063f, 0.0069f, -0.0077f, 0.0085f, 0.0095f, 0.0106f, 0.0119f, 0.0135f, 0.0153f, 0.0174f, -0.0199f, 0.0227f, 0.0259f, 0.0293f, 0.0330f, 0.0365f, 0.0395f, 0.0415f, -0.0423f -}; -static const float glorentz61[53] = { -0.0009f, 0.0010f, 0.0010f, 0.0011f, 0.0011f, 0.0011f, 0.0012f, 0.0012f, -0.0013f, 0.0014f, 0.0014f, 0.0015f, 0.0016f, 0.0016f, 0.0017f, 0.0018f, -0.0019f, 0.0020f, 0.0021f, 0.0023f, 0.0024f, 0.0025f, 0.0027f, 0.0029f, -0.0031f, 0.0033f, 0.0035f, 0.0038f, 0.0041f, 0.0044f, 0.0047f, 0.0051f, -0.0056f, 0.0061f, 0.0067f, 0.0073f, 0.0081f, 0.0089f, 0.0099f, 0.0110f, -0.0124f, 0.0139f, 0.0156f, 0.0176f, 0.0199f, 0.0225f, 0.0253f, 0.0283f, -0.0314f, 0.0343f, 0.0367f, 0.0383f, 0.0389f -}; -static const float glorentz62[57] = { -0.0009f, 0.0009f, 0.0009f, 0.0010f, 0.0010f, 0.0011f, 0.0011f, 0.0011f, -0.0012f, 0.0012f, 0.0013f, 0.0013f, 0.0014f, 0.0015f, 0.0015f, 0.0016f, -0.0017f, 0.0018f, 0.0019f, 0.0020f, 0.0021f, 0.0022f, 0.0023f, 0.0024f, -0.0026f, 0.0027f, 0.0029f, 0.0031f, 0.0033f, 0.0035f, 0.0038f, 0.0040f, -0.0043f, 0.0047f, 0.0050f, 0.0055f, 0.0059f, 0.0064f, 0.0070f, 0.0077f, -0.0085f, 0.0093f, 0.0103f, 0.0114f, 0.0127f, 0.0142f, 0.0158f, 0.0177f, -0.0198f, 0.0221f, 0.0246f, 0.0272f, 0.0297f, 0.0321f, 0.0340f, 0.0353f, -0.0357f -}; -static const float glorentz63[62] = { -0.0008f, 0.0008f, 0.0009f, 0.0009f, 0.0009f, 0.0010f, 0.0010f, 0.0010f, -0.0011f, 0.0011f, 0.0011f, 0.0012f, 0.0012f, 0.0013f, 0.0013f, 0.0014f, -0.0015f, 0.0015f, 0.0016f, 0.0017f, 0.0017f, 0.0018f, 0.0019f, 0.0020f, -0.0021f, 0.0022f, 0.0023f, 0.0025f, 0.0026f, 0.0028f, 0.0029f, 0.0031f, -0.0033f, 0.0035f, 0.0038f, 0.0040f, 0.0043f, 0.0046f, 0.0050f, 0.0053f, -0.0058f, 0.0062f, 0.0068f, 0.0074f, 0.0081f, 0.0088f, 0.0097f, 0.0106f, -0.0117f, 0.0130f, 0.0144f, 0.0159f, 0.0176f, 0.0195f, 0.0216f, 0.0237f, -0.0259f, 0.0280f, 0.0299f, 0.0315f, 0.0325f, 0.0328f -}; -static const float glorentz64[65] = { -0.0008f, 0.0008f, 0.0008f, 0.0009f, 0.0009f, 0.0009f, 0.0010f, 0.0010f, -0.0010f, 0.0011f, 0.0011f, 0.0012f, 0.0012f, 0.0012f, 0.0013f, 0.0013f, -0.0014f, 0.0014f, 0.0015f, 0.0016f, 0.0016f, 0.0017f, 0.0018f, 0.0019f, -0.0020f, 0.0021f, 0.0022f, 0.0023f, 0.0024f, 0.0025f, 0.0027f, 0.0028f, -0.0030f, 0.0031f, 0.0033f, 0.0035f, 0.0038f, 0.0040f, 0.0043f, 0.0046f, -0.0049f, 0.0052f, 0.0056f, 0.0061f, 0.0066f, 0.0071f, 0.0077f, 0.0084f, -0.0091f, 0.0100f, 0.0109f, 0.0120f, 0.0132f, 0.0145f, 0.0159f, 0.0175f, -0.0192f, 0.0209f, 0.0228f, 0.0246f, 0.0264f, 0.0279f, 0.0291f, 0.0299f, -0.0301f -}; -static const float *gptr_tab_lorentz[64] = { -glorentz1, glorentz2, glorentz3, glorentz4, -glorentz5, glorentz6, glorentz7, glorentz8, -glorentz9, glorentz10, glorentz11, glorentz12, -glorentz13, glorentz14, glorentz15, glorentz16, -glorentz17, glorentz18, glorentz19, glorentz20, -glorentz21, glorentz22, glorentz23, glorentz24, -glorentz25, glorentz26, glorentz27, glorentz28, -glorentz29, glorentz30, glorentz31, glorentz32, -glorentz33, glorentz34, glorentz35, glorentz36, -glorentz37, glorentz38, glorentz39, glorentz40, -glorentz41, glorentz42, glorentz43, glorentz44, -glorentz45, glorentz46, glorentz47, glorentz48, -glorentz49, glorentz50, glorentz51, glorentz52, -glorentz53, glorentz54, glorentz55, glorentz56, -glorentz57, glorentz58, glorentz59, glorentz60, -glorentz61, glorentz62, glorentz63, glorentz64 -}; diff --git a/lib/qra/qra64/main.c b/lib/qra/qra64/main.c deleted file mode 100644 index b685fd316..000000000 --- a/lib/qra/qra64/main.c +++ /dev/null @@ -1,746 +0,0 @@ -/* -main.c -QRA64 mode encode/decode tests - -(c) 2016 - Nico Palermo, IV3NWV - -Thanks to Andrea Montefusco IW0HDV for his help on adapting the sources -to OSs other than MS Windows - ------------------------------------------------------------------------------- -This file is part of the qracodes project, a Forward Error Control -encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes. - -Files in this package: - main.c - this file - qra64.c/.h - qra64 mode encode/decoding functions - - ../qracodes/normrnd.{c,h} - random gaussian number generator - ../qracodes/npfwht.{c,h} - Fast Walsh-Hadamard Transforms - ../qracodes/pdmath.{c,h} - Elementary math on probability distributions - ../qracodes/qra12_63_64_irr_b.{c,h} - Tables for a QRA(12,63) irregular RA - code over GF(64) - ../qracodes/qra13_64_64_irr_e.{c,h} - Tables for a QRA(13,64) irregular RA - code over GF(64) - ../qracodes/qracodes.{c,h} - QRA codes encoding/decoding functions - -------------------------------------------------------------------------------- - - qracodes is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - qracodes is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with qracodes source distribution. - If not, see <http://www.gnu.org/licenses/>. - ------------------------------------------------------------------------------ - -The code used by the QRA64 mode is the code: QRA13_64_64_IRR_E: K=13 -N=64 Q=64 irregular QRA code (defined in qra13_64_64_irr_e.{h,c}). - -This code has been designed to include a CRC as the 13th information -symbol and improve the code UER (Undetected Error Rate). The CRC -symbol is not sent along the channel (the codes are punctured) and the -resulting code is still a (12,63) code with an effective code rate of -R = 12/63. -*/ - -// OS dependent defines and includes ------------------------------------------ - -#if _WIN32 // note the underscore: without it, it's not msdn official! -// Windows (x64 and x86) -#include <windows.h> // required only for GetTickCount(...) -#include <process.h> // _beginthread -#endif - -#if __linux__ -#include <unistd.h> -#include <time.h> - -unsigned GetTickCount(void) { - struct timespec ts; - unsigned theTick = 0U; - clock_gettime( CLOCK_REALTIME, &ts ); - theTick = ts.tv_nsec / 1000000; - theTick += ts.tv_sec * 1000; - return theTick; -} -#endif - -#if __APPLE__ -#endif - -#include <stdlib.h> -#include <stdio.h> -#include <string.h> - -#include "qra64.h" -#include "../qracodes/normrnd.h" // gaussian numbers generator - -// ---------------------------------------------------------------------------- - -// channel types -#define CHANNEL_AWGN 0 -#define CHANNEL_RAYLEIGH 1 -#define CHANNEL_FASTFADE 2 - -#define JT65_SNR_EBNO_OFFSET 29.1f // with the synch used in JT65 -#define QRA64_SNR_EBNO_OFFSET 31.0f // with the costas array synch - -void printwordd(char *msg, int *x, int size) -{ - int k; - printf("\n%s ",msg); - for (k=0;k<size;k++) - printf("%2d ",x[k]); - printf("\n"); -} -void printwordh(char *msg, int *x, int size) -{ - int k; - printf("\n%s ",msg); - for (k=0;k<size;k++) - printf("%02hx ",x[k]); - printf("\n"); -} - -#define NSAMPLES (QRA64_N*QRA64_M) - -static float rp[NSAMPLES]; -static float rq[NSAMPLES]; -static float chp[NSAMPLES]; -static float chq[NSAMPLES]; -static float r[NSAMPLES]; - -float *mfskchannel(int *x, int channel_type, float EbNodB) -{ -/* -Simulate an MFSK channel, either AWGN or Rayleigh. - -x is a pointer to the transmitted codeword, an array of QRA64_N -integers in the range 0..63. - -Returns the received symbol energies (squared amplitudes) as an array of -(QRA64_M*QRA64_N) floats. The first QRA64_M entries of this array are -the energies of the first symbol in the codeword. The second QRA64_M -entries are those of the second symbol, and so on up to the last codeword -symbol. -*/ - const float No = 1.0f; // noise spectral density - const float sigma = (float)sqrt(No/2.0f); // std dev of noise I/Q components - const float sigmach = (float)sqrt(1/2.0f); // std dev of channel I/Q gains - const float R = 1.0f*QRA64_K/QRA64_N; - - float EbNo = (float)pow(10,EbNodB/10); - float EsNo = 1.0f*QRA64_m*R*EbNo; - float Es = EsNo*No; - float A = (float)sqrt(Es); - int k; - - normrnd_s(rp,NSAMPLES,0,sigma); - normrnd_s(rq,NSAMPLES,0,sigma); - - if(EbNodB>-15) - if (channel_type == CHANNEL_AWGN) - for (k=0;k<QRA64_N;k++) - rp[k*QRA64_M+x[k]]+=A; - else - if (channel_type == CHANNEL_RAYLEIGH) { - normrnd_s(chp,QRA64_N,0,sigmach); - normrnd_s(chq,QRA64_N,0,sigmach); - for (k=0;k<QRA64_N;k++) { - rp[k*QRA64_M+x[k]]+=A*chp[k]; - rq[k*QRA64_M+x[k]]+=A*chq[k]; - } - } - else { - return 0; // unknown channel type - } - - // compute the squares of the amplitudes of the received samples - for (k=0;k<NSAMPLES;k++) - r[k] = rp[k]*rp[k] + rq[k]*rq[k]; - - return r; -} - -// These defines are some packed fields as computed by JT65 -#define CALL_IV3NWV 0x7F85AE7 -#define CALL_K1JT 0xF70DDD7 -#define GRID_JN66 0x3AE4 // JN66 -#define GRID_73 0x7ED0 // 73 - -char decode_type[12][32] = { - "[? ? ?] AP0", - "[CQ ? ?] AP27", - "[CQ ? ] AP42", - "[CALL ? ?] AP29", - "[CALL ? ] AP44", - "[CALL CALL ?] AP57", - "[? CALL ?] AP29", - "[? CALL ] AP44", - "[CALL CALL G] AP72", - "[CQ CALL ?] AP55", - "[CQ CALL ] AP70", - "[CQ CALL G] AP70" -}; -char apmode_type[3][32] = { - "NO AP", - "AUTO AP", - "USER AP" -}; - -int test_proc_1(int channel_type, float EbNodB, int mode) -{ -/* -Here we simulate the following (dummy) QSO: - -1) CQ IV3NWV -2) IV3NWV K1JT -3) K1JT IV3NWV 73 -4) IV3NWV K1JT 73 - -No message repetition is attempted - -The QSO is counted as successfull if IV3NWV received the last message -When mode=QRA_AUTOAP each decoder attempts to decode the message sent -by the other station using the a-priori information derived by what -has been already decoded in a previous phase of the QSO if decoding -with no a-priori information has not been successful. - -Step 1) K1JT's decoder first attempts to decode msgs of type [? ? ?] -and if this attempt fails, it attempts to decode [CQ/QRZ ? ?] or -[CQ/QRZ ?] msgs - -Step 2) if IV3NWV's decoder is unable to decode K1JT's without AP it -attempts to decode messages of the type [IV3NWV ? ?] and [IV3NWV ?]. - -Step 3) K1JT's decoder attempts to decode [? ? ?] and [K1JT IV3NWV ?] -(this last decode type has been enabled by K1JT's encoder at step 2) - -Step 4) IV3NWV's decoder attempts to decode [? ? ?] and [IV3NWV K1JT -?] (this last decode type has been enabled by IV3NWV's encoder at step -3) - -At each step the simulation reports if a decode was successful. In -this case it also reports the type of decode (see table decode_type -above) - -When mode=QRA_NOAP, only [? ? ?] decodes are attempted and no a-priori -information is used by the decoder - -The function returns 0 if all of the four messages have been decoded -by their recipients (with no retries) and -1 if any of them could not -be decoded -*/ - - int x[QRA64_K], xdec[QRA64_K]; - int y[QRA64_N]; - float *rx; - int rc; - -// Each simulated station must use its own codec since it might work with -// different a-priori information. - qra64codec *codec_iv3nwv = qra64_init(mode); // codec for IV3NWV - qra64codec *codec_k1jt = qra64_init(mode); // codec for K1JT - -// Step 1a: IV3NWV makes a CQ call (with no grid) - printf("IV3NWV tx: CQ IV3NWV\n"); - encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_BLANK); - qra64_encode(codec_iv3nwv, y, x); - rx = mfskchannel(y,channel_type,EbNodB); - -// Step 1b: K1JT attempts to decode [? ? ?], [CQ/QRZ ? ?] or [CQ/QRZ ?] - rc = qra64_decode(codec_k1jt, 0, xdec,rx); - if (rc>=0) { // decoded - printf("K1JT rx: received with apcode=%d %s\n",rc, decode_type[rc]); - -// Step 2a: K1JT replies to IV3NWV (with no grid) - printf("K1JT tx: IV3NWV K1JT\n"); - encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_BLANK); - qra64_encode(codec_k1jt, y, x); - rx = mfskchannel(y,channel_type,EbNodB); - -// Step 2b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?] or [IV3NWV ?] - rc = qra64_decode(codec_iv3nwv, 0, xdec,rx); - if (rc>=0) { // decoded - printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]); - -// Step 3a: IV3NWV replies to K1JT with a 73 - printf("IV3NWV tx: K1JT IV3NWV 73\n"); - encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV, GRID_73); - qra64_encode(codec_iv3nwv, y, x); - rx = mfskchannel(y,channel_type,EbNodB); - -// Step 3b: K1JT attempts to decode [? ? ?] or [K1JT IV3NWV ?] - rc = qra64_decode(codec_k1jt, 0, xdec,rx); - if (rc>=0) { // decoded - printf("K1JT rx: received with apcode=%d %s\n",rc, decode_type[rc]); - -// Step 4a: K1JT replies to IV3NWV with a 73 - printf("K1JT tx: IV3NWV K1JT 73\n"); - encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_73); - qra64_encode(codec_k1jt, y, x); - rx = mfskchannel(y,channel_type,EbNodB); - -// Step 4b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?], or [IV3NWV ?] - rc = qra64_decode(codec_iv3nwv, 0, xdec,rx); - if (rc>=0) { // decoded - printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]); - return 0; - } - } - } - } - printf("no decode\n"); - return -1; -} - -int test_proc_2(int channel_type, float EbNodB, int mode) -{ -/* -Here we simulate the decoder of K1JT after K1JT has sent a msg [IV3NWV K1JT] -and IV3NWV sends him the msg [K1JT IV3NWV JN66]. - -If mode=QRA_NOAP, K1JT decoder attempts to decode only msgs of type [? ? ?]. - -If mode=QRA_AUTOP, K1JT decoder will attempt to decode also the msgs -[K1JT IV3NWV] and [K1JT IV3NWV ?]. - -In the case a decode is successful the return code of the qra64_decode function -indicates the amount of a-priori information required to decode the received -message according to this table: - - rc=0 [? ? ?] AP0 - rc=1 [CQ ? ?] AP27 - rc=2 [CQ ? ] AP42 - rc=3 [CALL ? ?] AP29 - rc=4 [CALL ? ] AP44 - rc=5 [CALL CALL ?] AP57 - rc=6 [? CALL ?] AP29 - rc=7 [? CALL ] AP44 - rc=8 [CALL CALL GRID] AP72 - rc=9 [CQ CALL ?] AP55 - rc=10 [CQ CALL ] AP70 - rc=11 [CQ CALL GRID] AP70 - -The return code is <0 when decoding is unsuccessful - -This test simulates the situation ntx times and reports how many times -a particular type decode among the above 6 cases succeded. -*/ - - int x[QRA64_K], xdec[QRA64_K]; - int y[QRA64_N]; - float *rx; - float ebnodbest, ebnodbavg=0; - int rc,k; - - int ndecok[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - int nundet = 0; - int ntx = 200,ndec=0; - - qra64codec *codec_iv3nwv = qra64_init(mode); // codec for IV3NWV - qra64codec *codec_k1jt = qra64_init(mode); // codec for K1JT - - printf("\nQRA64 Test #2 - Decoding with AP knowledge (SNR-Eb/No offset = %.1f dB)\n\n", - QRA64_SNR_EBNO_OFFSET); - -// This will enable K1JT's decoder to look for calls directed to him [K1JT ? ?/b] -// printf("K1JT decoder enabled for [K1JT ? ?/blank]\n"); -// qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_MYCALL); - -// This will enable K1JT's decoder to look for IV3NWV calls directed to him [K1JT IV3NWV ?/b] -// printf("K1JT decoder enabled for [K1JT IV3NWV ?]\n"); -// qra64_apset(codec_k1jt, CALL_CQ,CALL_IV3NWV,0,APTYPE_BOTHCALLS); - -// This will enable K1JT's decoder to look for msges sent by IV3NWV [? IV3NWV ?] -// printf("K1JT decoder enabled for [? IV3NWV ?/blank]\n"); -// qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_BLANK,APTYPE_HISCALL); - -// This will enable K1JT's decoder to look for full-knowledge [K1JT IV3NWV JN66] msgs - printf("K1JT decoder enabled for [K1JT IV3NWV JN66]\n"); - qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,GRID_JN66,APTYPE_FULL); - -// This will enable K1JT's decoder to look for calls from IV3NWV [CQ IV3NWV ?/b] msgs - printf("K1JT decoder enabled for [CQ IV3NWV ?/b/JN66]\n"); - qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_JN66,APTYPE_CQHISCALL); - - - // Dx station IV3NWV calls - printf("\nIV3NWV encoder sends msg: [K1JT IV3NWV JN66]\n\n"); - encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_JN66); - -// printf("\nIV3NWV encoder sends msg: [CQ IV3NWV JN66]\n\n"); -// encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_JN66); - -// printf("\nIV3NWV encoder sends msg: [CQ IV3NWV]\n\n"); -// encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_BLANK); - qra64_encode(codec_iv3nwv, y, x); - - printf("Simulating K1JT decoder up to AP72\n"); - - for (k=0;k<ntx;k++) { - printf("."); - rx = mfskchannel(y,channel_type,EbNodB); - rc = qra64_decode(codec_k1jt, &ebnodbest, xdec,rx); - if (rc>=0) { - ebnodbavg +=ebnodbest; - if (memcmp(xdec,x,12*sizeof(int))==0) - ndecok[rc]++; - else - nundet++; - } - } - printf("\n\n"); - - - printf("Transimtted msgs:%d\nDecoded msgs:\n\n",ntx); - for (k=0;k<12;k++) { - printf("%3d with %s\n",ndecok[k],decode_type[k]); - ndec += ndecok[k]; - } - printf("\nTotal: %d/%d (%d undetected errors)\n\n",ndec,ntx,nundet); - printf(""); - - ebnodbavg/=(ndec+nundet); - printf("Estimated SNR (average in dB) = %.2f dB\n\n",ebnodbavg-QRA64_SNR_EBNO_OFFSET); - - return 0; -} - -int test_fastfading(float EbNodB, float B90, int fadingModel, int submode, int apmode, int olddec, int channel_type, int ntx) -{ - int x[QRA64_K], xdec[QRA64_K]; - int y[QRA64_N]; - float *rx; - float ebnodbest, ebnodbavg=0; - int rc,k; - float rxolddec[QRA64_N*QRA64_M]; // holds the energies at nominal tone freqs - - int ndecok[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - int nundet = 0; - int ndec=0; - - qra64codec *codec_iv3nwv; - qra64codec *codec_k1jt; - - codec_iv3nwv=qra64_init(QRA_NOAP); - codec_k1jt =qra64_init(apmode); - - if (channel_type==2) { // fast-fading case - printf("Simulating the fast-fading channel\n"); - printf("B90=%.2f Hz - Fading Model=%s - Submode=QRA64%c\n",B90,fadingModel?"Lorentz":"Gauss",submode+'A'); - printf("Decoder metric = %s\n",olddec?"AWGN":"Matched to fast-fading signal"); - } - else { - printf("Simulating the %s channel\n",channel_type?"Rayleigh block fading":"AWGN"); - printf("Decoder metric = AWGN\n"); - } - - - printf("\nEncoding msg [K1JT IV3NWV JN66]\n"); - encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV,GRID_JN66); -// printf("["); -// for (k=0;k<11;k++) printf("%02hX ",x[k]); printf("%02hX]\n",x[11]); - - qra64_encode(codec_iv3nwv, y, x); - printf("%d transmissions will be simulated\n\n",ntx); - - if (apmode==QRA_USERAP) { - // This will enable K1JT's decoder to look for cq/qrz calls [CQ/QRZ ? ?/b] - printf("K1JT decoder enabled for [CQ ? ?/blank]\n"); - qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_CQQRZ); - - // This will enable K1JT's decoder to look for calls directed to him [K1JT ? ?/b] - printf("K1JT decoder enabled for [K1JT ? ?/blank]\n"); - qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_MYCALL); - - // This will enable K1JT's decoder to look for msges sent by IV3NWV [? IV3NWV ?] - printf("K1JT decoder enabled for [? IV3NWV ?/blank]\n"); - qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_BLANK,APTYPE_HISCALL); - - // This will enable K1JT's decoder to look for IV3NWV calls directed to him [K1JT IV3NWV ?/b] - printf("K1JT decoder enabled for [K1JT IV3NWV ?]\n"); - qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,0,APTYPE_BOTHCALLS); - - // This will enable K1JT's decoder to look for full-knowledge [K1JT IV3NWV JN66] msgs - printf("K1JT decoder enabled for [K1JT IV3NWV JN66]\n"); - qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,GRID_JN66,APTYPE_FULL); - - // This will enable K1JT's decoder to look for calls from IV3NWV [CQ IV3NWV ?/b] msgs - printf("K1JT decoder enabled for [CQ IV3NWV ?/b/JN66]\n"); - qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_JN66,APTYPE_CQHISCALL); - - } - - printf("\nNow decoding with K1JT's decoder...\n"); -/* - if (channel_type==2) // simulate a fast-faded signal - printf("Simulating a fast-fading channel with given B90 and spread type\n"); - else - printf("Simulating a %s channel\n",channel_type?"Rayleigh block fading":"AWGN"); -*/ - for (k=0;k<ntx;k++) { - - if ((k%10)==0) - printf(" %5.1f %%\r",100.0*k/ntx); -// printf("."); // work in progress - - if (channel_type==2) { - // generate a fast-faded signal - rc = qra64_fastfading_channel(&rx,y,submode,EbNodB,B90,fadingModel); - if (rc<0) { - printf("\nqra64_fastfading_channel error. rc=%d\n",rc); - return -1; - } - } - else // generate a awgn or Rayleigh block fading signal - rx = mfskchannel(y, channel_type, EbNodB); - - - if (channel_type==2) // fast-fading case - if (olddec==1) { - int k, j; - int jj = 1<<submode; - int bps = QRA64_M*(2+jj); - float *rxbase; - float *out = rxolddec; - // calc energies at nominal freqs - for (k=0;k<QRA64_N;k++) { - rxbase = rx + QRA64_M + k*bps; - for (j=0;j<QRA64_M;j++) { - *out++=*rxbase; - rxbase+=jj; - } - } - // decode with awgn decoder - rc = qra64_decode(codec_k1jt,&ebnodbest,xdec,rxolddec); - } - else // use fast-fading decoder - rc = qra64_decode_fastfading(codec_k1jt,&ebnodbest,xdec,rx,submode,B90,fadingModel); - else // awgn or rayleigh channel. use the old decoder whatever the olddec option is - rc = qra64_decode(codec_k1jt,&ebnodbest,xdec,rx); - - - - if (rc>=0) { - ebnodbavg +=ebnodbest; - if (memcmp(xdec,x,12*sizeof(int))==0) - ndecok[rc]++; - else { - fprintf(stderr,"\nUndetected error with rc=%d\n",rc); - nundet++; - } - } - - } - printf(" %5.1f %%\r",100.0*k/ntx); - - printf("\n\n"); - - printf("Msgs transmitted:%d\nMsg decoded:\n\n",ntx); - for (k=0;k<12;k++) { - printf("rc=%2d %3d with %s\n",k,ndecok[k],decode_type[k]); - ndec += ndecok[k]; - } - printf("\nTotal: %d/%d (%d undetected errors)\n\n",ndec,ntx,nundet); - printf(""); - - if (ndec>0) { - ebnodbavg/=(ndec+nundet); - printf("Estimated SNR (average in dB) = %.2f dB\n\n",ebnodbavg-QRA64_SNR_EBNO_OFFSET); - } - - return 0; -} - - - -void syntax(void) -{ - - printf("\nQRA64 Mode Tests\n"); - printf("2016, Nico Palermo - IV3NWV\n\n"); - printf("---------------------------\n\n"); - printf("Syntax: qra64 [-s<snrdb>] [-c<channel>] [-a<ap-type>] [-t<testtype>] [-h]\n"); - printf("Options: \n"); - printf(" -s<snrdb> : set simulation SNR in 2500 Hz BW (default:-27.5 dB)\n"); - printf(" -c<channel> : set channel type 0=AWGN (default) 1=Rayleigh 2=Fast-fading\n"); - printf(" -a<ap-type> : set decode type 0=NOAP 1=AUTOAP (default) 2=USERAP\n"); - printf(" -t<testtype>: 0=simulate seq of msgs between IV3NWV and K1JT (default)\n"); - printf(" 1=simulate K1JT receiving K1JT IV3NWV JN66\n"); - printf(" 2=simulate fast-fading/awgn/rayliegh decoders performance\n"); - printf(" -n<ntx> : simulate the transmission of ntx codewords (default=100)\n"); - - printf("Options used only for fast-fading simulations (-c2):\n"); - printf(" -b : 90%% fading bandwidth in Hz [1..230 Hz] (default = 2.5 Hz)\n"); - printf(" -m : fading model. 0=Gauss, 1=Lorentz (default = Lorentz)\n"); - printf(" -q : qra64 submode. 0=QRA64A,... 4=QRA64E (default = QRA64A)\n"); - printf(" -d : use the old awgn decoder\n"); - printf(" -h: this help\n"); - printf("Example:\n"); - printf(" qra64 -t2 -c2 -a2 -b50 -m1 -q2 -n10000 -s-26\n"); - printf(" runs the error performance test (-t2)\n"); - printf(" with USER_AP (-a2)\n"); - printf(" simulating a fast fading channel (-c2)\n"); - printf(" with B90 = 50 Hz (-b50), Lorentz Doppler (-m1), mode QRA64C (-q2)\n"); - printf(" ntx = 10000 codewords (-n10000) and SNR = -26 dB (-s-26)\n"); - -} - -int main(int argc, char* argv[]) -{ - int k, rc, nok=0; - float SNRdB = -27.5f; - unsigned int channel = CHANNEL_AWGN; - unsigned int mode = QRA_AUTOAP; - unsigned int testtype=0; - int nqso = 100; - float EbNodB; - float B90 = 2.5; - int fadingModel = 1; - int submode = 0; - int olddec = 0; - int ntx = 100; - -// Parse the command line - while(--argc) { - argv++; - - if (strncmp(*argv,"-h",2)==0) { - syntax(); - return 0; - } - else - if (strncmp(*argv,"-n",2)==0) { - ntx = ( int)atoi((*argv)+2); - if (ntx<100 || ntx>1000000) { - printf("Invalid -n option. ntx must be in the range [100..1000000]\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-a",2)==0) { - mode = ( int)atoi((*argv)+2); - if (mode>2) { - printf("Invalid decoding mode\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-s",2)==0) { - SNRdB = (float)atof((*argv)+2); - if (SNRdB>20 || SNRdB<-50) { - printf("SNR should be in the range [-50..20]\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-t",2)==0) { - testtype = ( int)atoi((*argv)+2); - if (testtype>2) { - printf("Invalid test type\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-c",2)==0) { - channel = ( int)atoi((*argv)+2); - if (channel>CHANNEL_FASTFADE) { - printf("Invalid channel type\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-b",2)==0) { - B90 = (float)atof((*argv)+2); - if (B90<1 || B90>230) { - printf("Invalid B90\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-m",2)==0) { - fadingModel = (int)atoi((*argv)+2); - if (fadingModel<0 || fadingModel>1) { - printf("Invalid fading model\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-q",2)==0) { - submode = (int)atoi((*argv)+2); - if (submode<0 || submode>4) { - printf("Invalid submode\n"); - syntax(); - return -1; - } - } - else - if (strncmp(*argv,"-d",2)==0) { - olddec = 1; - } - else { - printf("Invalid option\n"); - syntax(); - return -1; - } - } - - if (testtype<2) // old tests - if (channel==CHANNEL_FASTFADE) { - printf("Invalid Option. Test type 0 and 1 supports only AWGN or Rayleigh Channel model\n"); - return -1; - } - - EbNodB = SNRdB+QRA64_SNR_EBNO_OFFSET; - -#if defined(__linux__) || defined(__unix__) - srand48(GetTickCount()); -#endif - - if (testtype==0) { - for (k=0;k<nqso;k++) { - printf("\n\n------------------------\n"); - rc = test_proc_1(channel, EbNodB, mode); - if (rc==0) - nok++; - } - printf("\n\n%d/%d QSOs to end without repetitions\n",nok,nqso); - printf("Input SNR = %.1fdB channel=%s ap-mode=%s\n\n", - SNRdB, - channel==CHANNEL_AWGN?"AWGN":"RAYLEIGH", - apmode_type[mode] - ); - } - else if (testtype==1) { - test_proc_2(channel, EbNodB, mode); - printf("Input SNR = %.1fdB channel=%s ap-mode=%s\n\n", - SNRdB, - channel==CHANNEL_AWGN?"AWGN":"RAYLEIGH", - apmode_type[mode] - ); - } - else { - printf("Input SNR = %.1fdB ap-mode=%s\n\n", - SNRdB, - apmode_type[mode] - ); - test_fastfading(EbNodB,B90,fadingModel,submode,mode,olddec, channel, ntx); - } - return 0; -} diff --git a/lib/qra/qra64/qra64.c b/lib/qra/qra64/qra64.c deleted file mode 100644 index 6470cedbd..000000000 --- a/lib/qra/qra64/qra64.c +++ /dev/null @@ -1,1121 +0,0 @@ -/* -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)) - 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; -} diff --git a/lib/qra/qra64/qra64.h b/lib/qra/qra64/qra64.h deleted file mode 100644 index 3b2b5bd3d..000000000 --- a/lib/qra/qra64/qra64.h +++ /dev/null @@ -1,269 +0,0 @@ -// 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_ diff --git a/lib/qra/qra64/qra64_all.c b/lib/qra/qra64/qra64_all.c deleted file mode 100644 index 28f8ab928..000000000 --- a/lib/qra/qra64/qra64_all.c +++ /dev/null @@ -1,1050 +0,0 @@ -/* -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; -} diff --git a/lib/qra/qra64/qra64_subs.c b/lib/qra/qra64/qra64_subs.c deleted file mode 100644 index b60f9fafa..000000000 --- a/lib/qra/qra64/qra64_subs.c +++ /dev/null @@ -1,65 +0,0 @@ -// 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; - } -} - diff --git a/lib/qra/qra64/qra64example.txt b/lib/qra/qra64/qra64example.txt deleted file mode 100644 index 3add33d36..000000000 --- a/lib/qra/qra64/qra64example.txt +++ /dev/null @@ -1,88 +0,0 @@ -$ 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 diff --git a/lib/qra/qra64/qra64sim.f90 b/lib/qra/qra64/qra64sim.f90 deleted file mode 100644 index 2360ef7da..000000000 --- a/lib/qra/qra64/qra64sim.f90 +++ /dev/null @@ -1,201 +0,0 @@ -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 - complex c00(0:720000) !Analytic signal for dat() - real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq - character msg*22,fname*11,csubmode*1,arg*12,cd*1 - character msgsent*22 - logical lsync - data lsync/.false./ - - nargs=iargc() - if(nargs.ne.8) then - print*,'Usage: qra64sim "msg" A-E Nsigs fDop DT Nfiles Sync SNR' - print*,'Example qra64sim "K1ABC W9XYZ EN37" A 10 0.2 0.0 1 T -26' - print*,'Sync = T to include sync test.' - 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) - if(arg(1:1).eq.'T' .or. arg(1:1).eq.'1') lsync=.true. - call getarg(8,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('-')) - - nsync=0 - do ifile=1,nfiles !Loop over requested number of files - write(fname,1002) ifile !Output filename -1002 format('000000_',i4.4) - open(10,file=fname//'.wav',access='stream',status='unknown') - xnoise=0. - cdat=0. - if(snrdb.lt.90) then - do i=1,npts - xnoise(i)=gran() !Generate gaussian noise - enddo - endif - - 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 !Use truncated Lorenzian shape for fspread - 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) - - if(lsync) then - cd=' ' - if(ifile.eq.nfiles) cd='d' - nf1=200 - nf2=3000 - nfqso=nint(f0) - ntol=100 - minsync=0 - emedelay=0.0 - call ana64(dat,npts,c00) - call sync64(c00,nf1,nf2,nfqso,ntol,minsync,mode64,emedelay,xdt2,f02, & - jpk0,sync,sync2,width) - terr=1.01/(8.0*baud) - ferr=1.01*mode64*baud - if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 - open(40,file='sync64.out',status='unknown',position='append') - write(40,1030) ifile,64,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & - width,sync,sync2,nsync,cd -1030 format(i4,i3,1x,a1,2f7.1,f7.2,4f8.1,i5,1x,a1) - close(40) - endif - enddo - if(lsync) write(*,1040) snrdb,nfiles,nsync -1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5) - -999 end program qra64sim diff --git a/lib/qra_loops.f90 b/lib/qra_loops.f90 deleted file mode 100644 index cebd7036a..000000000 --- a/lib/qra_loops.f90 +++ /dev/null @@ -1,137 +0,0 @@ -subroutine qra_loops(c00,npts2,nsps,mode,mode64,nsubmode,nFadingModel, & - ndepth,nc1,nc2,ng2,naptype,jpk0,xdt,f0,width,snr2,irc,dat4) - - use packjt - use timer_module, only: timer - parameter (LN=2176*63) !LN=LL*NN; LL = 64*(mode64+2) - character*37 decoded - complex c00(0:npts2-1) !Analytic representation of dd(), 6000 Hz - complex ,allocatable :: c0(:) !Ditto, with freq shift - real a(3) !twkfreq params f,f1,f2 - real s3(LN) !Symbol spectra - real s3avg(LN) !Averaged symbol spectra - integer dat4(12),dat4x(12) !Decoded message (as 12 integers) - integer nap(0:11) !AP return codes - data nap/0,2,3,2,3,4,2,3,6,4,6,6/,nsave/0/ - save nsave,s3avg - - allocate(c0(0:npts2-1)) - irc=-99 - s3lim=20. - ibwmax=11 - if(mode64.le.4) ibwmax=9 - ibwmin=ibwmax - idtmax=3 - call qra_params(ndepth,maxaptype,idfmax,idtmax,ibwmin,ibwmax,maxdist) - LL=64*(mode64+2) - NN=63 - napmin=99 - ncall=0 - - do iavg=0,1 - if(iavg.eq.1) then - idfmax=1 - idtmax=1 - endif - do idf=1,idfmax - ndf=idf/2 - if(mod(idf,2).eq.0) ndf=-ndf - a=0. - a(1)=-(f0+0.4*ndf) - call twkfreq(c00,c0,npts2,6000.0,a) - do idt=1,idtmax - ndt=idt/2 - if(iavg.eq.0) then - if(mod(idt,2).eq.0) ndt=-ndt - jpk=jpk0 + 240*ndt !240/6000 = 0.04 s = tsym/32 - if(jpk.lt.0) jpk=0 - call timer('spec64 ',0) - call spec64(c0,nsps,mode,mode64,jpk,s3,LL,NN) - call timer('spec64 ',1) - call pctile(s3,LL*NN,40,base) - s3=s3/base - where(s3(1:LL*NN)>s3lim) s3(1:LL*NN)=s3lim - else - s3(1:LL*NN)=s3avg(1:LL*NN) - endif - do ibw=ibwmax,ibwmin,-2 - ndist=ndf**2 + ndt**2 + ((ibwmax-ibw)/2)**2 - if(ndist.gt.maxdist) cycle - b90=1.728**ibw - if(b90.gt.230.0) cycle - if(b90.lt.0.15*width) exit - ncall=ncall+1 - call timer('qra64_de',0) - call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90, & - nFadingModel,dat4,snr2,irc) - call timer('qra64_de',1) - if(irc.eq.0) go to 200 - if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2) - iirc=max(0,min(irc,11)) - if(irc.gt.0 .and. nap(iirc).lt.napmin) then - dat4x=dat4 - b90x=b90 - snr2x=snr2 - napmin=nap(iirc) - irckeep=irc - xdtkeep=jpk/6000.0 - 1.0 - f0keep=-a(1) - idfkeep=idf - idtkeep=idt - ibwkeep=ibw - ndistx=ndist - go to 100 !### - endif - enddo ! ibw (b90 loop) - !### if(iand(ndepth,3).lt.3 .and. irc.ge.0) go to 100 - enddo ! idt (DT loop) - enddo ! idf (f0 loop) -! if(iavg.eq.0 .and. abs(jpk0-4320).le.1300) then - if(iavg.eq.0) then - a=0. - a(1)=-f0 - call twkfreq(c00,c0,npts2,6000.0,a) - jpk=3000 !### These definitions need work ### -! if(nsps.ge.3600) jpk=4080 !### - if(nsps.ge.3600) jpk=6000 !### - call spec64(c0,nsps,mode,mode64,jpk,s3,LL,NN) - call pctile(s3,LL*NN,40,base) - s3=s3/base - where(s3(1:LL*NN)>s3lim) s3(1:LL*NN)=s3lim - s3avg(1:LL*NN)=s3avg(1:LL*NN)+s3(1:LL*NN) - nsave=nsave+1 - endif - if(iavg.eq.0 .and. nsave.lt.2) exit - enddo ! iavg - -100 if(napmin.ne.99) then - dat4=dat4x - b90=b90x - snr2=snr2x - irc=irckeep - xdt=xdtkeep - f0=f0keep - idt=idtkeep - idf=idfkeep - ibw=ibwkeep - ndist=ndistx - endif - -200 if(mode.eq.65 .and. nsps.eq.7200/2) xdt=xdt+0.4 !### Empirical -- WHY ??? ### - - if(irc.ge.0) then - navg=nsave - if(iavg.eq.0) navg=0 - !### For tests only: - open(53,file='fort.53',status='unknown',position='append') - call unpackmsg(dat4,decoded) !Unpack the user message - write(53,3053) idf,idt,ibw,b90,xdt,f0,snr2,ndist,irc,navg,decoded(1:22) -3053 format(3i5,f7.1,f7.2,2f7.1,3i4,2x,a22) - close(53) - !### - nsave=0 - s3avg=0. - irc=irc + 100*navg - endif - return -end subroutine qra_loops