Matlab: conv() -> fft() -> fft()

Question

祝大家今天愉快！

我试图解决的基本问题，从观察原始信号的卷积与一些已知的脉冲响应。

但是我得到的结果是完全错误的，而且我可能有不同的错误步骤。我已经在这里和其他网站上看过了类似的话题，比如developpez，但是没能找出原因。我很感谢你的帮助。

假设我真正的信号f是时间1的一个脉冲，脉冲响应g是高斯的。我用conv()计算它们的卷积h，然后，基本上，想要找到ifft( fft[h]./fft[g] )，期望它是f。

第一个问题是，conv()生成n+m-1元素数组，其中n，m是参数数组的长度。因此，要执行fft[h]./fft[g]，我需要执行长度为g的smth操作。这是我可能会出错的第一个可疑的地方(参见代码)。正确的方法是什么？

第二个问题是，我得到的信息与最初的真实信号非常不同。

第三个问题是，我不知道如何接收信号移位。在matlab中，我必须对正时信号进行操作，但是，例如，高斯脉冲响应有时间-负和时间-正两种元素，所以，为了在这里使用它，我需要将它“向前”( peek将移到右边)，而不是我需要“去移动”结果？

谢谢!

这是我的废话:)

close all;

TrueSignal = zeros( 101, 1 ); % impulse in t = 1.
TrueSignal( 1 ) = 1;
ImpulseResp = normpdf(-1:0.02:1)/normpdf( 0 ); % 101 elements array

figure;
subplot( 2,2,1 );
title('True signal')
plot( TrueSignal );
subplot( 2,2,2 );
title('Impulse response')
plot( ImpulseResp );

Conv = conv( TrueSignal, ImpulseResp ); % produces 201 elements array.
subplot( 2,2,3 );
title('Convolution')
plot( Conv );

% Wrong? I need a 201 elements array to represent the impulse response.
ImpulseResp_sparse = normpdf( -1:0.01:1 )/normpdf( 0 );
FIR = fft( ImpulseResp_sparse )/201;

Inverse = ifft( fft( Conv )./FIR ); % UPD Added fft() according to one of comments, bad mistake, but still not preventing.

subplot( 2,2,4 );
title('What is that???')
plot( abs( Inverse ) ); % It's weird! With no abs(), result is even more weird! 

Answer 1

直接使用fft进行卷积将导致圆卷积，而您想要的(以及conv所做的)是线性卷积。因此，要用fft实现这样的方案，您必须将信号压缩到长度为m+n-1。

下面是一个示例，显示基于conv和基于fft的线性卷积输出之间的等价性：

x=rand(4,1);y=rand(3,1); %sample data
out1=conv(x,y);          %output from conv()
X=fft(x,6);Y=fft(y,6);   %zero pad and compute fft
out2=ifft(X.*Y);         %output from fft based lin. conv.

您可以检查out1和out2是否相同(在FP精度范围内)。

用这种方式重新规划你的问题，你应该做得很好。我不知道您要求的是什么:轮班，但您可能需要查看fftshift和ifftshift。

Answer 2

1. 您可能希望使用filter(g, 1, f)而不是conv(g, f)来避免由此产生的长度问题。或者剪掉结果数组。
2. 看起来你好像忘了做一个快速傅立叶变换：Inverse = ifft( fft(Conv)./FIR );

这对我来说很管用：

   close all;

   TrueSignal = zeros( 101, 1 ); % impulse in t = 1.
   TrueSignal( 1 ) = 1;
   ImpulseResp = normpdf(-1:0.02:1)/normpdf( 0 ); % 101 elements array

   figure;
   subplot( 2,2,1 );
   title('True signal')
   plot( TrueSignal );
   subplot( 2,2,2 );
   title('Impulse response')
   plot( ImpulseResp );

   Conv = filter( TrueSignal, 1, ImpulseResp ); 
   subplot( 2,2,3 );
   title('Convolution')
   plot( Conv );

   fftConv = fft(Conv);

   FIR = fft( ImpulseResp );

   Inverse = ifft( fftConv./FIR );

   subplot( 2,2,4 );
   plot( abs( Inverse ) );

1. 正如所说，FFT假定信号是周期性的。如果要分析非周期信号，需要使用“加窗FFT”算法来提高质量。这几乎是一样的，但是您需要用一个特殊的窗口函数(例如布莱克曼)来乘以输入。添加填充也有效，但它是计算成本最高的方法。