Electronics > Projects, Designs, and Technical Stuff
Low frequencies Vector Network Analyzer, arduino based.
<< < (14/20) > >>
Mechatrommer:

--- Quote from: radioactive on September 24, 2018, 02:39:36 am ---I think at the point that you implement that, it probably won't have any gains over the single-bin DFT  (again correct me if I'm wrong).  I used Goertzel in an FPGA application where I didn't have enough cells for FFT in the past, but I was only interested in magnitude.  It worked quite well.

--- End quote ---
not sure if Goertzel is optimized single bin DFT for real part, cant differ, but i think even when calculating imaginary part in single-bin DFT (or maybe Goertzel?), the complexity should still near O(n) or O(2n)? still better than O(n.log(n)) if n > 100 or something.

as from the hardware side, there are interesting topology or some sort of impedance divider going on, but its not clear since there are few broken traces so i dont know where the DUT will be connected to. i agree i have to vote this device is more toward "Impendace Analyzer" or "Bode Plotter" as video in the OP, there is no direct S21 measurement let alone S11. so VNA is probably just a marketing hype from Analog Devices or the OP. i look up http://www.analog.com/media/en/technical-documentation/data-sheets/AD5933.pdf the OP circuit front end doesnt resemble anything like the diagrams in the datasheet nor this..


from http://www.analog.com/en/applications/markets/instrumentation-and-measurement-pavilion-home/electronic-test-and-measurement/rf-signal-analyzers-and-vector-network-analyzers.html
ogden:

--- Quote from: Mechatrommer on September 24, 2018, 02:32:32 am ---as others linked the Impedance Measurement Techique by Keysight (i think Goertzel algorithm is used in there but not sure didnt read)

--- End quote ---

They use DFT. Also it shall be noted that VNA's/LCR's usually does not do direct sampling but have heterodyne architecture with fixed, low IF frequency so it can be sampled using high resolution ADC's.

Agilent E5063A Network Analyzer Help:

The digital filter performs a discrete Fourier transformation (DFT) and picks
up IF signals. Each IF signal is then converted into a complex number that
has a real part and an imaginary part. The IF bandwidth of the analyzer is
equivalent to the bandwidth of the DFT filter.

Google for "e5063a_operation.pdf"
MasterT:

--- Quote from: radioactive on September 24, 2018, 12:28:27 am ---Don't bother.  I just realized you need the complex input/output for your application. 

--- End quote ---
No, pure real data in this project, I'm not doing reverse IFFT, and adc is always real  :)
Thanks again for posting a code, it has more value for me than a million words. I really appreciate, when conversation is going more specific, down to earth - digits in this particular case.

Running this piece of code up and down, I finally created the test procedure, that Magnifying the small prints in your insurance policy, I was talking about.   

Excuse for delay to everyone, who is reading this, was busy doing testing with a code. Here is what I came up to.

Test condition:
  //https://www.embedded.com/design/configurable-systems/4024443/The-Goertzel-Algorithm

  Goertzel algorithm is initialized to detect Target Frequency = 500.

  //FFT is my own, feel free to replace it with whatever you like, FFTW, CMSIS-DSP pack etc.

  FFT algorithm doesn't need to be initialized, magnitude is extracted at:
  Resolution = sampling freq. /256 = 31.25 Hz
  Bin = target /resolution = 500 / 31.25 = 16.

  No windowing function is applied in this test, not to obscure the purpose  of the testing.   

Input:
  Single tone sine wave, no phase data, no noise added in this test.
  Test tone is swapped from 0 to Nyquist = sampling /2.
  Swapped tone is stepping at 31.25 Hz, since there was no windowing function, stepping goes at exact "bin's frequencies" in
  order to prevent any side lobes that would be generated otherwise and ruin the objective of this test.
 
  When test tone "hits" a target frequency, detection happened and magnitude is recorded.
  Any other frequencies, except target, becomes in this case "interference", "noise" or "adjacent channel power ratio".

  Obviously, dsp core  (Goertzel or FFT)  Must eliminate any interference getting into the target frequency (detection band).
  The magnitude of the interference is preset to -40 dB, 20.47 riding over 2047 DC level -  12-bits ADC approximation.

  To precisely evaluate this "interference withstanding" or "noise reduction" capability  of the dsp-core SINAD is calculated :
   
   sinad = (20.0 * log10(temp /tot_magnt));

  where tot_magnt is RMS of all interference getting into detection band,
  and variable "temp" is assigned to magnitude at Target frequency.

