信息处理 - 使用 FFTW 在频域中执行卷积？ - 吾爱随笔录

使用 FFTW 在频域中执行卷积？

信息处理 fft 卷积频域 C

2022-02-01 02:36:13

我试图通过使用FFTW库的函数对C 中的两个信号和进行卷积，对每个信号执行傅里叶变换，将适当的复杂分量相乘，然后对结果乘积进行 IFFT。不幸的是，我的 IFFT 结果不正确，我不确定我做错了什么。 $x(n)$ $h(n)$

我正在使用我在网上找到的代码的修改版本（此处为 .c 和 Makefile ）进行测试，该版本使用 0 到 1 之间的随机数填充两个数组。我的过程如下：

分别为长度和长度的（音频信号）和（脉冲响应）创建double缓冲区。xhnnIR
用 0 到 1 之间的随机数填充每个缓冲区。
创建两个complex缓冲区X和H，每个缓冲区的长度nOut = n + nIR - 1
创建并执行两个从实数到复数的 FFT 计划，每个信号一个，都使用nOut元素（即plan_forward = fftw_plan_dft_r2c_1d ( nOut, x, X, FFTW_ESTIMATE );和plan_forwardIR = fftw_plan_dft_r2c_1d (nOut, h, H, FFTW_ESTIMATE);
创建复数数组fftMulti以保存两个数组的频率乘积。它的大小是nc = nOut / 2 + 1因为 FFTW 不存储 FFT 的冗余一半
循环for ( i = 0; i < nc; i++ )并执行复数乘法，即fftMulti[i][0] = X[i][0] * H[i][0] - X[i][1] * H[i][1];实部和fftMulti[i][1] = X[i][0] * H[i][1] + X[i][1] * H[i][0];虚部
创建大小double数组。该数组将保存 IFFTconvolvedSignOutfftMulti
创建并执行大小为的复数到实数的 FFT 计划nOut，即plan_backwardConv = fftw_plan_dft_c2r_1d(nOut, fftMulti, convolvedSig, FFTW_ESTIMATE);
convolvedSig通过将每个元素除以标准化nOut

中的值convolvedSig不正确（它们太高了），但我不确定我做错了什么。我还为我未修改的信号X和H（也是大小nOut）创建了数组和复数到实数的计划，并且这些 IFFT 工作得很好（即，值与我对它们执行 FFT 之前的值完全相同）。

使用我的流程和/或我的代码，有人可以帮我确定我做错了什么吗？

2个回答

您需要将 x 和 h 归零。

nZeros = 长度(x) + 长度(h) -1

MATLAB 中的示例：

clear all; close all; home

xSamples = 10;
hSamples = 7;

x = randn(1,xSamples);

h = randn(1,hSamples);

nZeros = xSamples + hSamples - 1;

X = fft(x,nZeros);
H = fft(h,nZeros);

x_conv_h = ifft(X.*H);


figure(1)
plot(real(x_conv_h))
hold on; grid on
plot(conv(x,h),'.r')
legend('fft based convolution','convolution')

真实信号的 FFT 卷积非常容易。频域中的乘法相当于时域中的卷积。但是要记住的重要一点是，您正在乘以复数，因此，您必须进行“复数乘法”。

我将在下面为此提供一些 C 源代码。

所以步骤是：

对过滤器内核进行 FFT，
对“干”信号进行 FFT。
做两个光谱的复数乘法
执行这个新频谱的逆 FFT。

当然，如果要对长信号进行连续处理，则需要使用重叠添加或重叠保存方法。

如果您只使用真实信号，在 Intel 格式（小端序）机器上，您可以使用 Surreall FFT 以及我在下面给出的乘法函数，该函数采用正确的数据顺序格式。

这是代码：

#define F_TYPE float;

void CompMulR(F_TYPE a[],F_TYPE b[],F_TYPE res[],long n)
  {
  long i;
  for(i=0;i<n;i++)
    {
    if(i<2)
      res[i]=a[i]*b[i]; // DC & NQ are not complex
    else
      {
      if((i&1)==0)
        {//real
        res[i]= a[i]*b[i] - a[i+1]*b[i+1];
        }
      else
        {//img.
        res[i]= a[i]*b[i-1] + a[i-1]*b[i];
        }
      }
    }
  }



#define N 1024
F_TYPE TimeDom [N];
F_TYPE FrqDom [N];
F_TYPE FrqDom2 [N];
F_TYPE Twidds[N/4];

InitFFT(N, Twidds); //Initialse the twiddle factors
FFT(TimeDom,FrqDom,N,Twidds); // do a forward transform
CompMulR(TimeDom,FrqDom,FrqDom2[],N);//do a complex multiply 
IFFT(TimeDom,FrqDom2,N,Twidds); // do an inverse transform 
//TimeDom now holds the convolution of the 2 signals

FFT 的代码在这里： http ://ravellescientific.co.uk/sureal-fft/ 希望这会有所帮助

其它你可能感兴趣的问题

上一篇结构化图像和纹理图像有什么区别下一篇如何在没有 DSP 工具箱的情况下设计窄带通滤波器