1
$\begingroup$

I have a small issue.

I'm going to do FFT2 onto a matrix $X$, and the only library I have is FFTPack. FFTPack is the oldest FFT library, but it's one of the fastest and 100% portable too.

The problem with FFTPack is that it has no FFT2 function. Only FFTPack version 5.1 has FFT2 functions, but those functions are written in Fortran 77. My FFTPack is written in ANSI C (C89) code and is a little bit older library than version 5.1. So I missing FFT2 functionality. Written in Fortran 77 and converted to ANSI C code. The library was first made in K&R-standard from 1982.

Anyway. When I do FFT on each row of the matrix $X$, then I get the same result as I do i MATLAB.

fft(X')

Problem:

The values from every row of $X$ when I do FFT, they are complex numbers. My FFTPack function cannot handle complex numbers because my compiler does not support C99 standard and above. Only C89 standard (Yes, I know....).

So is there any way to do FFT again onto the complex numbers by separate them so imaginary numbers will be in one part and the real numbers will be inside another part?

Here is the FFTPack 5.1 version: https://github.com/jlokimlin/fftpack5.1

$\endgroup$
12
  • $\begingroup$ So, port a more modern (and that is really not a hard to achieve requirement here) version of fftpack, or use a different library (fits, maybe) or learn to link against the original Fortran fftpack? $\endgroup$ Sep 3, 2023 at 11:48
  • $\begingroup$ And use a modern compiler. It worst case, transpile to c89. $\endgroup$ Sep 3, 2023 at 11:50
  • 2
    $\begingroup$ @euraad: why don't you just use a different library (FFTW for example) or roll your own? The core FFT algorithm isn't all that complicated and there are plenty of code examples out there. Most libraries have way more stuff in them than most people every need so I can restrict a port or write to your specific requirements. $\endgroup$
    – Hilmar
    Sep 3, 2023 at 13:20
  • 2
    $\begingroup$ What @Hilmar said. Probably everyone and their dog has a simple radix-2 Cooley-Tukey FFT in source code. I thought I might have posted mine here, but I don't see it. Lemme know if you want me to find my 4 decade old FFT code written in C. I know I posted it to comp.dsp. $\endgroup$ Sep 3, 2023 at 14:59
  • 1
    $\begingroup$ That. I meant to mention "FFTs" in my first comment, but my phone corrected that to "fits". Anyways, I mean, really, get a compiler from the last 30 years. Compilers had _Complex long before it appeared in C99. Of course, pre-standard implementations suffer lack of compatibility and portability (this very much bit us in a vector math library we're maintaining!), so, again, staying in a c89 compiler is generally undesirable for this very problem you're solving. $\endgroup$ Sep 3, 2023 at 15:04

2 Answers 2

2
$\begingroup$

Radix-2 FFT written in ANSI C. This is really ducking old.

/*  
    A set of utility programs to compute the Fast Fourier Transform (FFT):
    
                        N-1
                    X[k] = SUM { x[n]exp(-j2(pi)nk/N) }
                        n=0
    
    and inverse Fast Fourier Transform (iFFT):
    
                            N-1
                    x[n] = 1/N SUM { X[k]exp(+j2(pi)nk/N) }
                            k=0
    
    To speed things up, multiplication by 1 and j are avoided.  The input, x[],
    is an array of complex numbers (pairs of doubles) of length N = 2^p.  The
    calling program supplies p = log2(N) not the array length, N.  The output
    is placed in BIT REVERSED order in x[].  Call bitreverse(x, p) to swap back
    to normal order, if needed.  However, iFFT(X, p, trigtbl) requires its input,
    X[], to be in bit reversed order making bit reversing unnecessary in some
    cases, such as fast convolution.  trigtbl is an array of doubles of length
    N/4 containing the sin function from 0 to pi/2 used to compute the FFT.
    Call sintable(trigtbl, p) ONLY ONCE before either FFT(x, p, trigtbl) or
    iFFT(X, p, trigtbl) to set up a sin table for FFT computation.
    
    Written in Megamax C (for the Mac) by Robert Bristow-Johnson (1985).
*/

/*
#include <complex.h>
*/

typedef struct {
    double real;
    double imag;
} complex;

#define Re(z) (z).real
#define Im(z) (z).imag

#define PI (3.1415926535897932587)

double sin();

FFT (x, p, trigtbl)
complex x[];
int p;
register double *trigtbl;
{
    register long length, step, stepsize, size;
    register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
    complex temp;
    
    size = 1L<<p;
    end = x + size;
    
    length = size>>1;   
    size >>= 2;
    stepsize = 1L;
    while ( length > 1L) {
        top = x;
        while (top < end) {
            bottom = top + length;
            
            Re(temp) = Re(*top) - Re(*bottom);          /* butterfly: twiddle = 1 */
            Im(temp) = Im(*top) - Im(*bottom);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            *bottom = temp;
            top++;
            bottom++;
            
            for (step = stepsize; step < size; step += stepsize) {
                Re(temp) = Re(*top) - Re(*bottom);      /* butterfly: twiddle = exp(-j�) */
                Im(temp) = Im(*top) - Im(*bottom);
                Re(*top) += Re(*bottom);
                Im(*top) += Im(*bottom);
                Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
                Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
                top++;
                bottom++;
            }
            
            Re(temp) = Im(*top) - Im(*bottom);          /* butterfly: twiddle = -j */
            Im(temp) = Re(*bottom) - Re(*top);
            Re(*top) += Re(*bottom);
            Im(*top) += Im(*bottom);
            *bottom = temp;
            top++;
            bottom++;
            
            for (step = stepsize; step < size; step += stepsize) {
                Re(temp) = Im(*top) - Im(*bottom);      /* butterfly: twiddle = -j*exp(-j�) */
                Im(temp) = Re(*bottom) - Re(*top);
                Re(*top) += Re(*bottom);
                Im(*top) += Im(*bottom);
                Re(*bottom) = Re(temp)*trigtbl[size - step] + Im(temp)*trigtbl[step];
                Im(*bottom) = Im(temp)*trigtbl[size - step] - Re(temp)*trigtbl[step];
                top++;
                bottom++;
            }
            top = bottom;
        }
        stepsize <<= 1;
        length >>= 1;
    }
    
    top = x;
    bottom = x + 1;
    while (top <  end) {
        Re(temp) = Re(*top) - Re(*bottom);              /* butterfly: twiddle = 1 */
        Im(temp) = Im(*top) - Im(*bottom);
        Re(*top) += Re(*bottom);
        Im(*top) += Im(*bottom);
        *bottom = temp;
        top += 2;
        bottom += 2;
    }
}


iFFT (X, p, trigtbl)
complex X[];
int p;
register double *trigtbl;
{
    register long length, step, stepsize, size;
    double scale;
    register complex *top, *bottom, *end;               /* top & bottom of FFT butterfly */
    complex temp;
    
    size = 1L<<p;
    end = X + size;
    
    scale = 1.0/(double)size;
    top = X;
    bottom = X + 1;
    while (top <  end) {
        Re(temp) = (Re(*top) - Re(*bottom))*scale;      /* butterfly: twiddle = 1/N */
        Im(temp) = (Im(*top) - Im(*bottom))*scale;
        Re(*top) = (Re(*top) + Re(*bottom))*scale;
        Im(*top) = (Im(*top) + Im(*bottom))*scale;
        *bottom = temp;
        top += 2;
        bottom += 2;
    }
    
    length = 1L;
    size >>= 2;
    stepsize = size;
    while ( stepsize >= 1L) {
        length <<= 1;
        top = X;
        while (top < end) {
            bottom = top + length;
            
            temp = *bottom;                             /* butterfly: twiddle = 1 */
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
            
            for (step = stepsize; step < size; step += stepsize) {
                                                        /* butterfly: twiddle = exp(+j�) */
                Re(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
                Im(temp) = Im(*bottom)*trigtbl[size - step] + Re(*bottom)*trigtbl[step];
                Re(*bottom) = Re(*top) - Re(temp);
                Im(*bottom) = Im(*top) - Im(temp);
                Re(*top) += Re(temp);
                Im(*top) += Im(temp);
                top++;
                bottom++;
            }
            
            Re(temp) = -Im(*bottom);                    /* butterfly: twiddle = +j */
            Im(temp) = Re(*bottom);
            Re(*bottom) = Re(*top) - Re(temp);
            Im(*bottom) = Im(*top) - Im(temp);
            Re(*top) += Re(temp);
            Im(*top) += Im(temp);
            top++;
            bottom++;
            
            for (step = stepsize; step < size; step += stepsize) {
                                                        /* butterfly: twiddle = +j*exp(+j�) */
                Re(temp) = -Im(*bottom)*trigtbl[size - step] - Re(*bottom)*trigtbl[step];
                Im(temp) = Re(*bottom)*trigtbl[size - step] - Im(*bottom)*trigtbl[step];
                Re(*bottom) = Re(*top) - Re(temp);
                Im(*bottom) = Im(*top) - Im(temp);
                Re(*top) += Re(temp);
                Im(*top) += Im(temp);
                top++;
                bottom++;
            }
            top = bottom;
        }
        stepsize >>= 1;
    }
}


sintable(trigtbl, p)
register double *trigtbl;
int p;
{
    register long size, i;
    double theta;
    
    size = 1L<<(p-2);
    theta = (PI/2.0)/(double)size;
    
    for (i = 0; i < size; i++) 
        trigtbl[i] = sin(theta*(double)i);
}


bitreverse (x, p)
register complex *x;
int p;
{
    complex temp;
    register long k, i, r, size, count;
    
    size = (1L<<p) - 1L;
    for (k = 1L; k < size; k++) {
        i = k;
        r = 0;
        for (count = size; count > 0; count >>= 1) {
            r <<= 1;
            r += (i & 0x00000001L);
            i >>= 1;
        }
        if (r > k) {
            temp = x[r];
            x[r] = x[k];
            x[k] = temp;
        }
    }
}
$\endgroup$
6
  • 2
    $\begingroup$ Wow, this really is ducking old! Love the K&R variable declaring between function stub and function body! $\endgroup$ Sep 3, 2023 at 15:07
  • 1
    $\begingroup$ Ya know, I think you're right, @MarcusMüller . This is K&R C and precedes ANSI C. It was written very soon after I was married in 1985. (I know it's after because my name is hyphenated.) $\endgroup$ Sep 3, 2023 at 15:15
  • $\begingroup$ :) dating code that way definitely is cool; anyway, any compiler that actually supports ANSI C but not yet C99 will work happily with that code. $\endgroup$ Sep 3, 2023 at 15:19
  • $\begingroup$ Them Aussies just dunno how to cuss down there down under. $\endgroup$ Sep 4, 2023 at 4:01
  • 1
    $\begingroup$ Old school C code is so cool :D $\endgroup$
    – euraad
    Sep 4, 2023 at 10:13
1
$\begingroup$

Up until a few years ago, MS Visual Studio supported C only as C89. Which was a major pain in bigger projects locked into that IDE/compiler and multi-platform projects.

If possible, I would use a modern compiler and a modern library.

If fft size is moderate and speed is of little importance, you could just implement a complex DFT as a plain matrix multiplication using your own definition of complex on top of real numbers in C89?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.