Low frequencies Vector Network Analyzer, arduino based.

[Kleinstein]

It dependents on the way the sin / cos values are calculated. With an efficient method (e.g. table or cordic like (Görtzel algorithm or similar)) the DFT for 2 bins is usually more computational efficient than the full FFT and ignoring most of the data. The result is essentially the same, except for possible rounding errors, that should not be significant for 12 Bit raw data.  It is just a different way of calculation. So the effective BW is exactly the same.

The expected advantage in computation time is about about a factor of 10  (details depend on the µC). The DFT code can definitely do real time operation, that is doing the calculation as fast as the data come in at 1 MSPS. Depending on the µC and compiler used one might get something like 8-30 µC cycles per ADC sample.  There is a slight chance the FFT code could also operate real time with 1 MSPS data, though it would be hard and need a faster µC.

The memory use is much lower for the DFT - something like 2 or 4  numbers of about 32 Bit, maybe even 16 bit. For the same window length the result is essentially the same, so no difference in SNR ratio. The DFT solution is flexible to allow longer windows if needed, as there is essentially no memory limitation. In this case a slightly better SNR would be possible. However even with FFT data there is also a chance to use coherent averaging over longer time and thus the same SNR.

[MasterT]

With an efficient method (e.g. table or cordic like (Görtzel algorithm or similar)) the DFT for 2 bins is usually more computational efficient than the full FFT and ignoring most of the data. The result is essentially the same, except for possible rounding errors, that should not be significant for 12 Bit raw data.  It is just a different way of calculation. So the effective BW is exactly the same.
..///.

The memory use is much lower for the DFT - something like 2 or 4  numbers of about 32 Bit, maybe even 16 bit. For the same window length the result is essentially the same, so no difference in SNR ratio. The DFT solution is flexible to allow longer windows if needed, as there is essentially no memory limitation. In this case a slightly better SNR would be possible. However even with FFT data there is also a chance to use coherent averaging over longer time and thus the same SNR.

This is where you do a big mistake. BW is Not the same. SNR is not the same. Stop fooling yourself, there is no free cheese in the mouse trap.
 Reducing fft code down,  in any manner, gaining in memory or clock cycles is become NON relevant, since you are getting the product with lower quality. Lower freq. resolution per bin , lower SNR, missing imaginary part (DCT) etc.
It's same as small print on your insurance policy, bunch of crooks advertise you different variety of the fft, but any variation is
a crap, if you look more carefully into the output.
 Getting back  to the numbers, Görtzel algorithm is unable to select channel 244 Hz wide at 249 756 Hz. Never ever.

[Kleinstein]

The DFT result is the same (except minor rounding differences) as FFT. One needs the second bin to get the imaginary part. So using the 2 bins DFT instead of the FFT is just a different computational method, but essentially the same result.

The FFT algorithm is the most efficient known way if all bins are needed. However it is not the most efficient way anymore if one only needs 2 bins.  

The exact speed does not matter - the question is only if the FFT code could be fast enough to do the calculations in real time, I have some doubt it can with the STM32F3... Even if fast enough, there is the limitation with memory usage.

The possible frequency resolution with the Görtzel Algorithm only depends on the time used and the numerical accuracy. So with enough resolution (e.g. 32 Bit data) the resolution can be way better than 1 Hz in 1 MHz if needed. It just needs time (e.g. a few seconds).

[ogden]

The DFT result is the same (except minor rounding differences) as FFT.

This sentence is said like five times here in this thread? Your patience is unreal

This is where you do a big mistake. BW is Not the same. SNR is not the same. Stop fooling yourself, there is no free cheese in the mouse trap.

Right. BW of FFT and DFT. If such actually exists 

[MasterT]

The possible frequency resolution with the Görtzel Algorithm only depends on the time used and the numerical accuracy. So with enough resolution (e.g. 32 Bit data) the resolution can be way better than 1 Hz in 1 MHz if needed. It just needs time (e.g. a few seconds).

 Time frame is already defined, I have 2048 samples accumulated in data pull, please give me a favor - run Goertzel and see what frequency resolution you would get. 

If you are saying, that 1 000 000 samples has to be processed, to get EQUAL to fft-2048 frequency resolution (244 Hz) , than  calculate clock cycles you spend and compare to fft-2048 , I'm sure that would surprise you. Goertzel gonna to loose quite a lot to get same numbers.

Ratio :  " BW / Clock_Cycles", if you are pulling one (BW) or another "Clock_Cycles" out of equation,  you are trying to cheat.  It's not a  poker.
Both parameters has to be evaluated in the same time frame, numbers of samples, same test condition.  Otherwise you are manipulating data.

[radioactive]

