2
$\begingroup$

I would like to obtain FIR filter (low-pass) coefficients from a given frequency response by using C++ to develop algorithms. I have read at http://iowahills.com/B2PolynomialFIRFilters.html, that, if done properly, almost any filter response can be implemented as an FIR. In principle, one simply does an Inverse Fourier Transform on the filter's frequency response. Is there a C or C++ code that will do this? Any links or detailed guidelines how to solve this would be appreciated! I want to thank member hotpow2 for suggesting gradient descent routine, but I need more details on that as well.

$\endgroup$
  • 1
    $\begingroup$ The method you're referring to is called frequency sampling method, and this question and its answers may be helpful. $\endgroup$ – Matt L. May 13 '15 at 18:33
2
$\begingroup$

here's some very old source in C to do an FFT. it's nowhere near as good as FFTW, but it's free and self-contained.

if you're doing this inverse DFT and windowing (i would recommend the Kaiser window) to design an FIR, probably it doesn't need to be an FFT. you can write a very simple DFT (no Cooley-Tukey radix-2 optimizations) and run it. i don't think you'll care too much if it takes a couple of seconds longer to run.

/*  
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.14159265359)

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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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