设计一个分数时延滤波器（matlab教程翻译+Python代码实现）

设计一个分数时延滤波器（matlab教程翻译+Python代码实现）

原文：Design Fractional Delay FIR Filters

整数时延

假设数字信号的时延过程为$y[n]=x[n-D]$，其中$D$是一个整数。可以将该操作转化为一个卷积滤波器，$y=h\ast x$，其中$h$是一个有限冲激响应$h[n]=\delta (n-D)$，相应的传递函数表示为$H(z)=z^{-D}$，其频率响应则表达为$H(\omega )=e^{-i\omega D}$。
python版本代码：

D = 3
h = [0 for i in range(D),1]   ## [0,0,0,1]

x = np.array([i for i in range(1,20,1)])
print(x)
y = np.convolve(x.h)
print(y)

通过将序列经由有限脉冲响应（FIR）滤波器h进行滤波处理以实现序列的移位。输出结果中的前导零点表征了该类滤波器固有的初始状态。

通过D/A插值实现非整数延迟处理

对于序列$y[n]=x[n-D]$，当延迟$D$为非整数值时，其延迟效果并未明确定义。为了合理地实现这类分数延迟并确保输出在连续域上被采样，必须在处理过程中引入一个D/A插值环节。具体来说，有$y[n]=\widehat{x_n}(n-D)$，其中$\widehat{x_n}$代表对输入序列$x$进行D/A插值的结果。$\widehat{x_n}$可能与$n$有关，并能表征从中获得序列$x$的底层模拟信号模型。该策略同样适用于其他重采样问题，例如采样率转换。
本示例采用了数字信号处理系统工具箱中的两种插值模型来实现分数延迟滤波器。
基于$Sinc$函数的插值模型：采用带限重建方法来得到$\widehat{x_n}$。
基于$Lagrange$的插值模型：利用多项式重建方法来计算$\widehat{x_n}$。

有限带宽分数延迟滤波器研究

根据香农-惠特克插值公式，$\hat{x}(t)={\sum }_{k}x[k]sinc(t-k)$ 能够对限带信号进行精确描述。该公式表明，数字模拟转换过程 $\hat{x}$ 是对输入序列的限带重建。针对特定的延迟值 $D$，我们可以通过 $y[n]=\hat{x}(n-D)$ 来表达分数延迟，其中对于任意 $n$，都采用相同的 $\hat{x}$ 作为卷积滤波器的系数。这种滤波器被称作理想限带分数延迟滤波器，其冲激响应定义为 $h_D[k]=\text{sinc}(k-D)$。相应的频率响应（即离散时间傅里叶变换）可以表示为 $H_D(\omega )=e^{-i\omega D}$。

因果FIR近似的理想带限移位滤波器

前文所描述的理想sinc-shift滤波器属于全通滤波器（即$|H_d(\omega )|=1$），然而其脉冲响应$h_D$是无限且非因果的，无法在计算中实现。

在实际应用及计算处理中，通常需要将理想滤波器在一个有限索引窗口内进行截断，尽管这可能导致一定程度的带宽损耗。设目标延迟值为$D$，期望的滤波器长度为$N$，则满足$|D-k|\le N/2$ 的索引$k$所构成的窗口将以$D$为中心对称，并能够覆盖理想滤波器的主要通带。当$D=i_0+FD$（其中$0\le FD\le 1$，且$i_0$为整数）时，具体的窗口索引集合为 { i 0 − ⌊ ( N − 1 ) / 2 ⌋ , … , i 0 + ⌊ N / 2 ⌋ } \{i_0-\lfloor {(N-1)/2} \rfloor, \ldots, i_0+\lfloor N/2 \rfloor\} {i0​−⌊(N−1)/2⌋,…,i0​+⌊N/2⌋}。在此，整数$i_0$被称作整数延迟，其值可任意选定。为确保有限脉冲响应（FIR）滤波器的因果性，可将$i_0$设定为$\lfloor N-1/2\rfloor$，从而得到索引窗口 { 0 , … , N − 1 } \{0, \ldots, N-1\} {0,…,N−1}。（译者理解：如果截取【-2：3】,则-2、-1位置需要使用滤波卷积后数据，即滤波器的输出被送入输入，非因果索引。而截取【0：5】,则都可以使用原始数据，始终是因果位置）

