plugins/fourier/fourier.C

   1 #include <math.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4
   5 #include "clip.h"
   6 #include "fourier.h"
   7
   8 #define WINDOW_BORDER (window_size / 2)
   9
  10
  11 FFT::FFT()
  12 {
  13 }
  14
  15 FFT::~FFT()
  16 {
  17 }
  18
  19 int FFT::do_fft(unsigned int samples,  // must be a power of 2
  20         int inverse,         // 0 = forward FFT, 1 = inverse
  21         double *real_in,     // array of input's real samples
  22         double *imag_in,     // array of input's imag samples
  23         double *real_out,    // array of output's reals
  24         double *imag_out)
  25 {
  26     unsigned int num_bits;    // Number of bits needed to store indices
  27     unsigned int i, j, k, n;
  28     unsigned int block_size, block_end;
  29
  30     double angle_numerator = 2.0 * M_PI;
  31     double tr, ti;     // temp real, temp imaginary
  32
  33     if(inverse)
  34         angle_numerator = -angle_numerator;
  35
  36     num_bits = samples_to_bits(samples);
  37
  38 // Do simultaneous data copy and bit-reversal ordering into outputs
  39
  40     for(i = 0; i < samples; i++)
  41     {
  42         j = reverse_bits(i, num_bits);
  43         real_out[j] = real_in[i];
  44         imag_out[j] = (imag_in == 0) ? 0.0 : imag_in[i];
  45     }
  46
  47 // Do the FFT itself
  48
  49     block_end = 1;
  50     double delta_angle;
  51     double sm2;
  52     double sm1;
  53     double cm2;
  54     double cm1;
  55     double w;
  56     double ar[3], ai[3];
  57     double temp;
  58     for(block_size = 2; block_size <= samples; block_size <<= 1)
  59     {
  60         delta_angle = angle_numerator / (double)block_size;
  61         sm2 = sin(-2 * delta_angle);
  62         sm1 = sin(-delta_angle);
  63         cm2 = cos(-2 * delta_angle);
  64         cm1 = cos(-delta_angle);
  65         w = 2 * cm1;
  66
  67         for(i = 0; i < samples; i += block_size)
  68         {
  69             ar[2] = cm2;
  70             ar[1] = cm1;
  71
  72             ai[2] = sm2;
  73             ai[1] = sm1;
  74
  75             for(j = i, n = 0; n < block_end; j++, n++)
  76             {
  77                 ar[0] = w * ar[1] - ar[2];
  78                 ar[2] = ar[1];
  79                 ar[1] = ar[0];
  80
  81                 ai[0] = w * ai[1] - ai[2];
  82                 ai[2] = ai[1];
  83                 ai[1] = ai[0];
  84
  85                 k = j + block_end;
  86                 tr = ar[0] * real_out[k] - ai[0] * imag_out[k];
  87                 ti = ar[0] * imag_out[k] + ai[0] * real_out[k];
  88
  89                 real_out[k] = real_out[j] - tr;
  90                 imag_out[k] = imag_out[j] - ti;
  91
  92                 real_out[j] += tr;
  93                 imag_out[j] += ti;
  94             }
  95         }
  96
  97         block_end = block_size;
  98     }
  99
 100 // Normalize if inverse transform
 101
 102     if(inverse)
 103     {
 104         double denom = (double)samples;
 105
 106         for (i = 0; i < samples; i++)
 107         {
 108             real_out[i] /= denom;
 109             imag_out[i] /= denom;
 110         }
 111     }
 112         return 0;
 113 }
 114
 115
 116 unsigned int FFT::samples_to_bits(unsigned int samples)
 117 {
 118     unsigned int i;
 119
 120     for(i = 0; ; i++)
 121     {
 122         if(samples & (1 << i))
 123             return i;
 124     }
 125         return i;
 126 }
 127
 128 unsigned int FFT::reverse_bits(unsigned int index, unsigned int bits)
 129 {
 130     unsigned int i, rev;
 131
 132     for(i = rev = 0; i < bits; i++)
 133     {
 134         rev = (rev << 1) | (index & 1);
 135         index >>= 1;
 136     }
 137
 138     return rev;
 139 }
 140
 141 int FFT::symmetry(int size, double *freq_real, double *freq_imag)
 142 {
 143     int h = size / 2;
 144     for(int i = h + 1; i < size; i++)
 145     {
 146         freq_real[i] = freq_real[size - i];
 147         freq_imag[i] = -freq_imag[size - i];
 148     }
 149         return 0;
 150 }
 151
 152
 153
 154
 155
 156 CrossfadeFFT::CrossfadeFFT() : FFT()
 157 {
 158         reset();
 159         window_size = 4096;
 160 }
 161
 162 CrossfadeFFT::~CrossfadeFFT()
 163 {
 164         delete_fft();
 165 }
 166
 167 int CrossfadeFFT::reset()
 168 {
 169         input_buffer = 0;
 170         output_buffer = 0;
 171         freq_real = 0;
 172         freq_imag = 0;
 173         temp_real = 0;
 174         temp_imag = 0;
 175         first_window = 1;
 176         input_size = 0;          // samples in input_buffer and output_buffer
 177         output_size = 0;
 178         input_allocation = 0;
 179         output_allocation = 0;
 180         return 0;
 181 }
 182
 183 int CrossfadeFFT::delete_fft()
 184 {
 185         if(input_buffer) delete [] input_buffer;
 186         if(output_buffer) delete [] output_buffer;
 187         if(freq_real) delete [] freq_real;
 188         if(freq_imag) delete [] freq_imag;
 189         if(temp_real) delete [] temp_real;
 190         if(temp_imag) delete [] temp_imag;
 191         reset();
 192         return 0;
 193 }
 194
 195 int CrossfadeFFT::fix_window_size()
 196 {
 197 // fix the window size
 198 // window size must be a power of 2
 199         int new_size = 16;
 200         while(new_size < window_size) new_size *= 2;
 201         window_size = MIN(131072, window_size);
 202         window_size = new_size;
 203         return 0;
 204 }
 205
 206 int CrossfadeFFT::initialize(int window_size)
 207 {
 208         this->window_size = window_size;
 209         first_window = 1;
 210         reconfigure();
 211         return 0;
 212 }
 213
 214 long CrossfadeFFT::get_delay()
 215 {
 216         return window_size + WINDOW_BORDER;
 217 }
 218
 219 int CrossfadeFFT::reconfigure()
 220 {
 221         delete_fft();
 222         fix_window_size();
 223
 224
 225
 226         return 0;
 227 }
 228
 229 int CrossfadeFFT::process_fifo(long size, double *input_ptr, double *output_ptr)
 230 {
 231 //printf("CrossfadeFFT::process_fifo 1\n");
 232 // Load next input buffer
 233         if(input_size + size > input_allocation)
 234         {
 235                 double *new_input = new double[input_size + size];
 236 //printf("CrossfadeFFT::process_fifo 1\n");
 237                 if(input_buffer)
 238                 {
 239                         memcpy(new_input, input_buffer, sizeof(double) * input_size);
 240                         delete [] input_buffer;
 241                 }
 242 //printf("CrossfadeFFT::process_fifo 1\n");
 243                 input_buffer = new_input;
 244                 input_allocation = input_size + size;
 245 //printf("CrossfadeFFT::process_fifo 1\n");
 246         }
 247
 248 //printf("CrossfadeFFT::process_fifo 1\n");
 249         memcpy(input_buffer + input_size,
 250                 input_ptr,
 251                 size * sizeof(double));
 252 //printf("CrossfadeFFT::process_fifo 1\n");
 253         input_size += size;
 254
 255 //printf("CrossfadeFFT::process_fifo 1\n");
 256
 257
 258
 259
 260
 261
 262 // Have enough to do some windows
 263         while(input_size >= window_size)
 264         {
 265                 if(!freq_real) freq_real = new double[window_size];
 266                 if(!freq_imag) freq_imag = new double[window_size];
 267                 if(!temp_real) temp_real = new double[window_size];
 268                 if(!temp_imag) temp_imag = new double[window_size];
 269
 270
 271
 272                 do_fft(window_size,  // must be a power of 2
 273                 0,         // 0 = forward FFT, 1 = inverse
 274                 input_buffer,     // array of input's real samples
 275                 0,     // array of input's imag samples
 276                 freq_real,    // array of output's reals
 277                 freq_imag);
 278
 279                 int result = signal_process();
 280
 281                 if(!result)
 282                 {
 283                         do_fft(window_size,  // must be a power of 2
 284                         1,               // 0 = forward FFT, 1 = inverse
 285                         freq_real,     // array of input's real samples
 286                         freq_imag,     // array of input's imag samples
 287                         temp_real,     // array of output's reals
 288                         temp_imag);
 289                 }
 290
 291
 292 // Crossfade into the output buffer
 293                 long new_allocation = output_size + window_size;
 294                 if(new_allocation > output_allocation)
 295                 {
 296                         double *new_output = new double[new_allocation];
 297
 298                         if(output_buffer)
 299                         {
 300                                 memcpy(new_output, output_buffer, sizeof(double) * output_size);
 301 //printf("CrossfadeFFT::process_fifo 1 %p\n", output_buffer);
 302                                 delete [] output_buffer;
 303 //printf("CrossfadeFFT::process_fifo 2\n");
 304                         }
 305                         output_buffer = new_output;
 306                         output_allocation = new_allocation;
 307                 }
 308
 309                 if(output_size >= WINDOW_BORDER)
 310                 {
 311                         for(int i = 0, j = output_size - WINDOW_BORDER;
 312                                 i < WINDOW_BORDER;
 313                                 i++, j++)
 314                         {
 315                                 double src_level = (double)i / WINDOW_BORDER;
 316                                 double dst_level = (double)(WINDOW_BORDER - i) / WINDOW_BORDER;
 317                                 output_buffer[j] = output_buffer[j] * dst_level + temp_real[i] * src_level;
 318                         }
 319
 320                         memcpy(output_buffer + output_size,
 321                                 temp_real + WINDOW_BORDER,
 322                                 sizeof(double) * (window_size - WINDOW_BORDER));
 323                         output_size += window_size - WINDOW_BORDER;
 324                 }
 325                 else
 326                 {
 327 // First buffer has no crossfade
 328                         memcpy(output_buffer + output_size,
 329                                 temp_real,
 330                                 sizeof(double) * window_size);
 331                         output_size += window_size;
 332                 }
 333
 334
 335 // Shift input buffer forward
 336                 for(int i = window_size - WINDOW_BORDER, j = 0;
 337                         i < input_size;
 338                         i++, j++)
 339                         input_buffer[j] = input_buffer[i];
 340                 input_size -= window_size - WINDOW_BORDER;
 341         }
 342
 343
 344 //printf("CrossfadeFFT::process_fifo 1\n");
 345
 346 //printf("CrossfadeFFT::process_fifo 1 %d %d\n", output_size - WINDOW_BORDER, size);
 347
 348
 349 // Have enough to send to output
 350         int samples_rendered = 0;
 351         if(output_size - WINDOW_BORDER >= size)
 352         {
 353                 memcpy(output_ptr, output_buffer, sizeof(double) * size);
 354                 for(int i = size, j = 0; i < output_size; i++, j++)
 355                         output_buffer[j] = output_buffer[i];
 356                 output_size -= size;
 357                 samples_rendered = size;
 358         }
 359         else
 360         {
 361                 bzero(output_ptr, sizeof(double) * size);
 362         }
 363 //printf("CrossfadeFFT::process_fifo 2\n");
 364
 365         return samples_rendered;
 366 }
 367
 368
 369
 370
 371