卷积运算和程序实现

乌雷尼奥

已于 2023-10-22 21:26:56 修改

阅读量2.1k

点赞数 28

文章标签： c++ 算法

于 2023-08-20 19:50:16 首次发布

本文链接：https://blog.csdn.net/qq_45938795/article/details/132327493

版权

两个离散序列的卷积运算

对两个离散的数据序列进行卷积运算是学习和工作中常会遇见的事情。

这里首先给出卷积的多项式表示并给出卷积的直观概念，进一步从函数图像来解释“卷”和“积”。

从多项式乘法的系数来获得卷积结果

假设有两个离散的数据序列：

DataSquence_1 = {1, 2, 3};

DataSquence_2 = {4, 5, 6};

我们可以给出两个多项式，并将离散的数据序列视为多项式的系数，如下：

$Polynmial_1 = 1+2x+3x^{2}$

$Polynmial_2 = 4+5x+6x^{2}$

令两多项式的积为多项式 MultiPolynomial，则：

$MultiPolynomial = (1+2x+3x^{2})\cdot (4+5x+6x^{2})$

$=4+5x+6x^{2}+8x+10x^{2}+12x^{3}+12x^{2}+15x^{3}+18x^{4}\\$

$=4+13x+28x^{2}+27x^{3}+18x^{4}$

因此 DataSquence_1 ∗ DataSquence_2 = {1, 2, 3} ∗ {4, 5, 6} = { 4, 13, 28, 27, 18 }

同样，可以有两个不等长度的离散数据序列：

DataSquence_3 = {1, 2, 3, 4};

DataSquence_4 = {5, 6, 7};

同样使用多项式的积来获得卷积运算的结果：

$MultiPolynomial = (1+2x+3x^{2}+4x^{3})\cdot (5+6x+7x^{2})$

$=5x+6x+7x^{2}+10x+12x^{2}+14x^{3}+15x^{2}+x^{2}18x^{3}+21x^{4}$

$+20x^{3}+24x^{4}+28x^{5}$

$=5x+16x+34x^{2}+52x^{3}+45x^{4}+28x^{5}$

因此 DataSquence_3 ∗ DataSquence_4 = {1, 2, 3, 4} ∗ {5, 6, 7} = { 5, 16, 34, 52, 45, 28 }

在了解多项式的积可以获得卷积运算结果后，下面将基于此给出卷积的直观概念

首先，两个多项式及两多项式的乘积改写如下：

$Polynmial_1 = f(0)x^{0}+f(1)x+...+f(n-1)x^{n-1}$

$Polynmial_2 = h(0)x^{0}+h(1)x+...+h(m-1)x^{m-1}$

$MultiPolynomial = g(0)x^{0}+g(1)x+...+g(n+m-2)x^{n+m-2}$

也就是说：

