media/fft/fft-real-pair.c

   1 /*
   2  Free FFT and convolution (C)
   3
   4  Copyright (c) 2020 Project Nayuki. (MIT License)
   5  https://www.nayuki.io/page/free-small-fft-in-multiple-languages
   6
   7  Permission is hereby granted, free of charge, to any person obtaining a copy of
   8  this software and associated documentation files (the "Software"), to deal in
   9  the Software without restriction, including without limitation the rights to
  10  use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11  copies of the Software, and to permit persons to whom the Software is
  12  furnished to do so, subject to the following conditions:
  13  - The above copyright notice and this permission notice shall be included in
  14    all copies or substantial portions of the Software.
  15  - The Software is provided "as is", without warranty of any kind, express or
  16    implied, including but not limited to the warranties of merchantability,
  17    fitness for a particular purpose and noninfringement. In no event shall the
  18    authors or copyright holders be liable for any claim, damages or other
  19    liability, whether in an action of contract, tort or otherwise,
  20    arising from, out of or in connection with the Software or the use
  21    or other dealings in the Software.
  22 */
  23
  24 #include <math.h>
  25 #include <stdint.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28 #include "fft-real-pair.h"
  29
  30
  31 // Private function prototypes
  32 static size_t reverse_bits(size_t val, int width);
  33 static void *memdup(const void *src, size_t n);
  34
  35
  36 bool Fft_transform(float real[], float imag[], size_t n) {
  37     if (n == 0)
  38         return true;
  39     else if ((n & (n - 1)) == 0)  // Is power of 2
  40         return Fft_transformRadix2(real, imag, n);
  41 }
  42
  43
  44 bool Fft_inverseTransform(float real[], float imag[], size_t n) {
  45     return Fft_transform(imag, real, n);
  46 }
  47
  48
  49 /* looking to replace this with assembler, so limit it to n <= 8 */
  50 void Fft_transformRadixL8(float real[], float imag[],
  51                           float cos_table[], float sin_table[],
  52                           size_t n)
  53 {
  54     size_t size = (n / 2) * sizeof(float);
  55
  56     // Cooley-Tukey decimation-in-time radix-2 FFT
  57     for (size_t size = 2; size <= n; size *= 2) {
  58         size_t halfsize = size / 2;
  59         size_t tablestep = n / size;
  60         for (size_t i = 0; i < n; i += size) {
  61             for (size_t j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
  62                 size_t l = j + halfsize;
  63                 float tpre =  real[l] * cos_table[k] + imag[l] * sin_table[k];
  64                 float tpim = -real[l] * sin_table[k] + imag[l] * cos_table[k];
  65                 real[l] = real[j] - tpre;
  66                 imag[l] = imag[j] - tpim;
  67                 real[j] += tpre;
  68                 imag[j] += tpim;
  69             }
  70         }
  71         if (size == n)  // Prevent overflow in 'size *= 2'
  72             break;
  73     }
  74 }
  75
  76 bool Fft_transformRadix2(float real[], float imag[], size_t n) {
  77     // Length variables
  78     bool status = false;
  79     int levels = 0;  // Compute levels = floor(log2(n))
  80     for (size_t temp = n; temp > 1U; temp >>= 1)
  81         levels++;
  82     if ((size_t)1U << levels != n)
  83         return false;  // n is not a power of 2
  84
  85     // Trigonometric tables
  86     if (SIZE_MAX / sizeof(float) < n / 2)
  87         return false;
  88     size_t size = (n / 2) * sizeof(float);
  89     float *cos_table = malloc(size);
  90     float *sin_table = malloc(size);
  91     if (cos_table == NULL || sin_table == NULL)
  92         goto cleanup;
  93     for (size_t i = 0; i < n / 2; i++) {
  94         cos_table[i] = cos(2 * M_PI * i / n);
  95         sin_table[i] = sin(2 * M_PI * i / n);
  96     }
  97
  98     // Bit-reversed addressing permutation
  99     for (size_t i = 0; i < n; i++) {
 100         size_t j = reverse_bits(i, levels);
 101         if (j > i) {
 102             float temp = real[i];
 103             real[i] = real[j];
 104             real[j] = temp;
 105             temp = imag[i];
 106             imag[i] = imag[j];
 107             imag[j] = temp;
 108         }
 109     }
 110
 111     if (n <= 8) {
 112         Fft_transformRadixL8(real, imag, cos_table, sin_table, n);
 113     } else {
 114         // Cooley-Tukey decimation-in-time radix-2 FFT
 115         for (size_t size = 2; size <= n; size *= 2) {
 116             size_t halfsize = size / 2;
 117             size_t tablestep = n / size;
 118             for (size_t i = 0; i < n; i += size) {
 119                 for (size_t j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
 120                     size_t l = j + halfsize;
 121                     float tpre =  real[l] * cos_table[k] + imag[l] * sin_table[k];
 122                     float tpim = -real[l] * sin_table[k] + imag[l] * cos_table[k];
 123                     real[l] = real[j] - tpre;
 124                     imag[l] = imag[j] - tpim;
 125                     real[j] += tpre;
 126                     imag[j] += tpim;
 127                 }
 128             }
 129             if (size == n)  // Prevent overflow in 'size *= 2'
 130                 break;
 131         }
 132     }
 133     status = true;
 134
 135 cleanup:
 136     free(cos_table);
 137     free(sin_table);
 138     return status;
 139 }
 140
 141
 142
 143 static size_t reverse_bits(size_t val, int width) {
 144     size_t result = 0;
 145     for (int i = 0; i < width; i++, val >>= 1)
 146         result = (result << 1) | (val & 1U);
 147     return result;
 148 }
 149
 150
 151 static void *memdup(const void *src, size_t n) {
 152     void *dest = malloc(n);
 153     if (n > 0 && dest != NULL)
 154         memcpy(dest, src, n);
 155     return dest;
 156 }