一篇文章搞懂卷积

本文从一名工程师的角度讲解了数字图像处理领域的重要概念——卷积。本文力求使用最少的数学公式将卷积的物理含义以及计算方法讲解清楚。其中第四节和第五节图解卷积过程，第六节给出了卷积的代码实现。相信工程师朋友们可以通过这一篇文章轻松理解卷积。

Table of Contents

（一）线性时不变系统

（二）单位脉冲信号

（三）单位脉冲响应

（四）卷积计算（视角一）

（五）卷积计算(视角二)

（六）卷积计算代码实现

---

（一）线性时不变系统

如果一个系统具备齐次性，可加性和时不变性，则该系统为线性时不变系统。

齐次性：如果一个系统对输入信号的输出为，那么该系统对输入信号的输出为 。

可加性：如果一个系统对输入信号的输出为，该系统对输入信号的输出为 ，那么该系统对输入信号的输出为。

时不变性：如果一个系统对输入信号的输出为，那么该系统对的输出为

---

（二）单位脉冲信号

如果一个信号 ，只有在n=0时信号大小为1，对于其余所有n，信号大小均为0，则该信号叫做单位脉冲信号。

  

信号 ，在n=-k时信号大小为1，其余所有n信号大小均为0，这相当于将信号 ，沿着n轴向左平移k个单位。

任何离散信号都可以表示为单位脉冲信号的线性组合。如下图所示，信号x[n]可以表示为一组脉冲信号的线性组合：

---

（三）单位脉冲响应

如果已知一个线性时不变系统的单位脉冲响应（输入信号为单位脉冲信号时系统的输出），那我们就知道了线性时不变系统的一切。假设一个线性时不变系统的单位脉冲响应为。即输入信号为的时候，系统的输出为 。则给定任何离散输入信号，我们均能得到对应该输入信号的输出信号。

        （1）因为任何离散信号均可以表示成单位脉冲信号的线性组合。

                

        （2）由线性时不变系统的齐次性和可加性可知，线性时不变系统的输出为：

                

        （3）由线性时不变系统的时不变性可知：

                 

        （4）已知线性时不变系统的单位脉冲响应为 ，对于任意输入信号，系统的输出为：

                

上式即为卷积，即已知线性时不变系统的单位脉冲响应 ，则对于任意输入信号 ，系统的输出等于输入信号 和单位脉冲响应 的卷积，记做：

---

（四）卷积计算（视角一）

该视角从输入信号出发，对于输入信号中的每一个元素 ，将单位脉冲响应 进行相应的平移，得到 。然后使用元素 乘以 ​​​​​​​ 得到一个结果信号。对于输入信号中的每一个元素执行同样的操作。具体计算过程图解如下(信号​​​​​​​和​​​​​​​ 参考后文的图):

    

    

    

当输入信号的所有元素均完成该操作后，将所有的结果信号相加，即可得到​​​​​​​ 和​​​​​​​ 的卷积。将以上7个信号向加，即可得到​​​​​​​ 与​​​​​​​ 的卷积，计算结果如下：

视角一的每次迭代，输入信号中的一个元素就完成了它的使命，输入信号​​​​​​​ 与​​​​​​​ 相乘并将结果信号贡献给输出信号的若干元素。该视角更利于理解卷积的物理含义。

---

（五）卷积计算(视角二)

该视角从输出信号​​​​​​​ 出发，根据视角一中的分析，可以观察到输出信号的每个元素的计算方法如下图所示：

    

        

        

    

    

该视角每次迭代完成输出信号中的一个元素的计算，要注意该方法中信号​​​​​​​ 的反序操作，该视角更方便理解卷积的计算过程。

---

（六）卷积计算代码实现

本节从文中的视角一和视角二两个角度实现卷积的计算。卷积的计算类的代码粘贴如下，该类为C++模板类，之所以设计为模板类是为了使得该类支持任意数值类型的卷积计算。其中成员函数ExecuteInputView对应视角一的算法，成员函数ExecuteOutputView对应视角二的算法。想通过代码更加深入理解卷积计算过程的读者可以仔细研究这两个函数的函数体。

#pragma once

#include <vector>
#include <type_traits>
#include <stdio.h>

template <typename T>
class Convolution
{
public:

    Convolution();

public:

    std::vector<T> ExecuteInputView(const std::vector<T> &x, const std::vector<T> &h);
    std::vector<T> ExecuteOutputView(const std::vector<T> &x, const std::vector<T> &h);

private:

    T ReverseDotProduct(const T *x, const T *y, int len);
};

template <typename T>
Convolution<T>::Convolution()
{
    if (!std::is_arithmetic<T>::value) {
        std::exception e("Convolution template class only support arithmetic type.");
        throw e;
    }
}

template <typename T>
std::vector<T> Convolution<T>::ExecuteInputView(const std::vector<T> &x, const std::vector<T> &h)
{
    int nInputSignalNum = static_cast<int>(x.size());
    int nPulseResponseNum = static_cast<int>(h.size());
    int nOutputSignalNum = nInputSignalNum + nPulseResponseNum - 1;
    std::vector<T> y(nOutputSignalNum, static_cast<T>(0));
    for (int i = 0; i < nInputSignalNum; ++i) {
        for (int j = 0; j < nPulseResponseNum; ++j) {
            y[i + j] += x[i] * h[j];
        }
    }
    return y;
}

template <typename T>
std::vector<T> Convolution<T>::ExecuteOutputView(const std::vector<T> &x, const std::vector<T> &h)
{
    int nInputSignalNum = static_cast<int>(x.size());
    int nPulseResponseNum = static_cast<int>(h.size());
    int nOutputSignalNum = nInputSignalNum + nPulseResponseNum - 1;
    std::vector<T> xEx(nInputSignalNum + 2 * (nPulseResponseNum - 1), static_cast<T>(0));
    memcpy(xEx.data() + nPulseResponseNum - 1, x.data(), sizeof(T)*nInputSignalNum);
    std::vector<T> y(nOutputSignalNum, static_cast<T>(0));
    for (int i = 0; i < nOutputSignalNum; ++i) {
        y[i] = ReverseDotProduct(xEx.data() + i, h.data(), nPulseResponseNum);
    }

    return y;
}

template <typename T>
T Convolution<T>::ReverseDotProduct(const T *x, const T *y, int len)
{
    T res = static_cast<T>(0);
    for (int i = 0, j = len - 1; i < len; ++i, --j) {
        res += x[i] * y[j];
    }

    return res;
}

如下代码是卷积计算类的单元测试代码，感兴趣的读者可以在自己的电脑上运行。C++11以上即可使用。

#include <iostream>

#include "Convolution.h"

using namespace std;

int main()
{
    vector<int> x = { 1,2,1,2,4,2,1 };
    vector<int> h = { 1,2,2,3 };

    try {
        Convolution<int> conv;

        vector<int> yInput = conv.ExecuteInputView(x, h);
        vector<int> yOutput = conv.ExecuteOutputView(x, h);

        vector<int> ref = { 1, 4, 7, 11, 16, 17, 19, 18, 8, 3 };

        if (yInput != ref || yOutput != ref) {
            cerr << "UTConvolution failed." << endl;
            return 1;
        }
    }
    catch (const exception &e) {
        cerr << "UTConvolution exception:" << e.what() << endl;
        return 1;
    }

    cout << "UTConvolution pass." << endl;

    return 0;
}

以上内容如果存在任何错误和不足，欢迎各位朋友在下方留言，我定会认真接受！！！