2013-03-11 145 views
1

给大家的好日子!Matlab:conv() - > fft()* fft() - > ifft()

我试图通过观察与一些已知脉冲响应的卷积来获得原始信号的基本问题。

但我得到的结果是某种程度上完全错误的,可能我在这里混合了不同的错误步骤。我已经在这里和其他网站上看过类似的话题,比如developpez,但是没有弄清楚原因。我将不胜感激任何帮助。

比方说,我的真实信号˚F []只是当时1的冲动,和脉冲响应 []为高斯。我计算他们的卷积h [。] conv(),然后,基本上,想要找到ifft(fft[h]./fft[g]),期待这是f [。]。

的第一个问题是conv()使N + M-1元素,其中N,M都是参数阵列的长度的阵列。所以,要执行fft[h]./fft[g]我需要做长度为g的smth。这是我可能犯错的第一个可疑的地方(见代码)。什么是正确的做法?

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

第三个问题是我无法理解如何采取信号转换。在matlab中,我必须使用正时信号进行操作,但是,例如,高斯脉冲响应既有时间负向的,也有时间正向的元素,因此,在这里使用它时,我需要将它向前移动(偷看者会向右移动),而且我需要“移动”结果?

谢谢!

这里是我的废话说:)

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! 
+0

FFT假定您的信号是周期性的。你可以通过使用适当的填充来处理这个问题。 [Numerical Recipes](http://www.nr.com)中的FFT章节对此进行了讨论,在此可能会有所帮助。 – sfstewman 2013-03-11 21:31:49

+0

另请参阅维基百科对[圆形卷积](http://en.wikipedia.org/wiki/Circular_convolution)的解释,特别是图表。使用FFT/IFFT来计算卷积是可能的,但您必须注意细节。 http://dsp.stackexchange.com上的这个问题你可能会有更好的运气。 – mtrw 2013-03-11 21:36:07

+0

相关问题:[验证卷积定理](http://stackoverflow.com/questions/14025967/verify-the-convolution-theorem) – 2013-03-11 22:10:02

回答

2

卷积的直接使用fft将导致circular convolution,而你想要什么(什么conv一样)是linear convolution。所以要用fft实现这样的方案,你必须将信号零点填充到长度为m+n-1

这里的表示convfft基于线性卷积的输出之间的等效示例:

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. 

可以检查out1out2是(FP精度内)是相同的。

以这种方式重新表达你的问题,你应该很好去。我无法弄清楚你在问什么:换班,但你可能想看看fftshiftifftshift

2
  1. 您可能需要使用filter(g, 1, f)而不是conv(g, f)以避免导致长的问题。或者剪切结果数组。
  2. 看起来你已经忘记了做一个FFT: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”算法来提高质量。它几乎相同,但是您需要通过特殊窗口函数(例如Blackman)来乘以输入。添加填充也适用,但它是计算量最大的方式。
+0

+1:鉴于问题是关于'conv',而不是'filter',我会建议保留'conv'并使用命令'fft(Conv,N)将FFT结果后缀填充0“和'fft(ImpulseResp,N)',其中'N =长度(TrueSignal)+长度(ImpulseResp) - 1'。此外,没有必要将fft(Conv)和fft(ImpulseResp)除以信号的长度。 – 2013-03-11 23:14:08

+0

非常感谢!我会试试看,我绝对需要一些关于FFT的理论读物。 – agronskiy 2013-03-12 07:49:05