WDSP: Sonar corrections and cleanup

wdsp/CMakeLists.txt

@@ -134,14 +134,3 @@ if (MSVC)
 endif()
 
 install(TARGETS wdsp DESTINATION ${INSTALL_LIB_DIR})
-
-if (LINUX)
-    add_executable(wdsp_make_interface
-        make_interface.cpp
-    )
-
-    add_executable(wdsp_make_calculus
-        make_calculus.cpp
-    )
-    install(TARGETS wdsp_make_interface wdsp_make_calculus DESTINATION ${INSTALL_BIN_DIR})
-endif()

wdsp/cfcomp.cpp

@@ -450,9 +450,9 @@ void CFCOMP::SetCFCOMPprofile (TXA& txa, int nfreqs, float* F, float* G, float *
     delete[] (a->ep);
     delete[] (a->gp);
     delete[] (a->fp);
-    a->fp = new float[a->nfreqs]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
+    a->fp = new float[a->nfreqs + 2]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
-    a->gp = new float[a->nfreqs]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
+    a->gp = new float[a->nfreqs + 2]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
-    a->ep = new float[a->nfreqs]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
+    a->ep = new float[a->nfreqs + 2]; // (float *) malloc0 ((a->nfreqs + 2) * sizeof (float));
     calc_comp(a);
 }
 

wdsp/cfir.cpp

@@ -46,8 +46,23 @@ void CFIR::decalc_cfir (CFIR *a)
     FIRCORE::destroy_fircore (a->p);
 }
 
-CFIR* CFIR::create_cfir (int run, int size, int nc, int mp, float* in, float* out, int runrate, int cicrate,
+CFIR* CFIR::create_cfir (
-    int DD, int R, int Pairs, float cutoff, int xtype, float xbw, int wintype)
+    int run,
+    int size,
+    int nc,
+    int mp,
+    float* in,
+    float* out,
+    int runrate,
+    int cicrate,
+    int DD,
+    int R,
+    int Pairs,
+    double cutoff,
+    int xtype,
+    double xbw,
+    int wintype
+)
 //  run:  0 - no action; 1 - operate
 //  size:  number of complex samples in an input buffer to the CFIR filter
 //  nc:  number of filter coefficients
@@ -130,7 +145,20 @@ void CFIR::setOutRate_cfir (CFIR *a, int rate)
     calc_cfir (a);
 }
 