This might help with evaluating.

I started with a C implementation of the Goertzel algorithm from here:  https://www.embedded.com/design/configurable-systems/4024443/The-Goertzel-Algorithm

I hacked it up to add the FFT to compare against.  It depends on libfftw   http://www.fftw.org/

compile it
gcc goertzel_test.c -o goertzel_test -lm -lc -lfftw3

The magnitudes are almost identical over a narrow sweep around the "frequency bin of interest" just to make it visually interesting. 

This is a very good reference from Richard Lyons book
https://www.embedded.com/design/real-world-applications/4401754/Single-tone-detection-with-the-Goertzel-algorithm

[Dave]

Nothing of the kind. Is it engineering forum? I don't like to be involved in any philosophical / political / or abstract mathematics debates. And I 'm not worry if there is a life on a Venus. The point is, engineers should talk in different terms, leaving abstract "better", "less memory"   " faster".  Tell me in dB, uVolts, microseconds, or in this particular case, tell me the Bandwidth in Hz for single bin DFT.

Sure, let's pretend that "better", "less memory" and "faster" are open to any kind of interpretation.

Attached to this message is a MATLAB script that I've written to compare DFT and FFT. Save the file and change the file extension to .m (the forum does not allow .m files). Run it and carefully read through the code.

[MasterT]

Attached to this message is a MATLAB script that I've written to compare DFT and FFT. Save the file and change the file extension to .m (the forum does not allow .m files). Run it and carefully read through the code.

  I'm "self-educated" in C and Java, have EE degree from Russia. But we didn't study matlab in 90-s, sorry for disappointing, no I'm not a digital robot.

This might help with evaluating.

I started with a C implementation of the Goertzel algorithm from here:///

Thanks for sharing, I will test your code in Qt, though I already know the answer.

[Dave]

This is the output of my script.
Top graph is magnitude, bottom is phase. Red line with points is DFT, green dashed line is FFT.

[radioactive]

This might help with evaluating.
I started with a C implementation of the Goertzel algorithm from here:///

Thanks for sharing, I will test your code in Qt, though I already know the answer.

Don't bother.  I just realized you need the complex input/output for your application.  Goertzel algorithm will work almost exactly as FFT/DFT for calculating magnitude from real data, but you need complex input/output for your application.  I will add that I do agree with the others here suggesting that a single-bin DFT would be more optimal for N=2048 and completely identical in output compared to FFT  (including your bandwidth concerns).  They are the same thing.  Others have already said it, but FFT is better only if you need to calculate *all* of the bins.

Sorry for not taking a look at your project video before my earlier posts.  It looks really great!

[Bud]

What i think is fft or not is the least of the worries so to speak. I think some the device architecture needs a few clarifications. What calibrations do you perform to remove systematic errors?

[Mechatrommer]

...me the Bandwidth in Hz for single bin DFT. There is no issue with pure math, consider ratio: "Bandwidth / Clock_cycles_uCPU"....
I'm not argue, that single bin (6-8 math operation code) is faster than 2048 bins, it's all about bandwidth (second page in a row).
For single bin Bandwidth is 250 kHz. 

everyday i learn new thing.. today is Goertzel Algorithm if you just dont care to read here is the snip from Goertzel...

For covering a full spectrum, the Goertzel algorithm has a higher order of complexity than fast Fourier transform (FFT) algorithms, but for computing a small number of selected frequency components, it is more numerically efficient

when Goertzel said computation efficient for single or few bins, it means the same bin as we speak when we calculate DFT or FFT, the resolution for each bin will be = sample rate / number of points / 2. whether DFT, FFT, or Goertzel algorithm they are the same. Goertzel was not fooling around when formulating the single bin algorithm. nobody want to argue with Goertzel that his bin's resolution will be = (sample rate / 2) aka BW. i think the confusion happened (vocabulary) is between DFT and Goertzel. the competition of FFT is not between DFT (FFT is derivative or DFT we know that), but the competition should be between FFT and Goertzel (single bin DFT) and again, the resolution of the single bin of Goertzel will be just the same as bins for DFT or FFT. fwiw...

[radioactive]

...and again, the resolution of the single bin of Goertzel will be just the same as bins for DFT or FFT. fwiw...

