fft - Ron Mayer FFT 卷积算法

Question

有人知道有关 Mayer FFT 实现的任何信息吗(不需要我花大量时间来研究代码)？

我正在尝试执行卷积，并且 ifft 似乎产生了我所说的“镜像”输出。换句话说，我的内核+信号长度限制为 N/2，并且无论占用 n=0...N/2 ，都会镜像 n=N...N/2 。就负频率而言，它看起来有点像我对 FFT 的期望……只不过它就像负时间中的一面镜子。

这是我的卷积代码:

   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);
}

使用 gnu Octave(对于不熟悉的人来说，MATLAB 语法等效)做基本上相同的事情会产生预期的结果，包括允许我在信号输出中占据 M+N-1:

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))

等效信号和 FIR 内核综合是用 C 语言实现的(未显示)。我使用 gnuplot 来显示 C 实现的结果，因此我知道滤波器内核和信号的实现与我使用 Octave 所做的相同。

据我所知，正在做的事情的唯一区别是 FFT 实现。有人知道这些结果是否是我对 FFT 算法总体理解的根本问题，或者是 Ron Mayer 编写的这段陈旧代码中基于 FHT 的实现的根本问题吗？您可以访问他的存档网站来获取我正在使用的代码: http://reocities.com/researchtriangle/8869/fft_summary.html

现在，如果我对一个数据 block 执行 FFT，然后执行 ifft，我会按照预期返回原始数据。如果我以任何方式修改光谱数据，我会得到意想不到的结果。

我曾经尝试用这种mayer FFT算法来代替S. Bernsee的音调变换算法中使用的算法，但这根本不起作用。我使用了 fftw3，该代码按我的预期工作。我很好奇用 fftw3 尝试同样的基本算法，看看会发生什么。

我不知道的是，我是否误解了一些导致我错误应用 rmayer 实现的基本原理，或者这只是一个我必须解决的简单工件(即，使用 FFT 大小的两倍)我预计)。

Answer 1

哦!这是诸如忘记在行尾添加分号之类的事情之一。卷积是频域中的复杂乘法——我正在做一个简单的逐点乘法。这是更正后的代码，显示了复数乘法。当然，有 C 和 C++ 结构和宏/例程来执行此操作，但这里是用于教学目的的强力方法: 假设 struct cxType 定义为:

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);
    }

现在上面的代码片段在被调用时可以工作，并且不会产生我所说的这些奇怪的镜像效果。至于我的 Matlab 代码，它之所以有效，是因为 Octave/Matlab 将复杂乘法的细节隐藏在它们方便的语法 H.*X 后面。

我能够通过分别乘法和组合实部和虚部来模仿我在 C 中的错误，从而在 Octave 中重现我的问题，然后通过执行与上面类似的操作来修复它。也就是说，rmayer 的 FFT 实现没有什么独特之处……只是我的卷积实现有问题。