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

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

我曾经尝试用这个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中的错误，得以在八度中再现我的问题，然后通过做类似于上面的事情来修复它。这就是说，Rmayer对FFT的实现没有什么独特之处...只是我对卷积的实现有问题。