作者：桂。
時間：2017-03-07 22:33:37
鏈接：http://www.cnblogs.com/xingshansi/p/6517301.html
前言
信號時域、頻域對應關系，及其DFT、FFT等變換內容，在之前的文章1、文章2中已經給出相關的理論推導以及代碼實現，本文主要針對信號中常用的卷積進行介紹，內容主要包括：
1）卷積的物理意義；
2）卷積的直接實現；
3）卷積的FFT實現；
4）卷積的無延遲快速實現；
本文為自己的學習總結，內容多有參考他人，相關的參考文獻名稱最后一並給出。
一、卷積的物理意義
知乎上有一段話，非常經典，此處加以引用：
作者：中微子
鏈接：https://www.zhihu.com/question/20500497/answer/45708002
打個比方，往平靜的水面里面扔石頭。我們把水面的反應看作是一種沖擊響應。水面在t=0時刻石頭丟進去的時候會激起高度為h(0)的波紋，但水面不會立馬歸於平靜，隨着時間的流逝，波紋幅度會越來越小，在t=1時刻，幅度衰減為h(1), 在t=2時刻，幅度衰減為h(2)……直到一段時間后，水面重復歸於平靜。
從時間軸上來看，我們只在t=0時刻丟了一塊石頭，其它時刻並沒有做任何事，但在t=1,2….時刻，水面是不平靜的，這是因為過去（t=0時刻）的作用一直持續到了現在。
那么，問題來了：
如果我們在t=1時刻也丟入一塊石子呢？此時t=0時刻的影響還沒有消失（水面還沒有恢復平靜）新的石子又丟進來了，那么現在激起的波浪有多高呢？答案是當前激起的波浪與t=0時刻殘余的影響的疊加。那么t=0時刻對t=1時刻的殘余影響有多大呢？
為了便於說明，接下來我們作一下兩個假設：
1． 水面對於“單位石塊”的響應是固定的
2． 丟一個兩倍於的“單位石塊”的石塊激起的波紋高度是丟一個石塊的兩倍（即系統滿足線性疊加原理）
現在我們來計算每一時刻的波浪有多高：
- t=0時刻：
y(0)=x(0)*h(0);
- t=1時刻：
當前石塊引起的影響x(1)*h(0);
t=0時刻石塊x(0)引起的殘余影響x(0)*h(1);
y(1)=x(1)*h(0)+ x(0)*h(1);
- t=2時刻：
當前石塊引起的影響x(2)*h(0);
t=0時刻石塊x(0)引起的殘余影響x(0)*h(2);
t=1時刻石塊x(1)引起的殘余影響x(1)*h(1);
y(2)=x(2)*h(0)+ x(1)*h(1)+x(0)*h(2);
……
- t=N時刻：
當前石塊引起的影響x(N)*h(0);
t=0時刻石塊x(0)引起的殘余影響x(0)*h(N);
t=1時刻石塊x(1)引起的殘余影響x(1)*h(N-1);
y(N)=x(N)*h(0)+ x(N-1)*h(1)+x(N-2)*h(2)+…+x(0)*h(N);
這就是離散卷積的公式了。
二、卷積的直接實現（Direct convolution）
特點：無延遲、運算量大。
對應公式：
$y(n) = h(0)x(n) + h(1)x(n - 1) + ... + h({N_h} - 1)x(n - {N_h} + 1)$
對應結構圖：
對應代碼(MATLAB同樣有內置指令conv)：
function output_signal=my_direct_convolution(input_signal,impulse_response)
% Input:
% input_signal: the input signal
% impulse_response: the impulse response
% Output:
% output_signal:the convolution result
N=length(input_signal);%define length of signal
K=length(impulse_response);%define length of impulse response
output_signal=zeros(N+K-1,1);%initializing the output vector
xp=[zeros(K-1,1);input_signal;zeros(K-1,1)];
for i=1:N+K-1
output_signal(i)=xp(i+K-1:-1:i)'*impulse_response;
end
三、卷積的FFT實現
特點：運算量小、需要延遲至少$N_h+N_x$（$N_h$為濾波器長度）個數據點的采樣時間。
A-卷積與FFT的關系
這里用到一個重要性質：時域卷積等價於各自頻域變換相乘。但我們知道，fft或者ifft信號的長度不變，即fft實現的頻域相乘，在逆fft（ifft）到時域的結果，實現的是圓卷積，而不是線卷積。
$y(n) = x(n)*h(n)$
對應於頻域就是
$Y(k) = X(k) \cdot H(k)$
但$Y(k)、X(k)、H(k)$對應的長度必須一致，且與$y(n)$一致。
為了保持線卷積、圓卷積得到的結果一直，必須在FFT操作之前，適當補零。
舉個簡單的例子：
x(n) = [2,5,7]; h(n) = [3,4]，求二者的線卷積.
分析：線卷積的結果長度應該為：$N_y = N_x + N_h -1$，即3+2-1.
方法1：直接卷積
y(n) = x(n)*h(n),對應結果（會乘法運算的都能做卷積）：
方法2：圓卷積
因為利用FFT操作，對應圓卷積，且x、h長度需要一致，這里我們先取長度為3，看看效果，長度不足3的補零。
從結果看，結果為{34,23,41}，與結果不一致，可以看到28對應的位置有數值，因此至少取$N_x +N_h$個長度，再來看看結果：
此時已經得出了真實的結果，結論就是：做FFT實現卷積，每一個序列長度補零延長到$N_x +N_h$，則圓卷積與線卷積等價，即FFT運算便可以實現卷積的功能。
B-卷積的FFT代碼實現
給出對應的代碼：
function output_signal=my_fft_convolution(input_signal,impulse_response) % Input: % input_signal: the input signal % impulse_response: the impulse response % Output: % output_signal:the convolution result siglen=length(input_signal);%define the length of signal implen=length(impulse_response);%define the length of the impulse response P=siglen+implen-1;%define P P2=pow2(nextpow2(P));% Find smallest power of 2 that is > P %%Utilize the code in .pdf output_signal= real(ifft(fft(input_signal,P2).*fft(impulse_response,P2))); % we need to take the real part to avoid any % Floating point errors making the output complex output_signal(P+1:end)=[];%modify the length
給一張直接卷積與FFT卷積的框圖：
四、卷積的分段實現
特點：相比於直接卷積，運算量小；相比於FFT實現，延遲量小（$2N_h$，假設數據長度大於濾波器長度）。
A-重疊相加法
主要思想：
假設每一段長度$N_x$，濾波器長度$N_h$，則每一段線卷積的長度為$N_x+N_h-1$，而每一段長度為$N_x$，故每次移動$N_x$，從而每次的結果后段的$N_h-1$個點重疊相加。核心是利用線卷積的思想，如前文所示，線卷積也可借助FFT快速實現。
對應代碼：
function output_signal=my_fast_convolution_add(input_signal,impulse_response)
% Input:
% input_signal: the input signal
% impulse_response: the impulse response
% Output:
% output_signal:the convolution result
siglen=length(input_signal);%define the length of signal
implen=length(impulse_response);%define the length of the impulse response
Nfra=ceil(siglen/implen);%number of frames
%Initialize the output signal
output_signal=zeros(implen*(Nfra+1),1);
%Add zeros
xp=reshape([input_signal;zeros(Nfra*implen-siglen,1)],...
implen,Nfra);
h=kron(impulse_response,ones(1,Nfra));
%Multiply the ffts together and calculate the ifft
P=pow2(nextpow2(2*implen));%Find smallest power of 2 that is > 2*implen
out_mat=real(ifft(fft(xp,P).*fft(h,P)));
%Overlapping the output blocks and summing
for i=0:Nfra-1
ni=i*implen+1:i*implen+P;
output_signal(ni)=output_signal(ni)+out_mat(:,i+1);
end
%Modify the length
output_signal(siglen+implen:end)=[];
B-重疊保留法
主要思想：
濾波器長度$N_h$，每一段卷積的長度為$2N_h-1$，核心是利用圓卷積的思想，對應頻域就是FFT相乘。
舉例：
已知$x(n) = n+1,0 \le n \ge 9 $，$h(n) =$ {1,0,-1}，求卷積結果。
直接實現：
$x(n)*h(n) =$ {1 2 2 2 2 2 2 2 2 2 -9 -10} .
重疊相加：
已知$N_h = 3$，長度L最小值為$2N_h-1 = 5$.將$x(n)$分段，得
$x_1(n) =$ {1 2 3 4 5}，$x_2(n) =$ {6 7 8 9 10}，
將每段與$h(n)$卷積：
$y_1 =$ {1 2 2 2 2 -4 -5 },
$y_2 =$ {6 7 2 2 2 -9 -10 }，
重疊$N_h-1 = 2$個點，的
$y =$ {1 2 2 2 2 2 2 2 2 2 -9 -10}
重疊保留：
同樣的計算方法，長度最小為5，起始保留$N_h -1 =2$個點，起始補零。
$x_1 =$ {0 0 1 2 3};$x_2 =$ {2 3 4 5 6};
$x_3 =$ {5 6 7 8 9};$x_4 =$ {8 9 10 0 0};
分別與$h(n)$做圓周卷積：
$y_1 =$ {-2 -3 1 2 2};$y_2 =$ {-3 -3 2 2 2};
$y_1 =$ {-3 -3 2 2 2};$y_1 =$ {8 9 2 -9 -10};
丟棄每段最前面的$N_h -1 = 2$個點，為什么這么丟棄，參考前文線卷積與圓卷積的對應關系分析。
$y =$ { 1 2 2 2 2 2 2 2 2 2 -9 -10}.
三種方式完全等價。
這一方法也可以解決輸入、輸出數據量不匹配的問題：
假設x長度為16，h長度為16，則起始保留Nh-1 = 15個點，起始補零。圓卷積借助FFT實現。重疊相加法也可以實現該功能。
五、卷積的Without Input-output delay實現
特點：無延遲、運算量小。
該方法主要針對濾波器進行分層。
具體可以參考文章：《Efficient convolution without Input-Output Delay》.
下面給出一個自己制作的gif示意圖，之前的圖片沒有了，pdf截圖有大有小，拼湊起來湊合看吧，這個圖描寫了該算法的步驟。
參考：
鄭君里：《信號與系統》第二版