From be640ffe80feac17a7e3fc12f8989b1bd4753d5b Mon Sep 17 00:00:00 2001 From: Xerbo Date: Mon, 18 Jul 2022 23:24:14 +0100 Subject: [PATCH] Convert entire code base to single precision --- .github/workflows/build.yml | 4 +- CMakeLists.txt | 4 +- src/dsp.c | 148 ++++++++++++++---------------------- src/filter.c | 83 +++++++++----------- src/filter.h | 112 +++++---------------------- src/image.c | 34 ++++----- src/taps.h | 92 ++++++++++++++++++++++ src/util.h | 20 +++++ 8 files changed, 248 insertions(+), 249 deletions(-) create mode 100644 src/taps.h create mode 100644 src/util.h diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index a468906..f6b4659 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -18,7 +18,9 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 + with: + submodules: 'recursive' - name: Install dependencies # The user does not run as root diff --git a/CMakeLists.txt b/CMakeLists.txt index 36b1f66..f783b0e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,7 +41,7 @@ if (PNG_FOUND AND LIBSNDFILE_FOUND) if(CMAKE_BUILD_TYPE MATCHES "Release") target_compile_options(aptdec PRIVATE -Wall -Wextra -pedantic -Wno-missing-field-initializers) else() - target_compile_options(aptdec PRIVATE -Wall -Wextra -pedantic) + target_compile_options(aptdec PRIVATE -Wall -Wextra -pedantic -Wno-missing-field-initializers) endif() endif() else() @@ -59,7 +59,7 @@ else() if(CMAKE_BUILD_TYPE MATCHES "Release") target_compile_options(apt PRIVATE -Wall -Wextra -pedantic -Wno-missing-field-initializers) else() - target_compile_options(apt PRIVATE -Wall -Wextra -pedantic) + target_compile_options(apt PRIVATE -Wall -Wextra -pedantic -Wno-missing-field-initializers) endif() endif() diff --git a/src/dsp.c b/src/dsp.c index d62ec3b..bddcf3a 100755 --- a/src/dsp.c +++ b/src/dsp.c @@ -1,20 +1,19 @@ -/* - * This file is part of Aptdec. - * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2004-2009 Thierry Leconte (F4DWV) 2019-2022 Xerbo (xerbo@protonmail.com) * - * Aptdec 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 2 of the License, or - * (at your option) any later version. + * This program 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 2 of the License, or + * (at your option) any later version. * - * This program 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 this program. If not, see . + * This program 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 this program. If not, see . */ #include @@ -24,6 +23,7 @@ #include "apt.h" #include "filter.h" +#include "taps.h" // In case your C compiler is so old that Pi hadn't been invented yet #ifndef M_PI @@ -41,13 +41,13 @@ #define RSMULT 15 #define Fi (Fp * RSMULT) -static double Fe; +static float Fe; -static double offset = 0.0; -static double FreqLine = 1.0; +static float offset = 0.0; +static float FreqLine = 1.0; -static double FreqOsc; -static double K1, K2; +static float FreqOsc; +static float K1, K2; // Check the sample rate and calculate some constants int apt_init(double F) { @@ -63,56 +63,25 @@ int apt_init(double F) { return 0; } -/* Fast phase estimator - * Calculates the phase angle of a signal from a IQ sample - */ -static inline double Phase(double I, double Q) { - double angle, r; - int s; - - if(I == 0.0 && Q == 0.0) return 0.0; - - if (Q < 0) { - s = -1; - Q = -Q; - } else { - s = 1; - } - - if (I >= 0) { - r = (I - Q) / (I + Q); - angle = 0.25 - 0.25 * r; - } else { - r = (I + Q) / (Q - I); - angle = 0.75 - 0.25 * r; - } - - if(s > 0){ - return angle; - }else{ - return -angle; - } -} - /* Phase locked loop * https://arachnoid.com/phase_locked_loop/ * Model of this filter here https://www.desmos.com/calculator/m0uadgkoee */ -static double pll(double I, double Q) { +static float pll(float I2, float Q) { // PLL coefficient - static double PhaseOsc = 0.0; - double Io, Qo; - double Ip, Qp; - double DPhi; + static float PhaseOsc = 0.0; + float Io, Qo; + float Ip, Qp; + float DPhi; // Quadrature oscillator / reference Io = cos(PhaseOsc); Qo = sin(PhaseOsc); // Phase detector - Ip = I * Io + Q * Qo; - Qp = Q * Io - I * Qo; - DPhi = Phase(Ip, Qp); + Ip = I2 * Io + Q * Qo; + Qp = Q * Io - I2 * Qo; + DPhi = atan2f(Qp, Ip); // Loop filter PhaseOsc += 2.0 * M_PI * (K1 * DPhi + FreqOsc); @@ -131,16 +100,16 @@ static double pll(double I, double Q) { } // Convert samples into pixels -static int getamp(double *ampbuff, int count, apt_getsamples_t getsamples, void *context) { +static int getamp(float *ampbuff, int count, apt_getsamples_t getsamples, void *context) { static float inbuff[BLKIN]; static int idxin = 0; static int nin = 0; for (int n = 0; n < count; n++) { - double I, Q; + float I2, Q; // Get some more samples when needed - if (nin < IQFilterLen * 2 + 2) { + if (nin < HILBERT_FILTER_SIZE * 2 + 2) { // Number of samples read int res; memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); @@ -151,13 +120,15 @@ static int getamp(double *ampbuff, int count, apt_getsamples_t getsamples, void nin += res; // Make sure there is enough samples to continue - if (nin < IQFilterLen * 2 + 2) + if (nin < HILBERT_FILTER_SIZE * 2 + 2) return n; } // Process read samples into a brightness value - iqfir(&inbuff[idxin], iqfilter, IQFilterLen, &I, &Q); - ampbuff[n] = pll(I, Q); + float complex tmp = hilbert_transform(&inbuff[idxin], hilbert_filter, HILBERT_FILTER_SIZE); + I2 = crealf(tmp); + Q = cimagf(tmp); + ampbuff[n] = pll(I2, Q); // Increment current sample idxin++; @@ -167,25 +138,25 @@ static int getamp(double *ampbuff, int count, apt_getsamples_t getsamples, void return count; } -// Sub-pixel offsetting + FIR compensation +// Sub-pixel offsetting int getpixelv(float *pvbuff, int count, apt_getsamples_t getsamples, void *context) { // Amplitude buffer - static double ampbuff[BLKAMP]; + static float ampbuff[BLKAMP]; static int nam = 0; static int idxam = 0; - double mult; + float mult; // Gaussian resampling factor - mult = (double) Fi / Fe * FreqLine; - int m = (int)(RSFilterLen / mult + 1); + mult = (float) Fi / Fe * FreqLine; + int m = (int)(LOW_PASS_SIZE / mult + 1); for (int n = 0; n < count; n++) { int shift; if (nam < m) { int res; - memmove(ampbuff, &(ampbuff[idxam]), nam * sizeof(double)); + memmove(ampbuff, &(ampbuff[idxam]), nam * sizeof(float)); idxam = 0; res = getamp(&(ampbuff[nam]), BLKAMP - nam, getsamples, context); nam += res; @@ -193,8 +164,7 @@ int getpixelv(float *pvbuff, int count, apt_getsamples_t getsamples, void *conte return n; } - // Gaussian FIR compensation filter - pvbuff[n] = (float)(rsfir(&(ampbuff[idxam]), rsfilter, RSFilterLen, offset, mult) * mult * 256.0); + pvbuff[n] = interpolating_convolve(&(ampbuff[idxam]), low_pass, LOW_PASS_SIZE, offset, mult) * mult * 256.0; shift = ((int) floor((RSMULT - offset) / mult)) + 1; offset = shift * mult + offset - RSMULT; @@ -208,15 +178,15 @@ int getpixelv(float *pvbuff, int count, apt_getsamples_t getsamples, void *conte // Get an entire row of pixels, aligned with sync markers int apt_getpixelrow(float *pixelv, int nrow, int *zenith, int reset, apt_getsamples_t getsamples, void *context) { - static float pixels[PixelLine + SyncFilterLen]; - static int npv; + static float pixels[PixelLine + SYNC_PATTERN_SIZE]; + static size_t npv; static int synced = 0; - static double max = 0.0; - static double minDoppler = 1000000000, previous = 0; + static float max = 0.0; + static float minDoppler = 1000000000, previous = 0; if(reset) synced = 0; - double corr, ecorr, lcorr; + float corr, ecorr, lcorr; int res; // Move the row buffer into the the image buffer @@ -224,20 +194,20 @@ int apt_getpixelrow(float *pixelv, int nrow, int *zenith, int reset, apt_getsamp memmove(pixelv, pixels, npv * sizeof(float)); // Get the sync line - if (npv < SyncFilterLen + 2) { - res = getpixelv(&(pixelv[npv]), SyncFilterLen + 2 - npv, getsamples, context); + if (npv < SYNC_PATTERN_SIZE + 2) { + res = getpixelv(&(pixelv[npv]), SYNC_PATTERN_SIZE + 2 - npv, getsamples, context); npv += res; - if (npv < SyncFilterLen + 2) + if (npv < SYNC_PATTERN_SIZE + 2) return 0; } // Calculate the frequency offset - ecorr = fir(pixelv, Sync, SyncFilterLen); - corr = fir(&pixelv[1], Sync, SyncFilterLen - 1); - lcorr = fir(&pixelv[2], Sync, SyncFilterLen - 2); + ecorr = convolve(pixelv, sync_pattern, SYNC_PATTERN_SIZE); + corr = convolve(&pixelv[1], sync_pattern, SYNC_PATTERN_SIZE - 1); + lcorr = convolve(&pixelv[2], sync_pattern, SYNC_PATTERN_SIZE - 2); FreqLine = 1.0+((ecorr-lcorr) / corr / PixelLine / 4.0); - double val = fabs(lcorr - ecorr)*0.25 + previous*0.75; + float val = fabs(lcorr - ecorr)*0.25 + previous*0.75; if(val < minDoppler && nrow > 10){ minDoppler = val; *zenith = nrow; @@ -255,19 +225,19 @@ int apt_getpixelrow(float *pixelv, int nrow, int *zenith, int reset, apt_getsamp int mshift; static int lastmshift; - if (npv < PixelLine + SyncFilterLen) { - res = getpixelv(&(pixelv[npv]), PixelLine + SyncFilterLen - npv, getsamples, context); + if (npv < PixelLine + SYNC_PATTERN_SIZE) { + res = getpixelv(&(pixelv[npv]), PixelLine + SYNC_PATTERN_SIZE - npv, getsamples, context); npv += res; - if (npv < PixelLine + SyncFilterLen) + if (npv < PixelLine + SYNC_PATTERN_SIZE) return 0; } // Test every possible position until we get the best result mshift = 0; for (int shift = 0; shift < PixelLine; shift++) { - double corr; + float corr; - corr = fir(&(pixelv[shift + 1]), Sync, SyncFilterLen); + corr = convolve(&(pixelv[shift + 1]), sync_pattern, SYNC_PATTERN_SIZE); if (corr > max) { mshift = shift; max = corr; diff --git a/src/filter.c b/src/filter.c index 300fad4..fac6e40 100755 --- a/src/filter.c +++ b/src/filter.c @@ -1,65 +1,56 @@ -/* - * This file is part of Aptdec. - * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2004-2009 Thierry Leconte (F4DWV) 2019-2022 Xerbo (xerbo@protonmail.com) * - * Aptdec 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 2 of the License, or - * (at your option) any later version. + * This program 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 2 of the License, or + * (at your option) any later version. * - * This program 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 this program. If not, see . + * This program 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 this program. If not, see . */ #include +#include "filter.h" +#include "util.h" -// Finite impulse response -float fir(float *buff, const float *coeff, const int len) { - double r; - - r = 0.0; - for (int i = 0; i < len; i++) { - r += buff[i] * coeff[i]; +float convolve(const float *in, const float *taps, size_t len) { + float sum = 0.0; + for (size_t i = 0; i < len; i++) { + sum += in[i] * taps[i]; } - return (float)r; + + return sum; } -/* IQ finite impulse response - * Turn samples into a single IQ sample - */ -void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q) { - double i = 0.0, q = 0.0; +float complex hilbert_transform(const float *in, const float *taps, size_t len) { + float i = 0.0; + float q = 0.0; - for (int k = 0; k < len; k++) { - q += buff[2*k] * coeff[k]; - i += buff[2*k]; + for (size_t k = 0; k < len; k++) { + q += in[2*k] * taps[k]; + i += in[2*k]; } - i = buff[len-1] - (i / len); - *I = i, *Q = q; + i = in[len-1] - (i / len); + return i + q*I; } -/* Gaussian finite impulse responce compensation - * https://www.recordingblogs.com/wiki/gaussian-window - */ -float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta) { - double out; +float interpolating_convolve(const float *in, const float *taps, size_t len, float offset, float delta) { + float out = 0.0; + float n = offset; - out = 0.0; - double n = offset; - for (int i = 0; i < (len-1)/delta-1; n += delta, i++) { - int k; - double alpha; + for (size_t i = 0; i < (len-1)/delta-1; n += delta, i++) { + int k = (int)floor(n); + float alpha = n - k; - k = (int)floor(n); - alpha = n - k; - out += buff[i] * (coeff[k] * (1.0 - alpha) + coeff[k + 1] * alpha); + out += in[i] * (taps[k] * (1.0f-alpha) + taps[k + 1] * alpha); } - return (float)out; + return out; } diff --git a/src/filter.h b/src/filter.h index 62c0aa1..7ff0b61 100755 --- a/src/filter.h +++ b/src/filter.h @@ -1,100 +1,24 @@ -/* - * This file is part of Aptdec. - * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2019-2022 Xerbo (xerbo@protonmail.com) * - * Aptdec 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 2 of the License, or - * (at your option) any later version. + * This program 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 2 of the License, or + * (at your option) any later version. * - * This program 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 this program. If not, see . + * This program 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 this program. If not, see . */ -float fir(float *buff, const float *coeff, const int len); -void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q); -float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta); - -// IQ finite impulse response filter -#define IQFilterLen 32 -const float iqfilter[IQFilterLen] = { 0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, -0.0335063, 0.0374482, 0.0424413, 0.0489708, 0.0578745, 0.0707355, 0.0909457, 0.127324, 0.212207, 0.63662, --0.63662, -0.212207, -0.127324, -0.0909457, -0.0707355, -0.0578745, -0.0489708, -0.0424413, -0.0374482, --0.0335063, -0.0303152, -0.0276791, -0.0254648, -0.0235785, -0.0219524, -0.0205361 }; - -// Pattern of a NOAA sync line -#define SyncFilterLen 32 -const float Sync[SyncFilterLen] = { -14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14 }; +#include +#include -// Gaussian finite impulse response compensation filter -#define RSFilterLen 437 -const float rsfilter[RSFilterLen] = { -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, --3.75376e-04, -6.95337e-04, -5.93148e-04, -8.79730e-04, -8.15327e-04, -1.05669e-03, -1.01377e-03, --1.19836e-03, -1.15443e-03, -1.26937e-03, -1.20955e-03, -1.23904e-03, -1.15302e-03, -1.08660e-03, --9.64235e-04, -8.02450e-04, -6.46202e-04, -3.95376e-04, -2.18096e-04, 1.11906e-04, 2.89567e-04, -6.67167e-04, 8.19039e-04, 1.21725e-03, 1.30556e-03, 1.69365e-03, 1.68588e-03, 2.03277e-03, 1.90159e-03, -2.18455e-03, 1.90833e-03, 2.12100e-03, 1.69052e-03, 1.77484e-03, 1.42542e-03, 1.18292e-03, 8.66979e-04, -5.54161e-04, 2.15793e-04, -1.11623e-04, -4.35173e-04, -7.27194e-04, -9.91551e-04, -1.20407e-03, --1.37032e-03, -1.46991e-03, -1.51120e-03, -1.48008e-03, -1.39047e-03, -1.23115e-03, -1.02128e-03, --7.60099e-04, -4.68008e-04, -1.46339e-04, 1.80867e-04, 5.11244e-04, 8.19243e-04, 1.09739e-03, 1.32668e-03, -1.50632e-03, 1.61522e-03, 1.66246e-03, 1.62390e-03, 1.52430e-03, 1.34273e-03, 1.10736e-03, 8.10335e-04, -4.76814e-04, 1.13622e-04, -2.64150e-04, -6.26595e-04, -9.95436e-04, -1.27846e-03, -1.54080e-03, --1.74292e-03, -1.86141e-03, -1.89318e-03, -1.83969e-03, -1.69770e-03, -1.47938e-03, -1.18696e-03, --8.37003e-04, -4.39507e-04, -1.56907e-05, 4.19904e-04, 8.43172e-04, 1.23827e-03, 1.58411e-03, -1.86382e-03, 2.06312e-03, 2.17177e-03, 2.18121e-03, 2.08906e-03, 1.89772e-03, 1.61153e-03, -1.24507e-03, 8.13976e-04, 3.29944e-04, -1.74591e-04, -6.83619e-04, -1.17826e-03, -1.61659e-03, --2.00403e-03, -2.29070e-03, -2.49179e-03, -2.56546e-03, -2.53448e-03, -2.37032e-03, -2.10060e-03, --1.72140e-03, -1.24542e-03, -7.15425e-04, -1.24964e-04, 4.83736e-04, 1.08328e-03, 1.64530e-03, -2.14503e-03, 2.55400e-03, 2.85589e-03, 3.02785e-03, 3.06271e-03, 2.95067e-03, 2.69770e-03, -2.30599e-03, 1.79763e-03, 1.18587e-03, 5.04003e-04, -2.23591e-04, -9.57591e-04, -1.66939e-03, --2.31717e-03, -2.87636e-03, -3.31209e-03, -3.60506e-03, -3.73609e-03, -3.69208e-03, -3.44913e-03, --3.06572e-03, -2.50229e-03, -1.80630e-03, -1.00532e-03, -1.22305e-04, 7.83910e-04, 1.69402e-03, -2.53826e-03, 3.30312e-03, 3.91841e-03, 4.38017e-03, 4.63546e-03, 4.68091e-03, 4.50037e-03, -4.09614e-03, 3.47811e-03, 2.67306e-03, 1.70418e-03, 6.20542e-04, -5.36994e-04, -1.70981e-03, --2.84712e-03, -3.88827e-03, -4.78659e-03, -5.48593e-03, -5.95049e-03, -6.14483e-03, -6.05118e-03, --5.65829e-03, -4.97525e-03, -4.01796e-03, -2.82224e-03, -1.43003e-03, 1.00410e-04, 1.71169e-03, -3.31983e-03, 4.87796e-03, 6.23237e-03, 7.31013e-03, 8.20642e-03, 8.67374e-03, 8.77681e-03, -8.43444e-03, 7.66794e-03, 6.46827e-03, 4.87294e-03, 2.92923e-03, 6.98913e-04, -1.72126e-03, --4.24785e-03, -6.75380e-03, -9.13309e-03, -1.12532e-02, -1.30038e-02, -1.42633e-02, -1.49338e-02, --1.49145e-02, -1.41484e-02, -1.25761e-02, -1.01870e-02, -6.97432e-03, -2.97910e-03, 1.75386e-03, -7.11899e-03, 1.30225e-02, 1.93173e-02, 2.58685e-02, 3.24965e-02, 3.90469e-02, 4.53316e-02, -5.11931e-02, 5.64604e-02, 6.09924e-02, 6.46584e-02, 6.73547e-02, 6.90049e-02, 6.97096e-02, -6.90049e-02, 6.73547e-02, 6.46584e-02, 6.09924e-02, 5.64604e-02, 5.11931e-02, 4.53316e-02, -3.90469e-02, 3.24965e-02, 2.58685e-02, 1.93173e-02, 1.30225e-02, 7.11899e-03, 1.75386e-03, --2.97910e-03, -6.97432e-03, -1.01870e-02, -1.25761e-02, -1.41484e-02, -1.49145e-02, -1.49338e-02, --1.42633e-02, -1.30038e-02, -1.12532e-02, -9.13309e-03, -6.75380e-03, -4.24785e-03, -1.72126e-03, -6.98913e-04, 2.92923e-03, 4.87294e-03, 6.46827e-03, 7.66794e-03, 8.43444e-03, 8.77681e-03, -8.67374e-03, 8.20642e-03, 7.31013e-03, 6.23237e-03, 4.87796e-03, 3.31983e-03, 1.71169e-03, -1.00410e-04, -1.43003e-03, -2.82224e-03, -4.01796e-03, -4.97525e-03, -5.65829e-03, -6.05118e-03, --6.14483e-03, -5.95049e-03, -5.48593e-03, -4.78659e-03, -3.88827e-03, -2.84712e-03, -1.70981e-03, --5.36994e-04, 6.20542e-04, 1.70418e-03, 2.67306e-03, 3.47811e-03, 4.09614e-03, 4.50037e-03, -4.68091e-03, 4.63546e-03, 4.38017e-03, 3.91841e-03, 3.30312e-03, 2.53826e-03, 1.69402e-03, -7.83910e-04, -1.22305e-04, -1.00532e-03, -1.80630e-03, -2.50229e-03, -3.06572e-03, -3.44913e-03, --3.69208e-03, -3.73609e-03, -3.60506e-03, -3.31209e-03, -2.87636e-03, -2.31717e-03, -1.66939e-03, --9.57591e-04, -2.23591e-04, 5.04003e-04, 1.18587e-03, 1.79763e-03, 2.30599e-03, 2.69770e-03, -2.95067e-03, 3.06271e-03, 3.02785e-03, 2.85589e-03, 2.55400e-03, 2.14503e-03, 1.64530e-03, -1.08328e-03, 4.83736e-04, -1.24964e-04, -7.15425e-04, -1.24542e-03, -1.72140e-03, -2.10060e-03, --2.37032e-03, -2.53448e-03, -2.56546e-03, -2.49179e-03, -2.29070e-03, -2.00403e-03, -1.61659e-03, --1.17826e-03, -6.83619e-04, -1.74591e-04, 3.29944e-04, 8.13976e-04, 1.24507e-03, 1.61153e-03, -1.89772e-03, 2.08906e-03, 2.18121e-03, 2.17177e-03, 2.06312e-03, 1.86382e-03, 1.58411e-03, -1.23827e-03, 8.43172e-04, 4.19904e-04, -1.56907e-05, -4.39507e-04, -8.37003e-04, -1.18696e-03, --1.47938e-03, -1.69770e-03, -1.83969e-03, -1.89318e-03, -1.86141e-03, -1.74292e-03, -1.54080e-03, --1.27846e-03, -9.95436e-04, -6.26595e-04, -2.64150e-04, 1.13622e-04, 4.76814e-04, 8.10335e-04, -1.10736e-03, 1.34273e-03, 1.52430e-03, 1.62390e-03, 1.66246e-03, 1.61522e-03, 1.50632e-03, -1.32668e-03, 1.09739e-03, 8.19243e-04, 5.11244e-04, 1.80867e-04, -1.46339e-04, -4.68008e-04, --7.60099e-04, -1.02128e-03, -1.23115e-03, -1.39047e-03, -1.48008e-03, -1.51120e-03, -1.46991e-03, --1.37032e-03, -1.20407e-03, -9.91551e-04, -7.27194e-04, -4.35173e-04, -1.11623e-04, 2.15793e-04, -5.54161e-04, 8.66979e-04, 1.18292e-03, 1.42542e-03, 1.77484e-03, 1.69052e-03, 2.12100e-03, -1.90833e-03, 2.18455e-03, 1.90159e-03, 2.03277e-03, 1.68588e-03, 1.69365e-03, 1.30556e-03, -1.21725e-03, 8.19039e-04, 6.67167e-04, 2.89567e-04, 1.11906e-04, -2.18096e-04, -3.95376e-04, --6.46202e-04, -8.02450e-04, -9.64235e-04, -1.08660e-03, -1.15302e-03, -1.23904e-03, -1.20955e-03, --1.26937e-03, -1.15443e-03, -1.19836e-03, -1.01377e-03, -1.05669e-03, -8.15327e-04, -8.79730e-04, --5.93148e-04, -6.95337e-04, -3.75376e-04, -5.27511e-04, -1.78544e-04, -3.96418e-04, -8.80292e-06, --3.37279e-04 }; +float convolve(const float *in, const float *taps, size_t len); +float complex hilbert_transform(const float *in, const float *taps, size_t len); +float interpolating_convolve(const float *in, const float *taps, size_t len, float offset, float delta); diff --git a/src/image.c b/src/image.c index 875265f..274373b 100644 --- a/src/image.c +++ b/src/image.c @@ -33,8 +33,8 @@ static linear_t compute_regression(float *wedges) { return linear_regression(wedges, teleramp, 9); } -static double tele[16]; -static double Cs; +static float tele[16]; +static float Cs; void apt_histogramEqualise(float **prow, int nrow, int offset, int width){ // Plot histogram @@ -97,9 +97,9 @@ void calibrateImage(float **prow, int nrow, int offset, int width, linear_t regr } } -double teleNoise(float *wedges){ - double pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; - double noise = 0; +float teleNoise(float *wedges){ + float pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; + float noise = 0; for(int i = 0; i < 9; i++) noise += fabs(wedges[i] - pattern[i]); @@ -108,7 +108,7 @@ double teleNoise(float *wedges){ // Get telemetry data for thermal calibration apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { - double teleline[APT_MAX_HEIGHT] = { 0.0 }; + float teleline[APT_MAX_HEIGHT] = { 0.0 }; float wedge[16]; linear_t regr[APT_MAX_HEIGHT/APT_FRAME_LEN + 1]; int telestart, mtelestart = 0; @@ -156,7 +156,7 @@ apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { } // Find the least noisy frame - double minNoise = -1; + float minNoise = -1; int bestFrame = -1; for (int n = telestart, k = 0; n < nrow - APT_FRAME_LEN; n += APT_FRAME_LEN, k++) { int j; @@ -170,7 +170,7 @@ apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { wedge[j] /= 6; } - double noise = teleNoise(wedge); + float noise = teleNoise(wedge); if(noise < minNoise || minNoise == -1){ minNoise = noise; bestFrame = k; @@ -313,16 +313,16 @@ int apt_cropNoise(apt_image_t *img){ #include "satcal.h" typedef struct { - double Nbb; - double Cs; - double Cb; + float Nbb; + float Cs; + float Cb; int ch; } tempparam_t; // IR channel temperature compensation -static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) { - double Tbb, T[4]; - double C; +static void tempcomp(float t[16], int ch, int satnum, tempparam_t *tpr) { + float Tbb, T[4]; + float C; tpr->ch = ch - 4; @@ -352,9 +352,9 @@ static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) { } // IR channel temperature calibration -static double tempcal(float Ce, int satnum, tempparam_t * rgpr) { - double Nl, Nc, Ns, Ne; - double T, vc; +static float tempcal(float Ce, int satnum, tempparam_t * rgpr) { + float Nl, Nc, Ns, Ne; + float T, vc; Ns = satcal[satnum].cor[rgpr->ch].Ns; Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb); diff --git a/src/taps.h b/src/taps.h new file mode 100644 index 0000000..9d5f5cf --- /dev/null +++ b/src/taps.h @@ -0,0 +1,92 @@ +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2004-2009 Thierry Leconte (F4DWV) 2019-2022 Xerbo (xerbo@protonmail.com) + * + * This program 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 2 of the License, or + * (at your option) any later version. + * + * This program 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 this program. If not, see . + */ + +static const float hilbert_filter[] = { 0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, +0.0335063, 0.0374482, 0.0424413, 0.0489708, 0.0578745, 0.0707355, 0.0909457, 0.127324, 0.212207, 0.63662, +-0.63662, -0.212207, -0.127324, -0.0909457, -0.0707355, -0.0578745, -0.0489708, -0.0424413, -0.0374482, +-0.0335063, -0.0303152, -0.0276791, -0.0254648, -0.0235785, -0.0219524, -0.0205361 }; +#define HILBERT_FILTER_SIZE (sizeof(hilbert_filter)/sizeof(hilbert_filter[0])) + +static const float low_pass[] = { -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, +-3.75376e-04, -6.95337e-04, -5.93148e-04, -8.79730e-04, -8.15327e-04, -1.05669e-03, -1.01377e-03, +-1.19836e-03, -1.15443e-03, -1.26937e-03, -1.20955e-03, -1.23904e-03, -1.15302e-03, -1.08660e-03, +-9.64235e-04, -8.02450e-04, -6.46202e-04, -3.95376e-04, -2.18096e-04, 1.11906e-04, 2.89567e-04, +6.67167e-04, 8.19039e-04, 1.21725e-03, 1.30556e-03, 1.69365e-03, 1.68588e-03, 2.03277e-03, 1.90159e-03, +2.18455e-03, 1.90833e-03, 2.12100e-03, 1.69052e-03, 1.77484e-03, 1.42542e-03, 1.18292e-03, 8.66979e-04, +5.54161e-04, 2.15793e-04, -1.11623e-04, -4.35173e-04, -7.27194e-04, -9.91551e-04, -1.20407e-03, +-1.37032e-03, -1.46991e-03, -1.51120e-03, -1.48008e-03, -1.39047e-03, -1.23115e-03, -1.02128e-03, +-7.60099e-04, -4.68008e-04, -1.46339e-04, 1.80867e-04, 5.11244e-04, 8.19243e-04, 1.09739e-03, 1.32668e-03, +1.50632e-03, 1.61522e-03, 1.66246e-03, 1.62390e-03, 1.52430e-03, 1.34273e-03, 1.10736e-03, 8.10335e-04, +4.76814e-04, 1.13622e-04, -2.64150e-04, -6.26595e-04, -9.95436e-04, -1.27846e-03, -1.54080e-03, +-1.74292e-03, -1.86141e-03, -1.89318e-03, -1.83969e-03, -1.69770e-03, -1.47938e-03, -1.18696e-03, +-8.37003e-04, -4.39507e-04, -1.56907e-05, 4.19904e-04, 8.43172e-04, 1.23827e-03, 1.58411e-03, +1.86382e-03, 2.06312e-03, 2.17177e-03, 2.18121e-03, 2.08906e-03, 1.89772e-03, 1.61153e-03, +1.24507e-03, 8.13976e-04, 3.29944e-04, -1.74591e-04, -6.83619e-04, -1.17826e-03, -1.61659e-03, +-2.00403e-03, -2.29070e-03, -2.49179e-03, -2.56546e-03, -2.53448e-03, -2.37032e-03, -2.10060e-03, +-1.72140e-03, -1.24542e-03, -7.15425e-04, -1.24964e-04, 4.83736e-04, 1.08328e-03, 1.64530e-03, +2.14503e-03, 2.55400e-03, 2.85589e-03, 3.02785e-03, 3.06271e-03, 2.95067e-03, 2.69770e-03, +2.30599e-03, 1.79763e-03, 1.18587e-03, 5.04003e-04, -2.23591e-04, -9.57591e-04, -1.66939e-03, +-2.31717e-03, -2.87636e-03, -3.31209e-03, -3.60506e-03, -3.73609e-03, -3.69208e-03, -3.44913e-03, +-3.06572e-03, -2.50229e-03, -1.80630e-03, -1.00532e-03, -1.22305e-04, 7.83910e-04, 1.69402e-03, +2.53826e-03, 3.30312e-03, 3.91841e-03, 4.38017e-03, 4.63546e-03, 4.68091e-03, 4.50037e-03, +4.09614e-03, 3.47811e-03, 2.67306e-03, 1.70418e-03, 6.20542e-04, -5.36994e-04, -1.70981e-03, +-2.84712e-03, -3.88827e-03, -4.78659e-03, -5.48593e-03, -5.95049e-03, -6.14483e-03, -6.05118e-03, +-5.65829e-03, -4.97525e-03, -4.01796e-03, -2.82224e-03, -1.43003e-03, 1.00410e-04, 1.71169e-03, +3.31983e-03, 4.87796e-03, 6.23237e-03, 7.31013e-03, 8.20642e-03, 8.67374e-03, 8.77681e-03, +8.43444e-03, 7.66794e-03, 6.46827e-03, 4.87294e-03, 2.92923e-03, 6.98913e-04, -1.72126e-03, +-4.24785e-03, -6.75380e-03, -9.13309e-03, -1.12532e-02, -1.30038e-02, -1.42633e-02, -1.49338e-02, +-1.49145e-02, -1.41484e-02, -1.25761e-02, -1.01870e-02, -6.97432e-03, -2.97910e-03, 1.75386e-03, +7.11899e-03, 1.30225e-02, 1.93173e-02, 2.58685e-02, 3.24965e-02, 3.90469e-02, 4.53316e-02, +5.11931e-02, 5.64604e-02, 6.09924e-02, 6.46584e-02, 6.73547e-02, 6.90049e-02, 6.97096e-02, +6.90049e-02, 6.73547e-02, 6.46584e-02, 6.09924e-02, 5.64604e-02, 5.11931e-02, 4.53316e-02, +3.90469e-02, 3.24965e-02, 2.58685e-02, 1.93173e-02, 1.30225e-02, 7.11899e-03, 1.75386e-03, +-2.97910e-03, -6.97432e-03, -1.01870e-02, -1.25761e-02, -1.41484e-02, -1.49145e-02, -1.49338e-02, +-1.42633e-02, -1.30038e-02, -1.12532e-02, -9.13309e-03, -6.75380e-03, -4.24785e-03, -1.72126e-03, +6.98913e-04, 2.92923e-03, 4.87294e-03, 6.46827e-03, 7.66794e-03, 8.43444e-03, 8.77681e-03, +8.67374e-03, 8.20642e-03, 7.31013e-03, 6.23237e-03, 4.87796e-03, 3.31983e-03, 1.71169e-03, +1.00410e-04, -1.43003e-03, -2.82224e-03, -4.01796e-03, -4.97525e-03, -5.65829e-03, -6.05118e-03, +-6.14483e-03, -5.95049e-03, -5.48593e-03, -4.78659e-03, -3.88827e-03, -2.84712e-03, -1.70981e-03, +-5.36994e-04, 6.20542e-04, 1.70418e-03, 2.67306e-03, 3.47811e-03, 4.09614e-03, 4.50037e-03, +4.68091e-03, 4.63546e-03, 4.38017e-03, 3.91841e-03, 3.30312e-03, 2.53826e-03, 1.69402e-03, +7.83910e-04, -1.22305e-04, -1.00532e-03, -1.80630e-03, -2.50229e-03, -3.06572e-03, -3.44913e-03, +-3.69208e-03, -3.73609e-03, -3.60506e-03, -3.31209e-03, -2.87636e-03, -2.31717e-03, -1.66939e-03, +-9.57591e-04, -2.23591e-04, 5.04003e-04, 1.18587e-03, 1.79763e-03, 2.30599e-03, 2.69770e-03, +2.95067e-03, 3.06271e-03, 3.02785e-03, 2.85589e-03, 2.55400e-03, 2.14503e-03, 1.64530e-03, +1.08328e-03, 4.83736e-04, -1.24964e-04, -7.15425e-04, -1.24542e-03, -1.72140e-03, -2.10060e-03, +-2.37032e-03, -2.53448e-03, -2.56546e-03, -2.49179e-03, -2.29070e-03, -2.00403e-03, -1.61659e-03, +-1.17826e-03, -6.83619e-04, -1.74591e-04, 3.29944e-04, 8.13976e-04, 1.24507e-03, 1.61153e-03, +1.89772e-03, 2.08906e-03, 2.18121e-03, 2.17177e-03, 2.06312e-03, 1.86382e-03, 1.58411e-03, +1.23827e-03, 8.43172e-04, 4.19904e-04, -1.56907e-05, -4.39507e-04, -8.37003e-04, -1.18696e-03, +-1.47938e-03, -1.69770e-03, -1.83969e-03, -1.89318e-03, -1.86141e-03, -1.74292e-03, -1.54080e-03, +-1.27846e-03, -9.95436e-04, -6.26595e-04, -2.64150e-04, 1.13622e-04, 4.76814e-04, 8.10335e-04, +1.10736e-03, 1.34273e-03, 1.52430e-03, 1.62390e-03, 1.66246e-03, 1.61522e-03, 1.50632e-03, +1.32668e-03, 1.09739e-03, 8.19243e-04, 5.11244e-04, 1.80867e-04, -1.46339e-04, -4.68008e-04, +-7.60099e-04, -1.02128e-03, -1.23115e-03, -1.39047e-03, -1.48008e-03, -1.51120e-03, -1.46991e-03, +-1.37032e-03, -1.20407e-03, -9.91551e-04, -7.27194e-04, -4.35173e-04, -1.11623e-04, 2.15793e-04, +5.54161e-04, 8.66979e-04, 1.18292e-03, 1.42542e-03, 1.77484e-03, 1.69052e-03, 2.12100e-03, +1.90833e-03, 2.18455e-03, 1.90159e-03, 2.03277e-03, 1.68588e-03, 1.69365e-03, 1.30556e-03, +1.21725e-03, 8.19039e-04, 6.67167e-04, 2.89567e-04, 1.11906e-04, -2.18096e-04, -3.95376e-04, +-6.46202e-04, -8.02450e-04, -9.64235e-04, -1.08660e-03, -1.15302e-03, -1.23904e-03, -1.20955e-03, +-1.26937e-03, -1.15443e-03, -1.19836e-03, -1.01377e-03, -1.05669e-03, -8.15327e-04, -8.79730e-04, +-5.93148e-04, -6.95337e-04, -3.75376e-04, -5.27511e-04, -1.78544e-04, -3.96418e-04, -8.80292e-06, +-3.37279e-04 }; +#define LOW_PASS_SIZE (sizeof(low_pass)/sizeof(low_pass[0])) + +static const float sync_pattern[] = { -14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, +18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14 }; +#define SYNC_PATTERN_SIZE (sizeof(sync_pattern)/sizeof(sync_pattern[0])) diff --git a/src/util.h b/src/util.h new file mode 100644 index 0000000..edcfeb9 --- /dev/null +++ b/src/util.h @@ -0,0 +1,20 @@ +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2019-2022 Xerbo (xerbo@protonmail.com) + * + * This program 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 2 of the License, or + * (at your option) any later version. + * + * This program 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 this program. If not, see . + */ + +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define MAX(a,b) (((a) > (b)) ? (a) : (b))