Results:
 bin  2620.176270 tot_m     1.124978 sinad_grt    67.343732
 bin  2620.160156 tot_m     0.020573 sinad_fft   102.100667

You see, there is 35 dB (!) difference in the SNR.
  Goertzel is miserable failed on the objective.

NB1:
Not sure what is correct name for this test would be. I already used hardware analogy with pre-selector,  so "selectivity" comes to mind.

NB2:
I agree, that using "bandwidth" was wrong on my side, indeed BW  (defined in terms of sqrt(2))  is practically identical  between
single bin (Goertzel) and FFT.

  Can't explain any better.

Code:

--- Code: ---#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define SAMPLING_RATE       8000.0 //  8kHz
#define TARGET_FREQUENCY    500     // 1875//941.0 //941 Hz
#define N                   256     // 205 //Block size


//bin  2559.962402      tot_m 1.118038 sinad_grt    67.195540
//bin  2559.995117      tot_m 0.020119 sinad_fft   102.092452
        float coeff;
        float Q1;
        float Q2;
        float sine;
        float cosine;
        double testData[N];

/*
//bin 2559.989502       tot_m 0.182476      sinad_grt   82.940644
//bin 2559.995117       tot_m 0.020119      sinad_fft  102.092452
        double coeff;
        double Q1;
        double Q2;
        double sine;
        double cosine;
        double testData[N];
*/

#define FFT_SIZE      N
#define MIRROR        (FFT_SIZE / 2)

void ResetGoertzel(void)
{
  Q2 = 0;
  Q1 = 0;
}

void InitGoertzel(void)
{
  int k;
  float floatN;
  float omega;

  floatN = (float) N;
  k = (int) (0.5 + ((floatN * TARGET_FREQUENCY) / SAMPLING_RATE));
  omega = (2.0 * M_PI * k) / floatN;
  sine = sin(omega);
  cosine = cos(omega);
  coeff = 2.0 * cosine;

  printf("For SAMPLING_RATE = %f", SAMPLING_RATE);
  printf(" N = %d", N);
  printf(" and FREQUENCY = %d,\n", TARGET_FREQUENCY);
  printf("k = %d and coeff = %f, sine = %f, cosine = %f\n\n", k, coeff, sine, cosine);

  ResetGoertzel();
}

void Processfloat(float sample)
{
  float Q0;
  Q0 = coeff * Q1 - Q2 + (float) sample;
  Q2 = Q1;
  Q1 = Q0;
}

void GetRealImag(float *realPart, float *imagPart)
{
  *realPart = (Q1 - Q2 * cosine);
  *imagPart = (Q2 * sine);
}

float GetMagnitudeSquared(void)
{
  float result;

  result = Q1 * Q1 + Q2 * Q2 - Q1 * Q2 * coeff;
  return result;
}

void Generate(float frequency)
{
  int index;
  float step;

  step = frequency * ((2.0 * M_PI) / SAMPLING_RATE);
  for(  index = 0; index < N; index++) {
    testData[index] = 20.470 * sin(index * step) +2047.0;
    }
}

