Ron Mayer FFT convolution algorithm

Question

Does anybody know anything about the Mayer FFT implementation (without me having to take quite a bit of time to study the code)?

I am attempting to perform a convolution, and the ifft appears to produce what I would call a "mirrored" output. In other words, my kernel+signal length is limited to N/2 and whatever occupies n=0...N/2 is mirrored about n=N...N/2. It sort of looks like what I would expect from an FFT in terms of negative frequencies...except it's like a mirror in negative time.

Here is my convolution code:

   void convolve(struct cxType* data,  struct cxType* kernel, int size) 
    {
    int i,j;
    int wrksz = size;
    float gain = 1.0f/((float) wrksz);


    mayer_fft(wrksz, data->re, data->im);
    mayer_fft(wrksz, kernel->re, kernel->im);

    for(i=0;i<wrksz;i++)
    {
    data->re[i]*=kernel->re[i]*gain;
    data->im[i]*=kernel->im[i]*gain;
    }

    mayer_ifft(wrksz, data->re, data->im);
}

Doing essentially the same thing using gnu Octave (MATLAB syntax equivalent for you who are unfamiliar) produces the expected results, including allowing me to occupy M+N-1 in the signal output:

fs=48000;
Ts = 1/fs;
NN =  1024
sincsz = floor(NN/4);
sigstart = floor(NN/16);
sigend = floor(NN/2);
dpi=2*pi;

%window func
tau=(1:sincsz)/sncsz;
window=0.5*(1.0 - cos(dpi*tau));
%plot(tau,window)

%sinc filter kernel
fc=5050;
wc=dpi*fc;
Ts=1/fs;
segTime=Ts*sincsz;
t0=-segTime/2;
t=Ts*((1:sincsz) - 1) + t0 ;

s=zeros(1,sincsz);
s=window.*sin(wc*t)./(wc*t);
s(sincsz/2+1) = 1;

%plot(t,s)
fund = 1650;
tt = (1:NN)*Ts;
signal = sin(dpi*tt*fund) + sin(dpi*tt*2*fund) + sin(dpi*tt*3*fund) + sin(dpi*tt*4*fund) + sin(dpi*tt*5*fund);
signal(1:sigstart) = signal(1:sigstart)*0;
signal(sigend:NN) = signal(sigend:NN)*0;
%plot(tt,signal)

h=zeros(1,NN);
h(1:sincsz) = s(1:sincsz);

H=fft(h);
X=fft(signal);

Y=H.*X;

y=ifft(Y);

plot(real(y))

The equivalent signal and FIR kernel synthesis are implemented in C (not shown). I am using gnuplot to display the results of the C implementation of this, so I know the filter kernel and the signal are implemented identically to what I did with Octave.

The only difference in what is being done as far as I can tell is the FFT implementation. Does somebody know if these results are a fundamental problem with my understanding of the FFT algorithm in general, or with the FHT-based implementation from this archaic code written by Ron Mayer? You can go to his archived web site to get the code I am using: http://reocities.com/researchtriangle/8869/fft_summary.html

Now, if I perform an FFT on a block of data and then perform the ifft, I get back the original data as I would expect. If I modify the data spectrally in any way, I get results that I don't expect.

I once tried to use this mayer FFT algorithm to replace that used in S. Bernsee's pitch shift algorithm, and this did not work at all. I used fftw3, and that code worked as I expected. I am curious to try this same basic algorithm with fftw3 to see what happens.

What I don't know is if I misunderstand something fundamental which is causing me to mis-apply rmayer's implementation, or if this is just a simple artifact I have to work around (that is, using an FFT size double of what I expect).

RJB RJB · Accepted Answer · 2013-04-19T03:40:35

Doh! This is one of those things like forgetting to put a semicolon at the end of the line. A convolution is a complex multiplication in the frequency domain -- I was doing a naive point-by-point multiplication. Here is the corrected code, which shows the complex multiplication. Of course, there are C and C++ constructs and macros/routines for doing this, but here is the brute force method for pedagogical purposes: Assume struct cxType is defined as:

struct cxType {
float* re;
float* im;
};  //and such a struct should have mem allocated before sending it into convolve()

    void convolve(struct cxType* data,  struct cxType* kernel, int size) 
        {
        int i,j;
        int wrksz = size;
        float gain = 1.0f/((float) wrksz);
        float a,b,c,d;


        mayer_fft(wrksz, data->re, data->im);
        mayer_fft(wrksz, kernel->re, kernel->im);

        for(i=0;i<wrksz;i++)
        {
        a=data->re[i];
        b=data->im[i];
        c=kernel->re[i];
        d=kernel->im[i];
        data->re[i]=(a*c - b*d)*gain;
        data->im[i]=(a*d + b*c)*gain;
        }

        mayer_ifft(wrksz, data->re, data->im);
    }

And now the above code snippet works when it is called and does not produce these strange mirroring effects I was talking about. As for my Matlab code, it works because Octave/Matlab hide the detail of complex multiplication behind their handy syntax, H.*X.

I was able to reproduce my problem in Octave by multiplying and combining the real and imaginary parts separately to mimic my mistake in C, and then fixed it by doing something similar to above. That said, there is nothing unique about rmayer's implementation of the FFT...only my implementation of convolution was amiss.

Ron Mayer FFT convolution algorithm

1 Answers