ifft.alg

Audacity是一款用於錄音和編輯聲音的、免費的開放源碼軟體。它可以執行於Mac OS X、Microsoft Windows、GNU/Linux和其它作業系統

(IFFT-ALG
  (NAME "ifft")
  (ARGUMENTS ("time_type" "t0") ("rate_type" "sr")
             ("LVAL" "src") ("long" "stepsize")
             ("LVAL" "window"))
  (SUPPORT-FUNCTIONS "
/* index: index into outbuf whree we get output samples
 * length: size of the frame, window, and outbuf; half size of samples
 * array: spectral frame goes here (why not a local var?)
 * window_len: size of window, should equal length
 * outbuf: real part of samples are multiplied by window and added to
 *          outbuf (after shifting)
 * src: send :NEXT to this object to get next frame
 * stepsize: shift by this many and add each frame
 * samples: result of ifft goes here, real and imag
 * window: multiply samples by window if any 
 *
 * IMPLEMENTATION NOTE:
 * The src argument is an XLisp object that returns either an
 * array of samples or NIL. The output of ifft is simply the
 * concatenation of the samples taken from the array. Later,
 * an ifft will be plugged in and this will return overlapped
 * adds of the ifft's.
 *
 * OVERLAP: stepsize must be less than or equal to the length
 * of real part of the transformed spectrum. A transform step
 * works like this: 
 * (1) shift the output buffer by stepsize samples, filling
 *     the end of the buffer with zeros
 * (2) get and transform an array of spectral coefficients
 * (3) multiply the result by a window
 * (4) add the result to the output buffer
 * (5) output the first stepsize samples of the buffer
 * 
 * DATA FORMAT: the DC component goes in array elem 0
 * Cosine part is in elements 2*i-1
 * Sine part is in elements 2*i
 * Nyquist frequency is in element length-1
 */

#include \"samples.h\"
#include \"fftn.h\"

table_type get_window_samples(LVAL window, sample_type **samples, long *len)
{
    table_type result = NULL;
    if (soundp(window)) {
        sound_type window_sound = getsound(window);
        xlprot1(window); /* maybe not necessary */
        result = sound_to_table(window_sound);
        xlpop();
        *samples = result->samples;
        *len = (long) (result->length + 0.5);
    }
    return result;
}
")

  (SAMPLE-RATE "sr")
  (STATE
          ("long" "index" "stepsize") ; samples index
          ("long" "length" "0")       ; samples length
          ("LVAL" "array" "NULL")
          ("long" "window_len" "0")
          ("sample_type *" "outbuf" "NULL")
          ("LVAL" "src" "src")
          ("long" "stepsize" "stepsize")
          ("sample_type *" "window" "NULL") ; window samples
          ("sample_type *" "samples" "NULL")
          ("table_type" "table" 
           "get_window_samples(window, &susp->window, &susp->window_len)"))
  
  (OUTER-LOOP "
        if (susp->src == NULL) {
out:        togo = 0;   /* indicate termination */
            break;      /* we're done */
        }
        if (susp->index >= susp->stepsize) {
            long i;
            long half_i;
            long n;
            LVAL elem;
            susp->index = 0;
            susp->array = 
                xleval(cons(s_send, cons(susp->src, consa(s_next))));
            if (susp->array == NULL) {
                susp->src = NULL;
                goto out;
            } else if (!vectorp(susp->array)) {
                xlerror(\"array expected\", susp->array);
            } else if (susp->samples == NULL) {
                /* assume arrays are all the same size as first one;
                   now that we know the size, we just have to do this
                   first allocation.
                 */
                susp->length = getsize(susp->array);
                if (susp->length < 1) 
                    xlerror(\"array has no elements\", susp->array);
                if (susp->window && (susp->window_len != susp->length))
                    xlerror(\"window size and spectrum size differ\", 
                            susp->array);
                susp->samples = 
                    (sample_type *) calloc(susp->length * 2,
                                           sizeof(sample_type));
                susp->outbuf = (sample_type *) calloc(susp->length, 
                    sizeof(sample_type));
            } else if (getsize(susp->array) != susp->length) {
                xlerror(\"arrays must all be the same length\", susp->array);
            }

            /* at this point, we have a new array to put samples */
            /* real part will be susp->samples[0:n-1], */
            /* im part in samples[n:2*n-1] */
            n = susp->length;
            elem = getelement(susp->array, 0);
            if (ntype(elem) != FLONUM) {
                xlerror(\"flonum expected\", elem);
            }
            susp->samples[0] = (sample_type) getflonum(elem);
            susp->samples[n] = 0;
            half_i = 0;
            for (i = 1; i < n - 1; i += 2) {
                half_i++;
                elem = getelement(susp->array, i);
                if (ntype(elem) != FLONUM) {
                    xlerror(\"flonum expected\", elem);
                }
                susp->samples[half_i] = susp->samples[n - half_i] = 
                    (sample_type) (getflonum(elem) / 2.0);

                elem = getelement(susp->array, i + 1);
                if (ntype(elem) != FLONUM) {
                    xlerror(\"flonum expected\", elem);
                }
                susp->samples[n + half_i] =
                    -(susp->samples[2*n - half_i] =
                          (sample_type) (getflonum(elem) / 2.0));
            }
            if (n % 2 == 0) {
                elem = getelement(susp->array, n - 1);
                if (ntype(elem) != FLONUM) {
                    xlerror(\"flonum expected\", elem);
                }
                susp->samples[n / 2] = (sample_type) getflonum(elem);
                susp->samples[n + (n / 2)] = 0;
            }
            susp->array = NULL; /* free the array */

            /* here is where the IFFT and windowing should take place */
            fftnf(1, &n, susp->samples, susp->samples + n, -1, 1.0);
            if (susp->window) {
                n = susp->length;
                for (i = 0; i < n; i++) {
                    susp->samples[i] *= susp->window[i];
                }
            }

            /* shift the outbuf */
            n = susp->length - susp->stepsize;
            for (i = 0; i < n; i++) {
                susp->outbuf[i] = susp->outbuf[i + susp->stepsize];
            }

            /* clear end of outbuf */
            for (i = n; i < susp->length; i++) {
                susp->outbuf[i] = 0;
            }

            /* add in the ifft result */
            n = susp->length;
            for (i = 0; i < n; i++) {
                susp->outbuf[i] += susp->samples[i];
            }
/*
            temp_fft = (double *) malloc (susp->length * sizeof(double));
            if (temp_fft == 0) return;
            big_samples = (double *) malloc (susp->length * sizeof(double));
            if (big_samples == 0) return;
            for (i = 0; i < susp->length; i++) {
                big_samples[i] = (double) susp->samples[i];
            }
            rp = rfftw_create_plan(susp->length, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
            rfftw_one(rp, big_samples, temp_fft);
            rfftw_destroy_plan(rp);
            free(big_samples);
            for (i = 0; i < susp->length; i++) {
                setelement(result, i, cvflonum(temp_fft[i]));
            }
            free (temp_fft);
*/
        }
        togo = MIN(togo, susp->stepsize - susp->index);
")
  (INNER-LOOP "output = outbuf[index++];")
  (CONSTANT "length" "samples" "array" "src" "window")
  (TERMINATE COMPUTED)
  (FINALIZATION "    if (susp->samples) free(susp->samples);
    if (susp->table) table_unref(susp->table);
    if (susp->outbuf) free(susp->outbuf);
")

)