diff options
author | Søren Sandmann Pedersen <ssp@redhat.com> | 2010-07-02 06:43:33 -0400 |
---|---|---|
committer | Søren Sandmann Pedersen <ssp@redhat.com> | 2010-07-02 06:43:33 -0400 |
commit | 186e2632569df76f940653566fe6d8785b60ebc9 (patch) | |
tree | 656c51422537d4fa627b5b661658fe3e4a339b30 /fft.c | |
parent | 9e4d0aa74e017b5f18766755e5d6421e053c1271 (diff) |
Shift to make it symmetric around sample 0
Diffstat (limited to 'fft.c')
-rw-r--r-- | fft.c | 144 |
1 files changed, 21 insertions, 123 deletions
@@ -62,7 +62,7 @@ static void fft (complex_t *buffer, int n, complex_t principal) { complex_t *temp, *even, *odd; - complex_t omega; + complex_t alpha; int i; if (n == 1) @@ -81,15 +81,15 @@ fft (complex_t *buffer, int n, complex_t principal) fft (even, n/2, complex_mul (principal, principal)); fft (odd, n/2, complex_mul (principal, principal)); - omega.re = 1.0; - omega.im = 0; + alpha.re = 1.0; + alpha.im = 0; for (i = 0; i < n/2; ++i) { - buffer[i] = complex_add (even[i], complex_mul (omega, odd[i])); - buffer[i + n/2] = complex_sub (even[i], complex_mul (omega, odd[i])); + buffer[i] = complex_add (even[i], complex_mul (alpha, odd[i])); + buffer[i + n/2] = complex_sub (even[i], complex_mul (alpha, odd[i])); - omega = complex_mul (omega, principal); + alpha = complex_mul (alpha, principal); } free (temp); @@ -109,7 +109,8 @@ main () for (i = 0; i < N_SAMPLES; ++i) { double sample_extent = (last - first) / N_SAMPLES; - double t = first + i * sample_extent + sample_extent / 2; + double sn = (i + N_SAMPLES/2) % N_SAMPLES - N_SAMPLES/2; + double t = sn * sample_extent; double iv = 0.0; double v; @@ -119,126 +120,10 @@ main () v = 1.0; else { -#if 0 - /* slightly stretched lanczos 2 */ - if (fabs (t) > 2.5) - v = 0.0; - else - v = (sin (4/5.0 * M_PI * t) / (4/5.0 * M_PI * t)) * pow ((sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2)), 1); -#endif - -#if 0 if (fabs (t) > 0.5) v = 0.0; else v = 1.0; -#endif - -#if 0 - /* lanczos 3 */ - if (fabs (t) > 4) - v = 0.0; - else - v = (sin (3.0/4 * M_PI * t) / (3.0/4 * M_PI * t)) * (sin (M_PI * 3.0/4 * t / 3) / (M_PI * 3.0/4 * t / 3)); -#endif - - /* gaussian window */ - if (fabs (t) > 4) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)) * (exp (-0.5 * t * t)); -#if 0 - /* raised cosine */ - if (fabs (t) > 2.5) - v = 0.0; - else - v = (sin (4/5.0 * M_PI * t) / (4/5.0 * M_PI * t)) * (pow (0.5 + 0.5 * cos (M_PI * t / 2.5), 0.55)); -#endif -#if 0 - v = pow (0.5 + 0.5 * cos (M_PI * t / 2.5), 0.55); - if (fabs (t) > 1.5) - v = 0.0; - else - v = 1.0; -#endif - -#if 0 - v = sin (2 * M_PI * t * 1/100000.0 * pow (drand48(), 0.5)); -#endif -#if 0 - int j; - - v = 0.0; - t = fmod ((drand48() * 1024 - 512.0), 17) * (512.0/17); - v = t; //512 * pow (fabs (t)/512.0, 0.5 * drand48()) * sin (2 * M_PI * (fabs (t)/512.0) * t); -#endif - -#if 0 - if (fabs (t) > 10) - v = 0.0; - else - v = (sin (4/4.5 * M_PI * t) / (4/4.5 * M_PI * t)) * (sin (4/4.5 * M_PI * t / 2.0) / (4/4.5 * M_PI * t / 2.0)); -#endif -#if 0 - if (fabs (t) < 0.01) - 0.01; -#endif - - double lol = drand48(); - -#if 0 - if (fabs (t) > 4) - v = 0; - else - v = sin (2 * M_PI * t); -#endif - -#if 0 -#define EXPO 0.6 - - v = (sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2)); - if (v < 0) - { - v = cos (M_PI * EXPO) * pow (-v, EXPO); - - iv = sin (M_PI * EXPO); - - } - else - { - v = pow (v, EXPO); - } -#endif -#if 0 - if (fabs (t) > 2.5) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)) * (0.5 * (1 - cos (((t*M_PI)/2.5 - M_PI)))); -#endif -#if 0 - if (fabs (t) > 1.5) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)) * (sin (M_PI * (t/1.5)) / M_PI / (t/1.5)); -#endif -#if 0 - if (fabs (t) > 2.5) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)) * (sin (M_PI * (t/2.5)) / M_PI / (t/2.5)); -#endif -#if 0 - if (fabs (t) > 3) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)) * (sin (M_PI * (t/3)) / M_PI / (t/3)); -#endif -#if 0 - if (fabs (t) > 2) - v = 0.0; - else - v = (sin (M_PI * t) / (M_PI * t)); -#endif } test[i].re = v; @@ -253,8 +138,21 @@ main () p.im = sin ( 2 * M_PI / N_SAMPLES); fft (test, N_SAMPLES, p); +#if 0 + complex_t *swapbuffer = malloc (N_SAMPLES * sizeof (complex_t)); + for (i = 0; i < N_SAMPLES; ++i) + { + swapbuffer[i] = test[i]; + + if (i % 2 == 1) + swapbuffer[i].re = -swapbuffer[i].re; + } for (i = 0; i < N_SAMPLES; ++i) + test[i] = swapbuffer[i]; +#endif + + for (i = 0; i < N_SAMPLES; ++i) { double im, re; int idx = (i + N_SAMPLES/2) % N_SAMPLES; |