Create a twiddle table with a specific order?

Question

I was hoping to create a twiddle table with a specific ordering for a specific purpose, here is the initial code:

#define TWIDDLE_LIMIT 64
#define PI 3.1415927

float *twiddle_real;
float *twiddle_imag;

void main()
{
  int N = 256;
  int TW_size = TWIDDLE_LIMIT + (TWIDDLE_LIMIT>>2);

  twiddle_real  = malloc(TW_size * sizeof(float));
  twiddle_imag  = malloc(TW_size * sizeof(float));

  int i;
  for(i=0; i<TWIDDLE_LIMIT; i++)
  {
    twiddle_real[i] = (float)   cos((float)i * 2.0 * PI / (float)N);
    twiddle_imag[i] = (float) - sin((float)i * 2.0 * PI / (float)N);
  }
  for(int a=0; a<TWIDDLE_LIMIT; a++)
    printf("RE = %f \t IM = %f \n",twiddle_real[a],twiddle_imag[a]);
}

And i get this kind of result:

RE = 1.000000    IM = -0.000000 //64 lines
RE = 0.999699    IM = -0.024541 
RE = 0.998795    IM = -0.049068 
RE = 0.997290    IM = -0.073565 
RE = 0.995185    IM = -0.098017 
RE = 0.992480    IM = -0.122411 
RE = 0.989177    IM = -0.146730 
RE = 0.985278    IM = -0.170962 
RE = 0.980785    IM = -0.195090 
RE = 0.975702    IM = -0.219101 
RE = 0.970031    IM = -0.242980 
RE = 0.963776    IM = -0.266713 
RE = 0.956940    IM = -0.290285 
RE = 0.949528    IM = -0.313682 
RE = 0.941544    IM = -0.336890 
RE = 0.932993    IM = -0.359895 
RE = 0.923880    IM = -0.382683 
RE = 0.914210    IM = -0.405241 
RE = 0.903989    IM = -0.427555 
RE = 0.893224    IM = -0.449611 
RE = 0.881921    IM = -0.471397 
RE = 0.870087    IM = -0.492898 
RE = 0.857729    IM = -0.514103 
RE = 0.844854    IM = -0.534998 
RE = 0.831470    IM = -0.555570 
RE = 0.817585    IM = -0.575808 
RE = 0.803208    IM = -0.595699 
RE = 0.788346    IM = -0.615232 
RE = 0.773010    IM = -0.634393 
RE = 0.757209    IM = -0.653173 
RE = 0.740951    IM = -0.671559 
RE = 0.724247    IM = -0.689541 
RE = 0.707107    IM = -0.707107 
RE = 0.689541    IM = -0.724247 
RE = 0.671559    IM = -0.740951 
RE = 0.653173    IM = -0.757209 
RE = 0.634393    IM = -0.773010 
RE = 0.615232    IM = -0.788346 
RE = 0.595699    IM = -0.803208 
RE = 0.575808    IM = -0.817585 
RE = 0.555570    IM = -0.831470 
RE = 0.534998    IM = -0.844854 
RE = 0.514103    IM = -0.857729 
RE = 0.492898    IM = -0.870087 
RE = 0.471397    IM = -0.881921 
RE = 0.449611    IM = -0.893224 
RE = 0.427555    IM = -0.903989 
RE = 0.405241    IM = -0.914210 
RE = 0.382683    IM = -0.923880 
RE = 0.359895    IM = -0.932993 
RE = 0.336890    IM = -0.941544 
RE = 0.313682    IM = -0.949528 
RE = 0.290285    IM = -0.956940 
RE = 0.266713    IM = -0.963776 
RE = 0.242980    IM = -0.970031 
RE = 0.219101    IM = -0.975702 
RE = 0.195090    IM = -0.980785 
RE = 0.170962    IM = -0.985278 
RE = 0.146730    IM = -0.989177 
RE = 0.122411    IM = -0.992480 
RE = 0.098017    IM = -0.995185 
RE = 0.073565    IM = -0.997290 
RE = 0.049068    IM = -0.998795 
RE = 0.024541    IM = -0.999699

And this is only a minimal exaple that I can explain you as much as possibble.

What i want and tried to create is a table that begins with the earlier set of lines and the resdt will be as following:

(earlier set of lines) expressed as following(idx:0 --- > 64)

    re[idx] = (float) cos((float)i * (2*pi)/(float)N);
    im[idx] = (float)-sin((float)i * (2*pi)/(float)N);

(2nd set of lines) expressed as following repeated 4 times consequently

    re[idx] = (float) cos(4 * (float)i * (2*pi)/(float)N);
    im[idx] = (float)-sin(4 * (float)i * (2*pi)/(float)N);

(3nd set of lines) expressed as following repeated 16 times consequently

    re[idx] = (float) cos(16 * (float)i * (2*pi)/(float)N);
    im[idx] = (float)-sin(16 * (float)i * (2*pi)/(float)N);

The result is expected to be as following:

//set 1 as above  repeated just 1 time
// ....
//set 2 repeated 4 times
RE = 1.000000    IM = -0.000000 //1st ligne of set 1
RE = 0.995185    IM = -0.098017 //4th ligne of set 1
RE = 0.980785    IM = -0.195090 //8th ligne of set 1
RE = 0.956940    IM = -0.290285 //12th ligne of set 1 ...
RE = 0.923880    IM = -0.382683 
RE = 0.881921    IM = -0.471397 
RE = 0.831470    IM = -0.555570 
RE = 0.773010    IM = -0.634393 
RE = 0.707107    IM = -0.707107 
RE = 0.634393    IM = -0.773010 
RE = 0.555570    IM = -0.831470 
RE = 0.471397    IM = -0.881921 
RE = 0.382683    IM = -0.923880 
RE = 0.290285    IM = -0.956940 
RE = 0.195090    IM = -0.980785 
RE = 0.098017    IM = -0.995185
// set 3 repeated 16 times
RE = 1.000000    IM = -0.000000 //1st ligne of set 1
RE = 0.923880    IM = -0.382683 //16th ligne of set 1
RE = 0.707107    IM = -0.707107 //38th ligne of set 1
RE = 0.382683    IM = -0.923880 //64th ligne of set 1

I've tried several times but i keep getting wrong results I don't know if it's precison issue or sometimng else.

Ok just a minute i will make an edit to show you the expected result — A.nechi, Jul 04 '17 at 16:17
Unrelated to your problem, but why don't you use [`cosf`](http://en.cppreference.com/w/c/numeric/math/cos) if you want to use `float`? — Some programmer dude, Jul 04 '17 at 16:18
Likely your problem is that FFT interleaving is very involved. I wrote some code to do it for arbitrary window sizes once, but I couldn't replicate it now without a lot of effort. — Malcolm McLean, Jul 04 '17 at 16:18
As for your problem, you want the 64 first entries to be the same as now? And then should it be four or 256 lines with "2nd set of lines"? And so on for the last "set"? Can you perhaps show us a example of the wanted output for a *smaller* set of data? — Some programmer dude, Jul 04 '17 at 16:21
4th set: 64 (previous * 4) or 256 (previous * previous) lines? — Aconcagua, Jul 04 '17 at 16:35
Please have a look at the edit that I've made ... i hope it's clearer now — A.nechi, Jul 04 '17 at 16:35
Yes I'm using the split data technique for more flexibility Because I'm trying to enhance time performance of an FFT ... so using a struct will make the access and loads a little bit heavy ...also using the complex type which is a 2 double type digits is impossible due to the limited hardware. — A.nechi, Jul 04 '17 at 17:35

Aconcagua · Accepted Answer · 2017-07-05T08:01:01.377

You might maintain factor and set size in an additional (outer) loop:

// you will need more than you calculated previously!
float twiddle_real[TWIDDLE_LIMIT * 3];
float twiddle_imag[TWIDDLE_LIMIT * 3];

// pointer arithmetic...
float* real = twiddle_real;
float* imag = twiddle_imag;

double factor = 1.0;
for(int size = TWIDDLE_LIMIT; size > 1; size /= 4)
{
    for(int j = 0; j < 64; ++j)
    {
        *real++ = (float)  cos((j % size) * factor * 2.0 * PI / N);
        *imag++ = (float) -sin((j % size) * factor * 2.0 * PI / N);
    }
    factor *= 4.0;
}

You do not need all those casts by the way - since factor is a double, (j % size) is converted to implicitly, as well as N afterwards.

Recommendation: As re and img belong together, I would represent them as such:

struct Complex
{
    double re;
    double im;
};

You could then have an array of these:

struct Complex twiddle[TW_size]; // no need for malloc, by the way...

You might have noticed: I changed to double, too. There is no reason to use float (with more limited precision) unless you have limited memory available (micro controllers)...

Alternative (why re-inventing the wheel?): use complex.h.

Gene · Answer 2 · 2017-07-04T18:31:24.893

It doesn't work because your math is wrong. The first table covers only the first quadrant of the complex plane. When you multiply by 4 for the second segment of your desired output, you won't get 4 repeats of the first quadrant, you'll get all four quadrants.

Here is something closer to what you specified. And it fixes some problems with your use of C:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N_ROOTS_OF_UNITY 64
#define PI 3.14159265358979323846264338327950

float *twiddle_real;
float *twiddle_imag;

int main()
{
  size_t table_size = N_ROOTS_OF_UNITY * 3;

  twiddle_real  = malloc(table_size * sizeof(float));
  twiddle_imag  = malloc(table_size * sizeof(float));

  size_t i, incr, repeat, k = 0;
  for (incr = 1; incr <= 16; incr *= 4) {
    for (repeat = 0; repeat < incr; ++repeat) {
      for(i = 0; i < N_ROOTS_OF_UNITY; i += incr) {
        twiddle_real[k] = (float) cos(i * (PI / 2) / N_ROOTS_OF_UNITY);
        twiddle_imag[k] = (float) -sin(i * (PI / 2) / N_ROOTS_OF_UNITY);
        ++k;
      }
    }
  }
  for (int a = 0; a < table_size; a++)
    printf("%d: RE = %f\tIM = %f\n", a, twiddle_real[a], twiddle_imag[a]);
  return 0;
}

Create a twiddle table with a specific order?

2 Answers2