-float* CFIR::cfir_impulse (int N, int DD, int R, int Pairs, float runrate, float cicrate, float cutoff, int xtype, float xbw, int rtype, float scale, int wintype)
+float* CFIR::cfir_impulse (
+    int N,
+    int DD,
+    int R,
+    int Pairs,
+    double runrate,
+    double cicrate,
+    double cutoff,
+    int xtype,
+    double xbw,
+    int rtype,
+    double scale,
+    int wintype
+)
 {
     // N:       number of impulse response samples
     // DD:      differential delay used in the CIC filter
@@ -144,18 +172,18 @@ float* CFIR::cfir_impulse (int N, int DD, int R, int Pairs, float runrate, float
     // rtype:   0 for real output, 1 for complex output
     // scale:   scale factor to be applied to the output
     int i, j;
-    float tmp, local_scale, ri, mag, fn;
+    double tmp, local_scale, ri, mag, fn;
     float* impulse;
     float* A = new float[N]; // (float *) malloc0 (N * sizeof (float));
-    float ft = cutoff / cicrate;                                       // normalized cutoff frequency
+    double ft = cutoff / cicrate;                                       // normalized cutoff frequency
     int u_samps = (N + 1) / 2;                                          // number of unique samples,  OK for odd or even N
     int c_samps = (int)(cutoff / runrate * N) + (N + 1) / 2 - N / 2;    // number of unique samples within bandpass, OK for odd or even N
     int x_samps = (int)(xbw / runrate * N);                             // number of unique samples in transition region, OK for odd or even N
-    float offset = 0.5 - 0.5 * (float)((N + 1) / 2 - N / 2);          // sample offset from center, OK for odd or even N
+    double offset = 0.5 - 0.5 * (float)((N + 1) / 2 - N / 2);          // sample offset from center, OK for odd or even N
-    float* xistion = new float[x_samps + 1]; // (float *) malloc0 ((x_samps + 1) * sizeof (float));
+    double* xistion = new double[x_samps + 1]; // (float *) malloc0 ((x_samps + 1) * sizeof (float));
-    float delta = PI / (float)x_samps;
+    double delta = PI / (float)x_samps;
-    float L = cicrate / runrate;
+    double L = cicrate / runrate;
-    float phs = 0.0;
+    double phs = 0.0;
     for (i = 0; i <= x_samps; i++)
     {
         xistion[i] = 0.5 * (cos (phs) + 1.0);
@@ -171,7 +199,8 @@ float* CFIR::cfir_impulse (int N, int DD, int R, int Pairs, float runrate, float
             fn = ri / (L * (float)N);
             if (fn <= ft)
             {
-                if (fn == 0.0) tmp = 1.0;
+                if (fn == 0.0)
+                    tmp = 1.0;
                 else if ((tmp = DD * R * sin (PI * fn / R) / sin (PI * DD * fn)) < 0.0)
                     tmp = -tmp;
                 mag = pow (tmp, Pairs) * local_scale;
@@ -225,7 +254,8 @@ float* CFIR::cfir_impulse (int N, int DD, int R, int Pairs, float runrate, float
             A[i] = A[u_samps - j];
     impulse = FIR::fir_fsamp (N, A, rtype, 1.0, wintype);
     // print_impulse ("cfirImpulse.txt", N, impulse, 1, 0);
-    delete[] (A);
+    delete[] A;
+    delete[] xistion;
     return impulse;
 }
 

wdsp/cfir.hpp

@@ -49,10 +49,10 @@ public:
     int DD;
     int R;
     int Pairs;
-    float cutoff;
+    double cutoff;
-    float scale;
+    double scale;
     int xtype;
-    float xbw;
+    double xbw;
     int wintype;
     FIRCORE *p;
 
@@ -68,9 +68,9 @@ public:
         int DD,
         int R,
         int Pairs,
-        float cutoff,
+        double cutoff,
         int xtype,
-        float xbw,
+        double xbw,
         int wintype
     );
     static void destroy_cfir (CFIR *a);
@@ -85,13 +85,13 @@ public:
         int DD,
         int R,
         int Pairs,
-        float runrate,
+        double runrate,
-        float cicrate,
+        double cicrate,
-        float cutoff,
+        double cutoff,
         int xtype,
-        float xbw,
+        double xbw,
         int rtype,
-        float scale,
+        double scale,
         int wintype
     );
     // TXA Properties

wdsp/channel.hpp

@@ -1,90 +0,0 @@
-/*  channel.h
-
-This file is part of a program that implements a Software-Defined Radio.
-
-Copyright (C) 2013 Warren Pratt, NR0V
-Copyright (C) 2024 Edouard Griffiths, F4EXB Adapted to SDRangel
-
-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, write to the Free Software
-Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
-
-The author can be reached by email at
-
-warren@wpratt.com
-
-*/
-
-#ifndef wdsp_channel_h
-#define wdsp_channel_h
-
-#include <QThread>
-
-#include "export.h"
-
-class WDSP_API Channel
-{
-public:
-    int type;
-    bool run;                   // thread running
-    int in_rate;                // input samplerate
-    int out_rate;               // output samplerate
-    int in_size;                // input buffsize (complex samples) in a fexchange() operation
-    int dsp_rate;               // sample rate for mainstream dsp processing
-    int dsp_size;               // number complex samples processed per buffer in mainstream dsp processing
-    int dsp_insize;             // size (complex samples) of the output of the r1 (input) buffer
-    int dsp_outsize;            // size (complex samples) of the input of the r2 (output) buffer
-    int out_size;               // output buffsize (complex samples) in a fexchange() operation
-    int state;                  // 0 for channel OFF; 1 for channel ON
-    float tdelayup;
-    float tslewup;
-    float tdelaydown;
-    float tslewdown;
-    int bfo;                    // 'block_for_output', block fexchange until output is available
-    volatile long flushflag;
-    QThread channelThread;
-    long upslew;
-    // struct  //io buffers
-    // {
-    //     IOB pc, pd, pe, pf;     // copies for console calls, dsp, exchange, and flush thread
-    //     volatile long ch_upslew;
-    // } iob;
-
-    static void create_channel (
-        int channel,
-        int in_size,
-        int dsp_size,
-        int input_samplerate,
-        int dsp_rate,
-        int output_samplerate,
-        int type,
-        int state,
-        float tdelayup,
-        float tslewup,
-        float tdelaydown,
-        float tslewdown,
-        int bfo
-    );
-    static void destroy_channel (int channel);
-    static void flush_channel (int channel);
-    // static void set_type (int channel, int type);
-    // static void SetInputBuffsize (int channel, int in_size);
-    // static void SetDSPBuffsize (int channel, int dsp_size);
-    // static void SetInputSamplerate  (int channel, int samplerate);
-    // static void SetDSPSamplerate (int channel, int samplerate);
-    // static void SetOutputSamplerate (int channel, int samplerate);
-    // static void SetAllRates (int channel, int in_rate, int dsp_rate, int out_rate);
-    // static int SetChannelState (int channel, int state, int dmode);
-};
-
-#endif

wdsp/emnr.cpp

@@ -233,7 +233,7 @@ void EMNR::interpM (double* res, double x, int nvals, double* xvals, double* yva
         int idx = 0;
         double xllow, xlhigh, frac;
 
-        while (x >= xvals[idx])
+        while ((x >= xvals[idx]) && (idx < nvals - 1))
             idx++;
 
         xllow = log10 (xvals[idx - 1]);

wdsp/eq.cpp

@@ -44,46 +44,60 @@ int EQP::fEQcompare (const void * a, const void * b)
         return 1;
 }
 
-float* EQP::eq_impulse (int N, int nfreqs, float* F, float* G, float samplerate, float scale, int ctfmode, int wintype)
+float* EQP::eq_impulse (int N, int nfreqs, float* F, float* G, double samplerate, double scale, int ctfmode, int wintype)
 {
     float* fp = new float[nfreqs + 2]; // (float *) malloc0 ((nfreqs + 2)   * sizeof (float));
     float* gp = new float[nfreqs + 2]; // (float *) malloc0 ((nfreqs + 2)   * sizeof (float));
     float* A  = new float[N / 2 + 1]; // (float *) malloc0 ((N / 2 + 1) * sizeof (float));
     float* sary = new float[2 * nfreqs]; // (float *) malloc0 (2 * nfreqs * sizeof (float));
-    float gpreamp, f, frac;
+    double gpreamp, f, frac;
     float* impulse;
     int i, j, mid;
     fp[0] = 0.0;
     fp[nfreqs + 1] = 1.0;
     gpreamp = G[0];
+
     for (i = 1; i <= nfreqs; i++)
     {
         fp[i] = 2.0 * F[i] / samplerate;
-        if (fp[i] < 0.0) fp[i] = 0.0;
+
-        if (fp[i] > 1.0) fp[i] = 1.0;
+        if (fp[i] < 0.0)
+            fp[i] = 0.0;
+
+        if (fp[i] > 1.0)
+            fp[i] = 1.0;
+
         gp[i] = G[i];
     }
+
     for (i = 1, j = 0; i <= nfreqs; i++, j+=2)
     {
         sary[j + 0] = fp[i];
         sary[j + 1] = gp[i];
     }
+
     qsort (sary, nfreqs, 2 * sizeof (float), fEQcompare);
+
     for (i = 1, j = 0; i <= nfreqs; i++, j+=2)
     {
         fp[i] = sary[j + 0];
         gp[i] = sary[j + 1];
     }
+
     gp[0] = gp[1];
     gp[nfreqs + 1] = gp[nfreqs];
     mid = N / 2;
     j = 0;
+
     if (N & 1)
     {
         for (i = 0; i <= mid; i++)
         {
-            f = (float)i / (float)mid;
+            f = (double)i / (double)mid;
-            while (f > fp[j + 1]) j++;
+
+            while ((f > fp[j + 1]) && (j < nfreqs))
+                j++;
+
             frac = (f - fp[j]) / (fp[j + 1] - fp[j]);
             A[i] = pow (10.0, 0.05 * (frac * gp[j + 1] + (1.0 - frac) * gp[j] + gpreamp)) * scale;
         }
@@ -92,36 +106,44 @@ float* EQP::eq_impulse (int N, int nfreqs, float* F, float* G, float samplerate,
     {
         for (i = 0; i < mid; i++)
         {
-            f = ((float)i + 0.5) / (float)mid;
+            f = ((double)i + 0.5) / (double)mid;
-            while (f > fp[j + 1]) j++;
+
+            while ((f > fp[j + 1]) && (j < nfreqs))
+                j++;
+
             frac = (f - fp[j]) / (fp[j + 1] - fp[j]);
             A[i] = pow (10.0, 0.05 * (frac * gp[j + 1] + (1.0 - frac) * gp[j] + gpreamp)) * scale;
         }
     }
+
     if (ctfmode == 0)
     {
         int k, low, high;
-        float lowmag, highmag, flow4, fhigh4;
+        double lowmag, highmag, flow4, fhigh4;
+
         if (N & 1)
         {
             low = (int)(fp[1] * mid);
             high = (int)(fp[nfreqs] * mid + 0.5);
             lowmag = A[low];
             highmag = A[high];
-            flow4 = pow((float)low / (float)mid, 4.0);
+            flow4 = pow((double)low / (double)mid, 4.0);
-            fhigh4 = pow((float)high / (float)mid, 4.0);
+            fhigh4 = pow((double)high / (double)mid, 4.0);
             k = low;
+
             while (--k >= 0)
             {
-                f = (float)k / (float)mid;
+                f = (double)k / (double)mid;
                 lowmag *= (f * f * f * f) / flow4;
                 if (lowmag < 1.0e-20) lowmag = 1.0e-20;
                 A[k] = lowmag;
             }
+
             k = high;
+
             while (++k <= mid)
             {
-                f = (float)k / (float)mid;
+                f = (double)k / (double)mid;
                 highmag *= fhigh4 / (f * f * f * f);
                 if (highmag < 1.0e-20) highmag = 1.0e-20;
                 A[k] = highmag;
@@ -133,30 +155,35 @@ float* EQP::eq_impulse (int N, int nfreqs, float* F, float* G, float samplerate,
             high = (int)(fp[nfreqs] * mid - 0.5);
             lowmag = A[low];
             highmag = A[high];
-            flow4 = pow((float)low / (float)mid, 4.0);
+            flow4 = pow((double)low / (double)mid, 4.0);
-            fhigh4 = pow((float)high / (float)mid, 4.0);
+            fhigh4 = pow((double)high / (double)mid, 4.0);
             k = low;
+
             while (--k >= 0)
             {
-                f = (float)k / (float)mid;
+                f = (double)k / (double)mid;
                 lowmag *= (f * f * f * f) / flow4;
                 if (lowmag < 1.0e-20) lowmag = 1.0e-20;
                 A[k] = lowmag;
             }
+
             k = high;
+
             while (++k < mid)
             {
-                f = (float)k / (float)mid;
+                f = (double)k / (double)mid;
                 highmag *= fhigh4 / (f * f * f * f);
                 if (highmag < 1.0e-20) highmag = 1.0e-20;
                 A[k] = highmag;
             }
         }
     }
+
     if (N & 1)
         impulse = FIR::fir_fsamp_odd(N, A, 1, 1.0, wintype);
     else
         impulse = FIR::fir_fsamp(N, A, 1, 1.0, wintype);
+
     // print_impulse("eq.txt", N, impulse, 1, 0);
     delete[] (sary);
     delete[] (A);

wdsp/eq.hpp

@@ -56,7 +56,7 @@ public:
     float* G;
     int ctfmode;
     int wintype;
-    float samplerate;
+    double samplerate;
     FIRCORE *p;
 
     static EQP* create_eqp (
@@ -73,7 +73,7 @@ public:
         int wintype,
         int samplerate
     );
-    static float* eq_impulse (int N, int nfreqs, float* F, float* G, float samplerate, float scale, int ctfmode, int wintype);
+    static float* eq_impulse (int N, int nfreqs, float* F, float* G, double samplerate, double scale, int ctfmode, int wintype);
     static void destroy_eqp (EQP *a);
     static void flush_eqp (EQP *a);
     static void xeqp (EQP *a);

wdsp/fir.cpp

@@ -318,8 +318,8 @@ float *FIR::fir_read (int N, const char *filename, int rtype, float scale)
 void FIR::analytic (int N, float* in, float* out)
 {
     int i;
-    float inv_N = 1.0 / (float) N;
+    double inv_N = 1.0 / (double) N;
-    float two_inv_N = 2.0 * inv_N;
+    double two_inv_N = 2.0 * inv_N;
     float* x = new float[N * 2]; // (float *) malloc0 (N * sizeof (complex));
     fftwf_plan pfor = fftwf_plan_dft_1d (
         N,

wdsp/make_calculus.cpp

@@ -1,63 +0,0 @@
-/*
- * make_calculus
- *
- * This program reads the contents of the binary WDSP file "calculus"
- * and dumps the data as two arrays of floating-point numbers
- *
- * The output is intended to be part of the file "calculus.c" which
- * initializes these arrays (static data) for use with "memcpy"
- * in emnr.c.
- *
- * Should the WDSP file "calculus" be changed, "calculus.c" should
- * be re-generated using this program.
- *
- * return values of main()
- *
- *  0  all OK
- * -1  sizeof(float) is not 8
- * -2  error opening file "calculus"
- * -3  read error
- */
-
-#include <stdio.h>
-#include <unistd.h>
-#include <fcntl.h>
-
-int main() {
-    int fd;
-    int i,j;
-    float d;
-
-    if (sizeof(float) != 8) {
-        printf("Data type DOUBLE is not 8-byte. Please check!\n");
-        return -1;
-    }
-    fd=open ("calculus", O_RDONLY);
-    if (fd < 0) {
-        printf("Could not open file 'calculus'\n");
-        return -2;
-    }
-
-    for (j=0; j<2; j++) {
-        switch (j) {
-            case 0:
-            printf("float GG[241*241]={\n");
-            break;
-            case 1:
-            printf("float GGS[241*241]={\n");
-            break;
-        }
-        for (i=0; i< 241*241; i++) {
-            if (read(fd, &d, 8) != 8) {
-                printf("READ ERROR\n");
-                return -3;
-            }
-            if (i == 241*241 -1) {
-                printf("%30.25f};\n", d);
-            } else {
-                printf("%30.25f,\n", d);
-            }
-        }
-        return 0;
-    }
-}

wdsp/make_interface.cpp

@@ -1,101 +0,0 @@
-/*
-    The purpose of this file is to extract interfaces from the WDSP source code.
-    The interfaces have the following form:
-
-    PORT blabla
-    firstline
-    secondline
-    {
-
-    where there may be an arbitrary number of lines between the line
-    containing "PORT" and the line starting with "{". This has to be
-    converted to
-
-    extern blabla firstline
-    secondline;
-
-    That is, the first line is prepended by "extern", and the last line is closed
-    with a semicolon. Comments starting with '//' are omitted, and lines starting
-    with '//' are ignored.
-*/
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-void trimm(char *line, size_t maxlen);
-
-int main(int argc, char **argv)
-{
-    FILE *infile;
-    int i;
-    int first_in_file;
-    int first_in_decl;
-    char line[1000];
-    size_t linesize=999;
-    char *buffer=line;
-
-    for (i=1; i<argc; i++) {
-        infile=fopen(argv[i],"r");
-        if (infile == NULL) continue;
-        first_in_file=1;
-        for (;;) {
-            if (getline(&buffer, &linesize, infile) < 0) break;
-            trimm(line, linesize);
-            if (strncmp(line,"PORT", 4) != 0) continue;
-            // found an interface
-            if (first_in_file) {
-                printf("\n//\n// Interfaces from %s\n//\n\n", argv[i]);
-                first_in_file=0;
-            }
-            if (strlen(line) >4) {
-                printf("extern %s ", line+4);
-            } else {
-                printf("extern ");
-            }
-            first_in_decl=1;
-
-            for (;;) {
-                if (getline(&buffer, &linesize, infile) < 0) {
-                    fprintf(stderr,"! Found a PORT but found EOF while scanning interface.\n");
-                    return 8;
-                }
-                trimm(line, linesize);
-                if (line[0] == 0) continue;
-                if (line[0] == '{') {
-                    printf(";\n");
-                    break;
-                } else {
-                    if (first_in_decl) {
-                        printf("%s", line);
-                        first_in_decl=0;
-                    } else {
-                        printf("\n%s", line);
-                    }
-                }
-            }
-        }
-        fclose(infile);
-    }
-    return 0;
-}
-
-void trimm(char *line, size_t maxlen) {
-    int len;
-
-    //
-    // Remove comments starting with '//'
-    //
-    len=strnlen(line,maxlen);
-    for (int i=0; i< len-1; i++) {
-        if (line[i] == '/' && line[i+1] == '/') line[i]=0;
-    }
-
-    //
-    // Remove trailing white space and newlines
-    //
-    len=strnlen(line,maxlen);
-    line[len--]=0;
-    while (len >= 0 && (line[len] == ' ' || line[len] == '\t' || line[len]== '\n')) line[len--]=0;
-}

wdsp/nbp.cpp

@@ -132,6 +132,7 @@ int NBP::make_nbp (
     else
     {
         nbp = 0;
+        delete[] del;
         return nbp;
     }
     *havnotch = 0;

wdsp/resample.cpp

@@ -60,6 +60,9 @@ void RESAMPLE::calc_resample (RESAMPLE *a)
     a->L = a->out_rate / x;
     a->M = a->in_rate / x;
 
+    a->L <= 0 ? 1 : a->L;
+    a->M <= 0 ? 1 : a->M;
+
     if (a->in_rate < a->out_rate)
         min_rate = a->in_rate;
     else

wdsp/resamplef.cpp

@@ -64,6 +64,9 @@ RESAMPLEF* RESAMPLEF::create_resampleF ( int run, int size, float* in, float* ou
     a->L = out_rate / x;
     a->M = in_rate / x;
 
+    a->L <= 0 ? 1 : a->L;
+    a->M <= 0 ? 1 : a->M;
+
     if (in_rate < out_rate)
         min_rate = in_rate;
     else