void GenerateAndTest_GRT(float frequency, float *result)
{
  int index;

  float magnitudeSquared;
  float magnitude;
  float real;
  float imag;

  printf("F1=%7.2f, ", frequency);
  Generate(frequency);
/*
  // 2.Windowing
  double base_angle = 2.0 *M_PI / (FFT_SIZE -1);
  for( int16_t i = 0; i < FFT_SIZE; i++ ) {
    double angle = i *base_angle;
    float temp = 0.35875 -(0.48829*cos(angle)) +(0.14128*cos(2*angle)) -(0.01168*cos(3*angle));
//    testData[i] *= temp;
      testData[i] *= (0.5 - (0.5 * cos( 2 *M_PI *i /(FFT_SIZE -1))));
      }
*/

  for(  index = 0; index < N; index++) {
    //Processfloat(testData[index]);
    float Q0;
    Q0 = coeff * Q1 - Q2 + (float)testData[index];
    Q2 = Q1;
    Q1 = Q0;
    }
  GetRealImag(&real, &imag);
    //  printf("r= %15.5f i= %15.5f", real, imag);
  magnitudeSquared = real*real + imag*imag;
    //  printf("rel mag^2=%16.5f   ", magnitudeSquared);
  magnitude = sqrt(magnitudeSquared);
    //  printf("mag=%15.5f\n", magnitude);
  printf("mag=%15.5f", magnitude);
  ResetGoertzel();

  *result = magnitude;
}

// FFT
void rev_bin( float *fr, int16_t fft_n)
{
    int m, mr, nn, l;
    float tr;

    mr = 0;
    nn = fft_n - 1;

    for (m=1; m<=nn; ++m) {
            l = fft_n;
         do {
             l >>= 1;
            } while (mr+l > nn);

            mr = (mr & (l-1)) + l;

        if (mr <= m)
             continue;
            tr = fr[m];
            fr[m] = fr[mr];
            fr[mr] = tr;
    }
}

void get_MagnitF(float *fr)
{
  for(int  i = 1; i < MIRROR; i++) {
    float real = fr[i];
    float imag = fr[FFT_SIZE -i];
    float Magnit = (float) (sqrt((real * real) + (imag * imag)));
    fr[i] = Magnit;
  }
}

