作者：桂。
時間：2017-10-10 22:38:46
鏈接：http://www.cnblogs.com/xingshansi/p/7648274.html
前言
陣列信號處理中，經常用到小數延遲（fractional delay，FD）的思路，例如Beamforming、GSC等等，本文摘錄幾個小數延遲的實現方式，不打算做系統性的梳理，具體可參考課件。
一、問題模型
給出用到的資料：
1）部分code：網盤code.
2）stanford課件，對應鏈接：https://ccrma.stanford.edu/~jos/Interpolation/
3）Beamforming應用實例：http://www.labbookpages.co.uk/audio/beamforming/fractionalDelay.html
以均勻線陣為例：
設麥克風陣列共用M個陣元，中心為參考點，陣元間距為d，信號入射角為θ，聲音傳播速度為c，則根據幾何知識，第m(0≤m≤M-1)個陣元的時延為τm = (d/c) sinθ(m-(K-1)/2)。
麥克風采集的是數字信號，設采樣周期為T，則對時域離散的信號來說，時延為D = τ/T。通常D不是一個整數，而對離散信號來說，整數時延才有意義。對於非整數D，可以分解為整數部分和分數部分D = ⌊D⌋ + d，式中，⌊D⌋為D的向下取整，0≤d<1。對於非零的分數部分d，此時信號實際值介於兩個相鄰采樣點之間，即分數延遲。在實際處理中，可對d四舍五入取整，然后加上⌊D⌋，得到近似整數時延，但這種方法處理的結果不夠精確。為了得到精確的結果，通常借助小數延遲的思路。
二、小數延遲濾波器
A-一階FIR設計
思路主要來自Taloy一階近似：
從而
即：
其中η為對應的小數延遲，濾波器架構（低通信號有效）：
B-一階IIR設計
此時對應的濾波器為全通濾波器（All pass），近似逼近
濾波器響應：
濾波器架構：
對應時間延遲：
C-Sinc逼近
根據小數延遲濾波器的特性：
- 幅度響應：全通
- 相位響應：線性
得出濾波器：
對於采樣信號，需要限定在-fs/2 ~ fs/2之間，即相當於對濾波器進行了頻域截斷，截斷的濾波器特性：
求解該濾波器：
容易證明該濾波器是原型濾波器均方誤差最小的逼近。
%線陣為例 delay = d*sin(theta)*fs/c*(0:1:element_num); n = -64:63; i = 3;%延遲的陣元 h = sinc(n+delay(i)); H = zeros(1,256); H(1:128) = h; output = ifft((fft(sig)).*(fft(H)));%sig為輸入信號 %相比時域卷積，頻域點乘進一步節省資源
存疑：如果不是[-pi pi]，而是取[0 2*pi]，即sin(2*pi*f0*t)其中f0為大於fs/2的信號，簡單的逼近結果錯誤，延遲如何實現？
已解決，具體參考：過采樣的小數延遲實現。
D-Sinc加窗
由於C中濾波器存在截斷，防止能量泄露的一個常用思路就是：加窗截斷。
即：
α < 1 provides for a nonzero transition band。
對應code:
function h = hsincw2(L,d,wp,win)
% HSINCW2
% MATLAB m-file for sinc windowing method for FD filter design
% h = hsincw2(L,d,wp,win) designs an (L-1)th-order FIR
% filter to approximate a fractional delay of d samples,
% where 0 <= d < 1, wp is the passband of approximation and
% win is a length-L window function
% (e.g., win = chebwin(L,ripple), with sidelope ripple in dB).
% Output: length-L filter coefficient vector h
% Function Calls: standard MATLAB functions and sinc.m
%
% Timo Laakso 23.12.1992
% Revised by Vesa Valimaki 19.10.1995
% Last modified 14.01.1996
N = L-1; % filter order
M = N/2; % middle value
if (M-round(M))==0 % if L is even...
D = M + d; % D = M + d
else D = M + d -0.5; % ...otherwise
end;
b = (0:N)-D; % sample instants
h = sinc(wp*b); % shifted & sampled sinc function
h = h.*win; % windowing by the given window function
E-拉格朗日插值
該思路主要是借助拉格朗日插值，近似得出小數延遲位置的數值。關於拉格朗日插值法的介紹有很多。
濾波器逼近：
基於Taloy展開的性質：
這一特性符合拉格朗日插值的思路，利用此思路：
得出濾波器實現架構
對應code:
function h = lagrange(N, delay) %LAGRANGE h=lagrange(N,delay) returns order N FIR % filter h which implements given delay % (in samples). For best results, % delay should be near N/2 +/- 1. n = 0:N; h = ones(1,N+1); for k = 0:N index = find(n ~= k); h(index) = h(index) * (delay-k)./ (n(index)-k); end
不同插值個數對應的系數：
F-Farrow濾波器
該思路為多項式擬合，即將每一階h看作是多項式擬合：
濾波器可表述為：
簡化：
α可存在RAM里，利用查找表直接調用，便可以實現快速的小數延遲。
其中，參數求解：
對應實現架構：
