数字信号处理之卷积 (3)

在学习数字信号处理的时候,卷积是个让人又恨又爱的东西,一方面很难理解它的真实效果,初次接触非常难以建模,另一方面在信号处理中又十分常见。本文将介绍一下卷积的特点,以及离散系统的线性卷积和循环卷积运算,如何利用 FFT 计算卷积

1. 卷积 (convolution) 介绍

在之前我们初步介绍过卷积的概念,同时在线性时不变系统中介绍卷积在这种系统分析的巨大作用,这里先初步介绍下卷积的性质。

注: 数学上我们通常用 * 表示卷积,即连续系统的 \[
y(t) = x(t) \ast h(t) = \int_{-\infty} ^{\infty} {x(\tau) h(t-\tau) \mathrm d \tau}
\]
或者离散系统的 \[
y(n) = x(n) \ast h(n) = \sum_{m= – \infty}^{\infty} x(m) h(n-m)
\]
卷积运算的性质:

  • 交换律,\(f \ast g = g \ast f\)

  • 结合律,\(f \ast (g \ast h) = (f \ast g)\ast h\)

  • 分配律,\(f \ast (g+h) = (f \ast g) + (f \ast h)\)

  • 微积分性质 \[
    \cfrac {d^n} {d t^n} [f_1(t) \ast f_2(t)] = \cfrac {d^n f_1(t)} {d t^n} \ast f_2(t) = f_1(t) \ast \cfrac {d^n f_2(t)} {d t^n}
    \]

    \[
    \int_{-\infty} ^t [f_1(\tau) \ast f_2(\tau)] d \tau = [\int_{-\infty} ^t f_1(\tau) d \tau ] \ast f_2(\tau) = f_1(\tau) \ast [\int_{-\infty} ^t f_2(\tau) d \tau ]
    \]

最重要的性质:

时域的卷积运算对应频域上的乘积, 即 \[
\mathcal{F} (f \ast g) = \mathcal{F}(f) \cdot \mathcal{F}(g)
\]
这一定理对拉普拉斯变换和 \(z\) 变换同样适用

2. 线性卷积如何计算 ?

在离散系统中,我们将卷积计算也称为线性卷积,即 \[
y(n) = \sum_{m= – \infty}^{\infty} h(m) x(n-m)
\]
性质:

如果 \(h(m)\)\(x(n-m)\) 只在有限长范围取值,长度分别为 \(N\)\(M\), 那么线性卷积结果的长度一定为 \(N+M-1\)

在看懂移位加权法后,可以很清楚看到这一点

我们用列表法移位加权法计算线性卷积

2.1 列表法 – 逐点 (pointwise) 计算

列表法是\(m\) 作为变量,\(n\) 作为参数,预先设置具体的 \(n\) 值,可以获得这个 \(n\) 值对应的 \(y(n)\) ,形象来说,列表法可以直接获得单个 \(n\) 对应的卷积值,但如果想要获得所有 \(n\) 对应的值,那就要计算所有的 \(n\) ,这种方法的优势在于可以很清楚的看到动态改变 \(n\) ,每次都可以计算出卷积值

步骤如下:

  • 首先将 \(x(n)\)\(h(n)\) 的自变量从 \(n\) 变为 \(m\), 其他都不变
  • 然后将 \(h(m)\) 翻转,得到 \(h(-m)\)
  • \(n\) 是为常量参数,则 \(h(n-m)\) 是将 \(h(m)\) 右移 \(n\) 个单位得到的

如图所示,我们首先将 \(f(t)\)\(g(t)\) 的自变量变为 \(f(\tau)\)\(g(\tau)\), 然后对 \(g(\tau)\) 进行反褶,得到 \(g(-\tau)\), 此时平移 \(t\), 就可以得到随之移动的函数 \(g(t-\tau)\) ,这里要注意的是,增加 \(t\), 函数是向右移动的,不是向左移动

2.2 移动加权法

移位加权法是\(n\) 作为变量,\(m\) 作为参数 , 此时卷积就是信号 \(h(n)\)\(x(n)\) 的每个时延版本 (时延参数为 \(m\) ) 的加权叠加,这种方法的特点非常符合人的理解方式。

比如我们生病了,在一定时间间隔内打了5次针,那么每个时刻我们体内的残留药效是每次打针效果的时延叠加 (时延是指打针的时间距离现在时间的时间间隔,显然间隔越长,药效越小)。

再比如说,FIR 滤波器,我们是如何计算输出的呢 ? 我们会把输入的多个时延和我们设计的滤波器系数相乘求和得到的 (数学上叫“加权平均”)。所以设计不同的滤波器的效果就是好比控制每次打针的效果权重。一种最常见的FIR滤波就是均值滤波。