void rfft33(float X[], int Nd)
{
    /****************************************************************************
    *   rfft(float X[],int N)                                                   *
    *     A real-valued, in-place, split-radix FFT program                      *
    *     Decimation-in-time, cos/sin in second loop                            *
    *     Length is N=2**M (i.e. N must be power of 2--no error checking)       *
    *                                                                           *
    * Original Fortran code by Sorensen; published in H.V. Sorensen, D.L. Jones,*
    * M.T. Heideman, C.S. Burrus (1987) Real-valued fast fourier transform      *
    * algorithms.  IEEE Trans on Acoustics, Speech, & Signal Processing, 35,    *
    * 849-863.  Adapted to C by Bill Simpson, 1995  wsimpson@uwinnipeg.ca       *
    * --------------------------------------------------------------------------*
    * C/C++ language bugs fixing & correction:                                  *
    * 1. C-style array start adress [0];                                        *
    * 2. SIN() and COS tweed factors change place;                              *
    * 3. Rev Bin - replaced by new fastest version.                             *
    * made by Anatoly Kuzmenko, 2017,  anatolyk69@gmail.com                     *
    ****************************************************************************/

  int I, I0, I1, I2, I3, I4, I5, I6, I7, I8, IS, ID;
  int J, K, M, N2, N4, N8;
  float A, A3, CC1, SS1, CC3, SS3, E, R1;
  float T1, T2, T3, T4, T5, T6;

  M = (int)(log(N) / log(2.0));             /* N=2^M */

  /* ----Length two butterflies--------------------------------------------- */
  IS = 0;
  ID = 4;
  do {
    for (I0 = IS; I0 < Nd; I0 += ID) {
      I1    = I0 + 1;
      R1    = X[I0];
      X[I0] = R1 + X[I1];
      X[I1] = R1 - X[I1];
    }
    IS = 2 * ID - 2;
    ID = 4 * ID;
  } while (IS < Nd);

  /* ----L shaped butterflies----------------------------------------------- */
  N2 = 2;
  for (K = 1; K < M; K++) {
    N2    = N2 * 2;
    N4    = N2 / 4;
    N8    = N2 / 8;
    E     = (float) 6.2831853071719586f / N2;
    IS    = 0;
    ID    = N2 * 2;
    do {
      for ( I = IS; I < Nd; I += ID) {
        I1 = I; // + 1;
        I2 = I1 + N4;
        I3 = I2 + N4;
        I4 = I3 + N4;
        T1 = X[I4] + X[I3];
        X[I4] = X[I4] - X[I3];
        X[I3] = X[I1] - T1;
        X[I1] = X[I1] + T1;
        if (N4 != 1) {
          I1 += N8;
          I2 += N8;
          I3 += N8;
          I4 += N8;
          T1 = (X[I3] + X[I4]) * .7071067811865475244f;
          T2 = (X[I3] - X[I4]) * .7071067811865475244f;
          X[I4] = X[I2] - T1;
          X[I3] = -X[I2] - T1;
          X[I2] = X[I1] - T2;
          X[I1] = X[I1] + T2;
        }
      }
      IS = 2 * ID - N2;
      ID = 4 * ID;
    } while ( IS < Nd );
    A = E;
    for ( J = 1; J < N8; J++) {
      A = (float)J * E;
      A3 = 3.0 * A;
      CC1   = cos( A);
      SS1   = sin( A);
      CC3   = cos(A3);
      SS3   = sin(A3);
      IS = 0;
      ID = 2 * N2;
      do {
        for ( I = IS; I < Nd; I += ID) {
          I1 = I + J;
          I2 = I1 + N4;
          I3 = I2 + N4;
          I4 = I3 + N4;
          I5 = I + N4 - J;// + 2;
          I6 = I5 + N4;
          I7 = I6 + N4;
          I8 = I7 + N4;
          T1 = X[I3] * CC1 + X[I7] * SS1;
          T2 = X[I7] * CC1 - X[I3] * SS1;
          T3 = X[I4] * CC3 + X[I8] * SS3;
          T4 = X[I8] * CC3 - X[I4] * SS3;
          T5 = T1 + T3;
          T6 = T2 + T4;
          T3 = T1 - T3;
          T4 = T2 - T4;
          T2 = X[I6] + T6;
          X[I3] = T6 - X[I6];
          X[I8] = T2;
          T2    = X[I2] - T3;
          X[I7] = -X[I2] - T3;
          X[I4] = T2;
          T1    = X[I1] + T5;
          X[I6] = X[I1] - T5;
          X[I1] = T1;
          T1    = X[I5] + T4;
          X[I5] = X[I5] - T4;
          X[I2] = T1;
        }
        IS = 2 * ID - N2;
        ID = 4 * ID;
      } while ( IS < Nd );
    }
  }
}

void GenerateAndTest_FFT(float frequency, float *result)
{
  float magnitude;

  Generate(frequency);

  float fr[FFT_SIZE];

  for( int16_t i = 0; i < FFT_SIZE; i++ ) {
    fr[i] = (float)testData[i];
    }
/*
  printf(" X[i] Input Data \n");
  for ( int  i = 0; i < FFT_SIZE; i++)
  {  printf("%8.1f", fr[i]);
     if ((i+1)%16 == 0) printf("\n");
   }
  printf("\n.........................................\n");
*/

/*
  // 2.Windowing
       double base_angle = 2.0 *M_PI / (FFT_SIZE -1);
       for( int16_t i = 0; i < FFT_SIZE; i++ ) {
          double angle = i *base_angle;
          float temp = 0.35875 -(0.48829*cos(angle)) +(0.14128*cos(2*angle)) -(0.01168*cos(3*angle));
//        fr[i] *= temp;
         fr[i] *= (0.5 - (0.5 * cos( 2 *M_PI *i /(FFT_SIZE -1))));
       }
*/
/*
       printf(" X[i] Input Data \n");
       for ( int  i = 0; i < FFT_SIZE; i++)
       {  printf("%8d", fr[i]);
          if ((i+1)%16 == 0) printf("\n");
        }
       printf("\n.........................................\n");
*/

       rev_bin( fr, (int16_t) FFT_SIZE);
/*
       printf(" X[i] Input Data \n");
       for ( int  i = 0; i < FFT_SIZE; i++)
       {  printf("%8.0f", fr[i]);
          if ((i+1)%16 == 0) printf("\n");
        }
       printf("\n.........................................\n");
*/

       rfft33( fr, FFT_SIZE);
/*
       printf(" X[i] Output Data \n");
       for ( int  i = 0; i < FFT_SIZE; i++)
       {  printf("%8.0f", fr[i]);
          if ((i+1)%16 == 0) printf("\n");
        }
       printf("\n.........................................\n");
*/

       get_MagnitF(fr);
/*
       printf(" X[i] Output Data \n");
       for ( int  i = 0; i < MIRROR; i++)
       {  printf("%8.0f", fr[i]);
          if ((i+1)%16 == 0) printf("\n");
        }
       printf("\n.........................................\n");
*/
       magnitude = fr[16];
  printf("  mag=%15.5f\n", magnitude);

  *result = magnitude;
}