{ f(0), f(1), ... , f(n - 1) } ∗ { h(0), h(1), ... , h(m - 1) } = { g(0), g(1), ... , g(n + m - 2 }

卷积结果的通式为：

$g(t)=\sum_{k=0}^{k=max(n,m)}h(k)\cdot f(t-k)$

对此需要注意三点：

1. g(t) 序列，也就是卷积运算的结果序列，其长度是参与卷积运算的两序列长度（n、m）的和再减去

1，即 n + m - 1;

2. 当涉及到并不存在的序列索引时，如上式中 t = 0，k = 1 时，f(-1) 实际并不存在，此时认为序列中

该索引的值为零，即认为 f(-1) 为0；

3. 上式中 n 和 m 的数值可以不同，但其表达依旧与卷积表达式的定义十分类似，区别在于上式为离散

形式。

从函数图像来了解卷积运算

以上论述，由多项式乘积的角度介绍了 “卷积” 的计算方法。实际上，在 g(t) 的计算中就体现了“卷”和“积”的特性，而下面将在连续函数上进一步解释。

下式为函数 f(x) 和 h(x) 的卷积公式：

$f(x)\ast h(x)=\int_{-\infty }^{+\infty}f(x^{'})\cdot h(x-x^{'})dx^{'}$

$=g(x)$

图一展示了上面卷积公式中的函数 f(x’) 和 h(x-x') 的示意图像。

图1. 函数 f(x) 和 h(x)示意图。左：f(x‘) 和 h(x’)的示意图；中：f(x‘) 和 h(-x’)的示意图；右：f(x‘) 和 h(x-x’)的示意图

函数卷积公式中的积分符号表示：图一中右图所示两函数图像在x’轴上对应的x’位置处函数值的乘积在 x’ 有效区间内的和。

图一直观表示了卷积运算中 “卷” 的概念，也就是翻转、平移。

对函数图像进行采样，可以将函数卷积过程可视化为离散序列的卷积如下。

其中的连线表示采样的函数值相乘，绿色虚线表示有限区间的离散序列无法索引的值，同上一小节一样认为无法索引的函数值为0：

图2. 对 f(x') 和 h(x-x') 采样序列的卷积。左： $X_{1}$ 处的卷积值；右： $X_{2}$ 处的卷积值

图二绿色连线表示了卷积运算中 “积” 的概念，也就是乘积、求和。

与此同时，可以将 f(x) 和 h(x) 的原始图像在X1位置处卷积的可视化进行对比以加深对卷积运算的理解，如下：

图3. 对 f(x) 和 h(x) 采样序列的卷积

就通过函数图像了解卷积来说，需要注意两点：

1. 函数的卷积为连续域上的积分；

2. 由图2和3可以看到，卷积运算结果在某个 x 位置处的值与参与卷积运算的两个函数在整个积分区间

上的函数值都相关。

卷积运算的程序实现

这部分将展示卷积运算的程序实现。

CT的卷积反投影涉及到投影数据与滤波函数的卷积运算，以此为例，将进行展示如下：

图4. 投影数据（ $p(t)$ ）与滤波函数（ $h_{SL}(t)$ ）

其中，滤波函数选取为S-L滤波函数，其离散化的序列表达如下，其中Δ为投影数据的采样间隔：

$h_{SL}(t )=\frac{-2}{\pi ^{2} \Delta^{2}(4t^{2}-1)}$

卷积运算代码如下：

// 卷积运算函数
void conv(double *pfSequence_1, double *pfSequence_2, int nSequence_Len, double *pfResult_Squence)
{
    int nTotal_Length = nSequence_Len + nSequence_Len - 1;
    int nCore_Length = nSequence_Len;
    int nCore_Start;
    if (((nTotal_Length - nCore_Length) % 2) == 0)
    {
        nCore_Start = (nTotal_Length - nCore_Length) / 2;
    }
    else
    {
        nCore_Start = 1 + (nTotal_Length - nCore_Length) / 2;
    }
    
    // 卷积核心算法部分
    double *pfTarget = new double[nTotal_Length]();
    for(int i = 0; i < nTotal_Length; i++)
    {
        for(int j = max(0, i + 1 - nSequence_Len); j <= min(i, nSequence_Len -1 ); j++)
        {
            pfTarget[i] += pfSequence_1[j] * pfSequence_2[i-j];
            
        }
    }

    // 选取截断部分
    for(int i = 0; i < nCore_Length; i++)
    {
        pfResult_Squence[i] = pfTarget[i + nCore_Start];
    }
}


const int nProj_Channels = 511;
const double fDerta = 0.6416;
const double PI = 3.1415926;

// 投影数据
double *pfProjection = new double[nProj_Channels];
for (int i = 0; i < nProj_Channels; i++)
{
    pfProjection[i] = Projection.Getdata(0, i);  // 投影数据序列

}

// S-L滤波函数的时间域表达（序列）
double *pfFilter_Func = new double[nProj_Channels]();
for (int i = 0; i < nProj_Channels; i++)
{
    int nSample = -(nProj_Channels -1 )/ 2.0 + i;
    pfFilter_Func[i] = -2 / (PI * PI * fDerta * fDerta * (4 * nSample * nSample - 1));
}

// 检查两个离散序列的长度
cout<<_msize(pfFilter_Func)/sizeof(pfFilter_Func[0])<<endl;   // 511
cout<<_msize(pfProjection)/sizeof(pfProjection[0])<<endl;     // 511

// 卷积
double *pfFiltered_Projection = new double[nProj_Channels](); // 卷积运算结果序列
conv(pfFilter_Func, pfProjection, nProj_Channels, pfFiltered_Projection);

// 卷积运算结果序列的长度
cout<<pfFiltered_Projection<<endl;                           // 511

卷积运算结果如下，g(t) 为卷积运算结果，即滤波投影数据：

图5. 投影数据与滤波函数卷积运算得到的滤波投影数据

就卷积运算的程序实现而言，需要注意三点：

1. 本例所展示的情况为：长度相同的两个离散序列的卷积。在需要进行不同长度乃至不同维度序列卷

积运算时，则不适用；

2. 卷积运算结果的长度是参与卷积运算的两个序列长度之和减1；

3. 在参与卷积运算的两个序列长度相同的情况下，一般希望卷积运算结果的长度也相同。

可以将卷积运算结果进行截断，所选取的截断部分通常位于序列的中心位置。本例展示的即为截断

后的卷积运算结果，nCore_Start 为截断起始位置。

//***//

//***//

以 DataSquence_1 ∗ DataSquence_2 为例：

DataSquence_1 ∗ DataSquence_2 = {1, 2, 3} ∗ {4, 5, 6} = { 4, 13, 28, 27, 18 }，

与参与卷积运算序列同长度的卷积结果序列为：{13, 28, 27}。

卷积运算与双倍采样定理

时域相乘后再做傅里叶变换相当于频域卷积。

对某个函数f(x) 进行傅里叶变换，本质上是使用采样函数与f(x) 相乘，接着对采样结果进行离散傅里叶变换，该过程实际相当于f(x)的频域形式F(ω)与采样函数频域形式H(ω)进行卷积。

设采样频率为5，则H(ω)在频域中为下图所示：

图6. 采样频率为5的采样函数频域形式

假设两种情况

1. 设F(ω)包含两个频率分量（频率为1和2的分量，外加直流分量），考虑到傅里叶空间的对称性，F(ω)表示为[2, 3, 4, 3, 2]；

2. 设F(ω)包含三个频率分量，即F(ω)表示为[1, 2, 3, 4, 3, 2, 1]；

H(ω)分别与这两种F(ω)卷积后的结果如图7所示，可以看到F(ω)=[2, 3, 4, 3, 2]时，各频率分量都被完好的保留了下来，而F(ω)=[1, 2, 3, 4, 3, 2, 1]时，发生了混叠现象：

图7. 采样函数H(ω)分别与F(ω)进行卷积的结果

F(ω)=[1, 2, 3, 4, 3, 2, 1]时发生混叠的原因在于F(ω)与H(ω)卷积时，采样函数H(ω)的其他周期也参与了卷积运算，如图8所示。由图8能轻松观察到，如果想要不发生混叠，则需采样函数的频率间隔 f 大于等于F(ω)所包含最大频率的两倍。这就是双倍采样定理，也就是奈奎斯特采样定理。

即F(ω)=[1, 2, 3, 4, 3, 2, 1]时，最大频率为3，f =5 < 3 · 2，因此发生混叠；

F(ω)=[ 2, 3, 4, 3, 2]时，最大频率为2，f =5 > 2 · 2，因此不发生混叠。