截断sinc滤波器会导致频率响应出现波动，为了解决这个问题，可以对FIR系数施加权重$w_k$（例如$Kaiser$或$Hamming$）。理想带限分数延迟滤波器的FIR近似模型可表示为以下形式：
$$h[k]=w_k{h}_d[k]=\{\begin{array}{ll}{w}_k\cdot \text{sinc}\left(k-FD-i_0\right)& 0\le k\le N-1\\ 0& \text{otherwise}\end{array}$$

python版本代码：

def sinc_base_FD_FIR(x,N,delay):
    FD = delay % 1
    int_delay = int(delay // 1)
    
    ## 分数延时滤波
    if 0 == FD:
        i0 = 0
        delayed_signal = x
    else:
        i0 = int(np.ceil(N/2) - 1)
        win_index_list = np.arange(0,N) - i0
        win_fun = np.hamming(N)
        sinc_filter = np.sinc(win_index_list - FD)
        sinc_filter = sinc_filter * win_fun 
        delayed_signal = np.convolve(x, sinc_filter, mode='FULL')    # 使用卷积来应用sinc滤波器

    
    ## 整数延迟
    t_delay = int_delay - i0
    res = np.zeros(len(delayed_signal) + 2 * np.abs(t_delay))
    if 0 == t_delay:
        res = delayed_signal
    elif 0 < t_delay:
        res[t_delay:t_delay + len(delayed_signal)] = delayed_signal
    else:
        t_delay = np.abs(t_delay)
        res[:len(delayed_signal) - t_delay] = delayed_signal[t_delay:]
        
    return res[:len(x)]

拉格朗日差值的分数时延滤波算法（最大平坦度准则）

分数时延滤波器的另一种设计方法是利用最大平坦度准则在频域去逼近理想的分数时延滤波器。

理想分数时延滤波器的频率响应可以表示为幅度响应和相位响应两个部分。对于理想的分数时延滤波器，其幅度响应在所有频率上都是常数1，即无失真地通过所有频率分量；而相位响应则与频率成线性关系，表示了不同频率成分经历的延迟。
幅度响应$A(f)$： $A(f)=1$
相位响应$\phi (f)$，其中$D$是理想的分数时延（以样本为单位）： $\phi (f)=2\pi fD$
综合起来，理想分数时延滤波器的频率响应H(f) 可以写作复数形式： $H(e^{j\omega })=e^{-j\omega D}$

假设设计的滤波器具有如下形式：$H(e^{j\omega })={\sum }_N^{n=0}h(n)e^{-j\omega n}$。
定义误差函数：$E(e^{j\omega })=H(e^{j\omega })-H(e^{j\omega })$
对误差函数求N阶导数，并使其导数在特定频率处为0，则${\frac{d^{N}E(x)}{d{x}^N}|}_{x=1}=0$
当$\omega =0$时，上式可以表示为：${\sum }_{n=0}^N{n}^{k}h(n)=D^k,k=0,1,2,...,N$
解得$h(n)={\Pi }_{k=0,k\ne n}^N\frac{D-k}{n-k},n=0,1,2,...,N$
将时延$D$以及阶数$N$代入公式，即可得到FIR分数时延滤波器的系数。最大平坦准则逼近法同拉格朗日差值法相同。

def lagrange_interpolation(x, F, j):
    '''拉格朗插值基函数
    x：x轴信号列表
    F: 小数时延
    j：
    '''
    result = 1.0
    for m in range(len(x)):
        if m != j and x[m] != x[j]:
            result *= (F - x[m]) / (x[j] - x[m])
    return result

def fractional_delay_filter(delay, L):
    N = int(np.floor(delay))
    F = delay - N
    k = int((L - 1) / 2)

    h = np.zeros(L)
    for i in range(L):
        h[i] = lagrange_interpolation(np.arange(L) - k, F, i)

    return h

def lagrange_inter_base_FD_FIR(x,N,delay):
    FD = delay % 1
    int_delay = int(delay // 1)
    
    ## 分数延时滤波
    if 0 == FD:
        i0 = 0
        delayed_signal = x
    else:
        i0 = int(np.ceil(N/2) - 1)
        win_index_list = np.arange(0,N) - i0
        filter_coeffs = fractional_delay_filter(FD,N)
        delayed_signal = np.convolve(x, filter_coeffs, mode='FULL')
    
    ## 整数延迟
    t_delay = int_delay - i0
    res = np.zeros(len(delayed_signal) + 2 * np.abs(t_delay))
    if 0 == t_delay:
        res = delayed_signal
    elif 0 < t_delay:
        res[t_delay:t_delay + len(delayed_signal)] = delayed_signal
    else:
        t_delay = np.abs(t_delay)
        res[:len(delayed_signal) - t_delay] = delayed_signal[t_delay:]
        
    return res[:len(x)]

基于Farrow结构的分数时延

假设滤波器的频响曲线为$H(e^{j\omega })={\sum }_{n=0}^{N-1}{C}_n(D)e^{-j\omega n}$
其中$D$为时延参数。由于分数时延滤波器系数会随着时延参数D的变化而变化，因此可以写为D的函数$C_n(D)$，$C_n(D)$可以用D的多项式来逼近。
$$C_n(D)=C_{n,0}{D}^0+C_{n,1}{D}^1+\cdots +C_{n,M-1}{D}^{M-1}=\sum _{m=0}^{M-1}{C}_{n,m}{D}^m$$
设计系数$C_{n,m}$，使其逼近理想的分数时延滤波器，满足下式：
$$min{\int }_{{\omega }_0}^{{\omega }_1}{\int }_{D_0}^{D_1}|\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}{C}_{n,m}{D}^m{e}^{-j\omega n}-e^{-j\omega D}{|}^{2}dDd\omega$$
其中，$\omega \in [{\omega }_0{\omega }_1]$表示数字频率上线性时延带宽，$D\in [D_0{D}_1]$表示分数时延参数的范围。
Farrow结构有M组N阶FIR滤波器组成。因此，Farrow结构所需要的滤波器阶数远大于前两种滤波器结构，并且运算量也大。但是，他的优点是当时延变化时，仅改变时延参数D，就可以获得不同的分数时延，不需要重新加载系数，节省存储空间，降低硬件实现的复杂度。