2013-03-11 143 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