int main(void)
{
  float freq;

  InitGoertzel();
  float test;
  float snr_grt[MIRROR] = { 0.0};
  float snr_fft[MIRROR] = { 0.0};

  int   index = 0;
  float step = SAMPLING_RATE / FFT_SIZE;

    //  for (freq = TARGET_FREQUENCY - 300; freq <= TARGET_FREQUENCY + 300; freq += 15)
    //  for( freq = 0.0; freq < SAMPLING_RATE /2.0; freq += 31.25) {
  for( freq = 0.0; freq < SAMPLING_RATE /2.0; freq += step) {

    GenerateAndTest_GRT(freq, &test);
    snr_grt[index] = test;

    GenerateAndTest_FFT(freq, &test);
    snr_fft[index] = test;

    if(index >= MIRROR) break;
    index++;
    }

/*
  printf(" snr_grt \n");
  for( int  i = 0; i < MIRROR; i++) {
    printf("%15.6f", snr_grt[i]);
    if ((i+1)%8 == 0) printf("\n");
    }
  printf("\n.........................................\n");
  printf(" snr_fft \n");
  for( int  i = 0; i < MIRROR; i++) {
    printf("%15.6f", snr_fft[i]);
    if ((i+1)%8 == 0) printf("\n");
    }
  printf("\n.........................................\n");
*/

    int    bin_numbr    = 16;
    double tot_magnt    = 0;
    double temp         = 0.0;
    double sinad        = 0.0;

    for( int i = 0; i < MIRROR; i++) {
      if(i == bin_numbr) continue;
      temp = snr_grt[i];
      tot_magnt += (temp * temp);
      }

    tot_magnt = sqrt(tot_magnt);
    if (tot_magnt == 0) tot_magnt = 0.000001;
    temp = snr_grt[bin_numbr];

    sinad = (20.0 * log10(temp /tot_magnt));
    printf("\n bin %12.6f", temp);
    printf(" tot_m %12.6f", tot_magnt);
    printf(" sinad_grt %12.6f", sinad);

    tot_magnt    = 0;
    temp         = 0.0;
    sinad        = 0.0;

    for( int i = 0; i < MIRROR; i++) {
      if(i == bin_numbr) continue;
      temp = snr_fft[i];
      tot_magnt += (temp * temp);
      }

    tot_magnt = sqrt(tot_magnt);
    if (tot_magnt == 0) tot_magnt = 0.000001;
    temp = snr_fft[bin_numbr];

    sinad = (20.0 * log10(temp /tot_magnt));
    printf("\n bin %12.6f", temp);
    printf(" tot_m %12.6f", tot_magnt);
    printf(" sinad_fft %12.6f", sinad);


  return 0;
}
--- End code ---


Raw data:


--- Code: ---For SAMPLING_RATE = 8000.000000 N = 256 and FREQUENCY = 500,
k = 16 and coeff = 1.847759, sine = 0.382683, cosine = 0.923880