这种方法的缺点是,计算卷积值是一次性计算出来的,所有点的值是所有时延最后一次性叠加出来的,我们不能只算出一个点的卷积。这种方式在实际计算时,如果是无穷个时延叠加,明显是不现实的。

举个离散系统的例子来理解, \[
x(n) = \{1, 1, 1\}, 0 \le n \le 2; \;\; h(n) = \{4, 3, 2, 1\}, 0 \le n \le 3
\]
求解线性卷积: \(y(n) = x(n) \ast h(n)​\)

2.3 总结

总结: 第一种方法适合单点卷积值的计算,第二种方法适合理解, 这也是卷积计算的本质。

同时在离散系统也适合计算线性卷积值 (注意是手动计算,不是机器计算)。

可以明显感受到,如果我们采用计算机计算卷积和,只能用第一种方法,因为第二种方式需要用相当大的内存来存储临时变量 (每次移位并乘以系数),这明显不现实。同时,我们在计算实时卷积的时候,比如实时滤波,我们也是实时更新最新点的滤波值。

因此,推荐第一种方法计算,第二种方法理解。

下面给出 C++ 版本的逐点卷积计算,样本数组和核进行卷积。

注意:

  1. 这里没有考虑 数组 x[i] 左右越过边界的问题,在实际计算可以将超过边界的点的值设为 0,或者跳过计算
  2. 卷积结果可能溢出。比如这里 x[i] 是 8-bit 的无符号整型,加权求和的计算结果非常有可能会溢出,比如小于 0 或者大于 255,我们必须检查计算结果
for ( i = 0; i < sampleCount; i++ )
{
    y[i] = 0;                      // set to zero before sum
    for ( j = 0; j < kernelCount; j++ )
    {
        y[i] += x[i - j] * h[j];   // convolve: multiply and accumulate
    }
}

3. 循环卷积如何计算 ?

注意循环卷积和线性卷积不同之处在于,线性卷积是无穷个时延的叠加,而循环卷积要指定卷积的点数,\(L\) 点的循环卷积如下: \[
y_L(n) = \sum_{m=0}^{L-1} x(m) (h(n-m))_{L} R_{L}(n) = \sum_{m=0}^{L-1} h(m) (x(n-m))_{L} R_{L}(n)
\]

还是相同的例子,

\[
x(n) = \{1, 1, 1\}, 0 \le n \le 2; \;\; h(n) = \{4, 3, 2, 1\}, 0 \le n \le 3
\]
求解 4 点和 8 点的循环卷积: \(y(n) = x(n) \ast h(n)\)

3.1 矩阵方程法

\(m\) 为变量

  1. 首先将 \(x(n)\)\(h(n)\) 补零到 \(L\) 点长;
  2. 将其中一个序列周期延拓、翻褶、取主值、循环右移构成方阵;
  3. 将另一个序列写成列矩阵
  4. 二者做矩阵乘法运算

4 点和 8 点的循环卷积如下所示:

3.2 循环移位加权和法

\(n\) 为变量

  1. \(h(n)\) 补零到 \(L\) 点长;
  2. \(h(n)\) 循环右移 \(m\) 个单位,并乘以加权系数 \(x(m)\)
  3. 累计求和所有的 \(x(m) (h(n-m))_L\)
  4. 二者做矩阵乘法运算

4 点和 8 点的循环卷积如下所示:

从上面图示的移位加权累计求和法,可以清晰看到

\(L\) 点循环卷积的结果,当 \(L \ge N+M-1\) 时(\(N\)\(M\) 分别是 \(h(n)\)\(x(n)\) 的长度),循环卷积周期延拓无混叠,此时循环卷积和线性卷积的结果完全一致

4. 用 FFT 计算卷积

我们也可以利用 FFT 计算卷积,因为两个序列时域卷积等于频域相乘。因此我们可以将两个序列分别用 FFT 计算得到频域值,然后相乘,再做 IFFT 即可得到卷积结果。

对于两个长度为 \(N\) 的序列:

  1. 直接卷积的代价为 \(\approx N^2\) 乘法和加法运算
  2. 而一次 FFT 运算为 \(N \cdot \log (N)\) 次乘法和加法

下表为两个方式的运算次数比较:

N FFT Direct Convolution
4 176 16
32 2560 1024
64 5888 4096
128 13,312 16,384
256 29,696 65,536
2048 311,296 4,194,304

发表评论

电子邮件地址不会被公开。 必填项已用*标注