True,  I showed that in an example a couple of posts ago, but the issue in this instance is that the Goertzel algorithm is good for real data input.  I think the OP needs complex input for his application (feel free to correct me if I'm wrong here).  The single-bin DFT, as others have pointed out many times, should be more efficient and identical to FFT in this case.  Goertzel is good for applications where you want to optimally detect a single frequency in magnitude  (e.g. DTMF detection).  The computational  savings would be very good for N=2048...

[Mechatrommer]

True,  I showed that in an example a couple of posts ago, but the issue in this instance is that the Goertzel algorithm is good for real data input.  I think the OP needs complex input for his application (feel free to correct me if I'm wrong here).  The single-bin DFT, as others have pointed out many times, should be more efficient and identical to FFT in this case.  Goertzel is good for applications where you want to optimally detect a single frequency in magnitude  (i.e. DTMF detection).  The computational  savings would be very good for N=2048...

there is no mentioning that Goertzel is limited to real part only, (istr there is even faster algortihm than FFT when calculating real part only) so Goertzel  should be capable of calculating both the magnitude and phase, as others linked the Impedance Measurement Techique by Keysight (i think Goertzel algorithm is used in there but not sure didnt read) its just a matter of the OP reading carefully and understanding how Goertzel formulation works, and start coding it in complex number, he got ee degree aye? i dont. he should be able to implement it in C/C++

[radioactive]

there is no mentioning that Goertzel is limited to real part only, (istr there is even faster algortihm than FFT when calculating real part only) so Goertzel  should be capable of calculating both the magnitude and phase,

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.

[edit]   fixed quote block

[Mechatrommer]

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.

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]

as others linked the Impedance Measurement Techique by Keysight (i think Goertzel algorithm is used in there but not sure didnt read)

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]

Don't bother.  I just realized you need the complex input/output for your application. 

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: [Select]

#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;
}

Raw data:

Code: [Select]

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


[ogden]

  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

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?

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

This comes into mind:

[radioactive]

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

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: [Select]

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

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.

[MasterT]

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

//

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.

Right, floating point type is generating rounding noise. 
1. The problem with a Goertzel (as any IIR) it demands 64-bits math, same time fft runs smooth and easy on 32-bits. SAM3X8E is 32-bits, w/o FPU, so difference in the time execution has to be verified on my DUE. Difference in 64-bits Goertzel vs 32-bits FFT.
2. Next, in the initial code you omitted windowing function, and this comparison 3x in timing is not correct, if

fft-time / grt-time = 3, than

(fft-time + windowing) / (grt-time  + windowing)  != 3.

3. Phase data I calculated out of Goertzel is not correct, error seems consistent 22.5 degree "off" - that may be corrected or not, anyway needs more testing to be verified in all circumstances.

//
All this hussle was started by people persistently pointing me that fft for single bin doesn't make any sense.   

[Bud]

This exersize in comparing sampling algorithm is pointless without a calibration procedure. You can fft until cows come home and get meaningless results because of the systematic errors not being removed.

[ogden]

All this hussle was started by people persistently pointing me that fft for single bin doesn't make any sense.

Hell, no! It was started by you suggesting that nobody but you knows The answer to the ultimate question of life, the universe and DSP:

Manny people don't understand the absolute supremacy of the FFT on any other method of data processing. The question of "why using FFT instead of DCT?" is the same as "why to sample data by ADC and than process in DSP core, if we can drive a couple of MUXes by quadrature shifted clock?"  Because FFT provides processing gain, increasing SNR ratio by >30 dB(!).

All this  delusion about DFT /DCT /Wavelets  arises from misunderstanding basic fact, that FFT is the most efficient, accurate and fastest algorithm.

There are things that not everyone would ever understand, despite "correct" vocabulary they use.

[radioactive]

2. Next, in the initial code you omitted windowing function, and this comparison 3x in timing is not correct, if
...
All this hussle was started by people persistently pointing me that fft for single bin doesn't make any sense.

I didn't omit the windowing function, you did. 

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

I agree that 3x timing difference is not correct.  This was just the timing difference for calling full-sweep with your SNR functions for grt vs fft  1000 times..  The Goertzel function would be thousands of times faster for your original specs:  N=2048,  single-bin.

I think the real point of all this was showing that the FFT, DFT, Goertzel are all going to give you the same results for a single bin.  In my opinion, the real take-away should be this:   

For optimal efficiency (sometimes this doesn't matter),

1) use FFT if you need results across the entire channel.   
2)  use a single-bin DFT if you need a single bin with complex input, 
3) use Goertzel if you need a single bin with real input only. 

[edit]

Btw,  thank you for the FFT code.  It is nice looking code.

[MasterT]

I didn't omit the windowing function, you did. 

I'm referring to your first post, where you comparing fftw vs goertzel:
https://www.eevblog.com/forum/projects/low-frequencies-vector-network-analyzer-arduino-based/msg1844921/#msg1844921

Why you didn't apply windowing, as you should?  Than, would be a huge difference in the two charts, would not it?

And you are ignoring the fact, that Goertzel is meaningless w/o 64-bits architecture.