F1=   0.00, mag=        0.10611  mag=        0.00000
F1=  31.25, mag=        0.09423  mag=        0.00000
F1=  62.50, mag=        0.10754  mag=        0.00000
F1=  93.75, mag=        0.08793  mag=        0.00000
F1= 125.00, mag=        0.11068  mag=        0.00000
F1= 156.25, mag=        0.07982  mag=        0.00000
F1= 187.50, mag=        0.11203  mag=        0.00000
F1= 218.75, mag=        0.07678  mag=        0.00000
F1= 250.00, mag=        0.09029  mag=        0.00000
F1= 281.25, mag=        0.10285  mag=        0.00000
F1= 312.50, mag=        0.08015  mag=        0.00000
F1= 343.75, mag=        0.10285  mag=        0.00000
F1= 375.00, mag=        0.09063  mag=        0.00000
F1= 406.25, mag=        0.09256  mag=        0.00000
F1= 437.50, mag=        0.08764  mag=        0.00000
F1= 468.75, mag=        0.09484  mag=        0.00000
F1= 500.00, mag=     2620.17627  mag=     2620.16016
F1= 531.25, mag=        0.10418  mag=        0.00000
F1= 562.50, mag=        0.10534  mag=        0.00000
F1= 593.75, mag=        0.10568  mag=        0.00000
F1= 625.00, mag=        0.09694  mag=        0.00000
F1= 656.25, mag=        0.08705  mag=        0.00000
F1= 687.50, mag=        0.09377  mag=        0.00000
F1= 718.75, mag=        0.09259  mag=        0.00000
F1= 750.00, mag=        0.10259  mag=        0.00195
F1= 781.25, mag=        0.08477  mag=        0.00000
F1= 812.50, mag=        0.10585  mag=        0.00000
F1= 843.75, mag=        0.10384  mag=        0.00000
F1= 875.00, mag=        0.10069  mag=        0.00000
F1= 906.25, mag=        0.10469  mag=        0.00000
F1= 937.50, mag=        0.11233  mag=        0.00000
F1= 968.75, mag=        0.09861  mag=        0.00000
F1=1000.00, mag=        0.11733  mag=        0.00391
F1=1031.25, mag=        0.10588  mag=        0.00000
F1=1062.50, mag=        0.09900  mag=        0.00000
F1=1093.75, mag=        0.09996  mag=        0.00000
F1=1125.00, mag=        0.09220  mag=        0.00000
F1=1156.25, mag=        0.11435  mag=        0.00000
F1=1187.50, mag=        0.09662  mag=        0.00000
F1=1218.75, mag=        0.10372  mag=        0.00000
F1=1250.00, mag=        0.09398  mag=        0.00276
F1=1281.25, mag=        0.09437  mag=        0.00000
F1=1312.50, mag=        0.09865  mag=        0.00000
F1=1343.75, mag=        0.11403  mag=        0.00000
F1=1375.00, mag=        0.10098  mag=        0.00000
F1=1406.25, mag=        0.08598  mag=        0.00000
F1=1437.50, mag=        0.10105  mag=        0.00000
F1=1468.75, mag=        0.09867  mag=        0.00000
F1=1500.00, mag=        0.10565  mag=        0.00485
F1=1531.25, mag=        0.09424  mag=        0.00000
F1=1562.50, mag=        0.09212  mag=        0.00000
F1=1593.75, mag=        0.10664  mag=        0.00000
F1=1625.00, mag=        0.11636  mag=        0.00195
F1=1656.25, mag=        0.08445  mag=        0.00000
F1=1687.50, mag=        0.10084  mag=        0.00000
F1=1718.75, mag=        0.09308  mag=        0.00000
F1=1750.00, mag=        0.10372  mag=        0.00249
F1=1781.25, mag=        0.09571  mag=        0.00000
F1=1812.50, mag=        0.09840  mag=        0.00000
F1=1843.75, mag=        0.10246  mag=        0.00000
F1=1875.00, mag=        0.11236  mag=        0.00195
F1=1906.25, mag=        0.11003  mag=        0.00000
F1=1937.50, mag=        0.10503  mag=        0.00000
F1=1968.75, mag=        0.11180  mag=        0.00000
F1=2000.00, mag=        0.09783  mag=        0.00254
F1=2031.25, mag=        0.09688  mag=        0.00000
F1=2062.50, mag=        0.09493  mag=        0.00000
F1=2093.75, mag=        0.09357  mag=        0.00000
F1=2125.00, mag=        0.09151  mag=        0.00000
F1=2156.25, mag=        0.10432  mag=        0.00000
F1=2187.50, mag=        0.09754  mag=        0.00000
F1=2218.75, mag=        0.10882  mag=        0.00195
F1=2250.00, mag=        0.10045  mag=        0.00202
F1=2281.25, mag=        0.11314  mag=        0.00000
F1=2312.50, mag=        0.09258  mag=        0.00000
F1=2343.75, mag=        0.09818  mag=        0.00000
F1=2375.00, mag=        0.11400  mag=        0.00195
F1=2406.25, mag=        0.09444  mag=        0.00000
F1=2437.50, mag=        0.10033  mag=        0.00000
F1=2468.75, mag=        0.07606  mag=        0.00000
F1=2500.00, mag=        0.10632  mag=        0.00465
F1=2531.25, mag=        0.11568  mag=        0.00000
F1=2562.50, mag=        0.11631  mag=        0.00000
F1=2593.75, mag=        0.10031  mag=        0.00000
F1=2625.00, mag=        0.10132  mag=        0.00000
F1=2656.25, mag=        0.10407  mag=        0.00000
F1=2687.50, mag=        0.10097  mag=        0.00000
F1=2718.75, mag=        0.07938  mag=        0.00195
F1=2750.00, mag=        0.09015  mag=        0.00383
F1=2781.25, mag=        0.11763  mag=        0.00195
F1=2812.50, mag=        0.09597  mag=        0.00000
F1=2843.75, mag=        0.10604  mag=        0.00000
F1=2875.00, mag=        0.11056  mag=        0.00195
F1=2906.25, mag=        0.09866  mag=        0.00000
F1=2937.50, mag=        0.09622  mag=        0.00000
F1=2968.75, mag=        0.09654  mag=        0.00361
F1=3000.00, mag=        0.08598  mag=        0.00149
F1=3031.25, mag=        0.09443  mag=        0.00289
F1=3062.50, mag=        0.10529  mag=        0.00383
F1=3093.75, mag=        0.10426  mag=        0.00195
F1=3125.00, mag=        0.07579  mag=        0.00361
F1=3156.25, mag=        0.08571  mag=        0.00000
F1=3187.50, mag=        0.11280  mag=        0.00000
F1=3218.75, mag=        0.07765  mag=        0.00000
F1=3250.00, mag=        0.09780  mag=        0.00195
F1=3281.25, mag=        0.08297  mag=        0.00000
F1=3312.50, mag=        0.09670  mag=        0.00195
F1=3343.75, mag=        0.11009  mag=        0.00000
F1=3375.00, mag=        0.11433  mag=        0.00345
F1=3406.25, mag=        0.10358  mag=        0.00000
F1=3437.50, mag=        0.09293  mag=        0.00195
F1=3468.75, mag=        0.11226  mag=        0.00868
F1=3500.00, mag=        0.12093  mag=        0.00379
F1=3531.25, mag=        0.10341  mag=        0.00000
F1=3562.50, mag=        0.10189  mag=        0.00076
F1=3593.75, mag=        0.11852  mag=        0.00195
F1=3625.00, mag=        0.08856  mag=        0.00217
F1=3656.25, mag=        0.11255  mag=        0.00000
F1=3687.50, mag=        0.10105  mag=        0.00195
F1=3718.75, mag=        0.11565  mag=        0.00518
F1=3750.00, mag=        0.08955  mag=        0.00293
F1=3781.25, mag=        0.09275  mag=        0.00195
F1=3812.50, mag=        0.10781  mag=        0.00000
F1=3843.75, mag=        0.10727  mag=        0.00000
F1=3875.00, mag=        0.09191  mag=        0.00546
F1=3906.25, mag=        0.09398  mag=        0.00656
F1=3937.50, mag=        0.09004  mag=        0.00000
F1=3968.75, mag=        0.10833  mag=        0.00000

 bin  2620.176270 tot_m     1.124978 sinad_grt    67.343732
 bin  2620.160156 tot_m     0.020573 sinad_fft   102.100667

--- End code ---
ogden:

--- Quote from: MasterT on September 25, 2018, 04:41:54 pm ---  Obviously, dsp core  (Goertzel or FFT)  Must eliminate any interference getting into the target frequency (detection band).
  The magnitude of the interference is preset to -40 dB, 20.47 riding over 2047 DC level -  12-bits ADC approximation.

  To precisely evaluate this "interference withstanding" or "noise reduction" capability  of the dsp-core SINAD is calculated

--- End quote ---

Author (who does not read my posts anymore) "conveniently forgot" that only FFT have granularity and bins. When you use DFT in VNA/LCR_meter - you "tune" it precisely to the test frequency and such thing as "interference withstanding" or "noise reduction" is irrelevant. Of course Goertzel with float coefficients can't win float FFT implementation that does not even use tables but plain sin() cos(). Why would someone need to write such a test to prove that?


--- Quote ---Not sure what is correct name for this test would be. I already used hardware analogy with pre-selector,  so "selectivity" comes to mind.

--- End quote ---

This comes into mind: :bullshit:
radioactive:

--- Quote from: MasterT on September 25, 2018, 04:41:54 pm ---...
  To precisely evaluate this "interference withstanding" or "noise reduction" capability  of the dsp-core SINAD is calculated :
   
   sinad = (20.0 * log10(temp /tot_magnt));

  where tot_magnt is RMS of all interference getting into detection band,
  and variable "temp" is assigned to magnitude at Target frequency.

Results:
 bin  2620.176270 tot_m     1.124978 sinad_grt    67.343732
 bin  2620.160156 tot_m     0.020573 sinad_fft   102.100667

You see, there is 35 dB (!) difference in the SNR.
  Goertzel is miserable failed on the objective.

--- End quote ---

I don't understand the need to calculate SNR over the entire channel if you are only interested in the results of a single bin.  I think the difference in results are because of floating point precision.   Just to test that idea,  I changed every float type in your test code to double and then implemented a test to run on a Cortex M7 MCU (results are the same on PC).  I also added a test to compare the execution time of Goertzel vs your implementation of FFT.  With doubles it appears that Goertzel has SNR improvement of 94 dB over FFT (using your calculations of SNR).  It executes about 3x faster with your current configuration of N=256.   I didn't try to optimize the initGoertzel() function either, so might gain a little there too.


--- Code: ---static int grt_count;
static int fft_count;
static uint32_t grt_st_time;
static uint32_t fft_st_time;
static uint32_t grt_end_time;
static uint32_t fft_end_time;


//in main while loop
//

while(1) {

    current_time = HAL_GetTick(); //milliseconds


    if(grt_count<1000) {
        if(grt_count==0) grt_st_time=current_time;
        test_grt(0);
        grt_count++;
        if(grt_count==1000) {
          grt_end_time=current_time;
          printf("\r\ngrt total time for 1k runs %d", grt_end_time-grt_st_time);
          test_grt(1);  //print the SNR results
        }
    }
    else if(fft_count<1000) {
        if(grt_count==0) fft_st_time=current_time;
        test_fft(0);
        fft_count++;
        if(fft_count==1000) {
          fft_end_time=current_time;
          printf("\r\nfft total time for 1k runs %d", fft_end_time-fft_st_time);
          test_fft(1);  //print the SNR results
        }
    }
}
--- End code ---


Results:

grt total time for 1k runs 21173
 bin  2620.160000 tot_m     0.000000 sinad_grt   251.172854
fft total time for 1k runs 63307
 bin  2620.159989 tot_m     0.000040 sinad_fft   156.307001

Even so, I guess that savings in time might not be worth messing with in your case.  If the FFT is working great and there is no need to improve the timing,  you have already validated the FFT.  I don't see a need to change it.  It gets the job done.
Navigation
Message Index
Next page
Previous page
There was an error while thanking
Thanking...

Go to full version
Powered by SMFPacks Advanced Attachments Uploader Mod