卷积和相关

语言：

经常遇到的 FFT 的两个应用（尤其是在信号处理区域中）是离散卷积和离散相关运算。

卷积运算

假设有两个函数 x(t) 和 y(t)，x(t) 和 y(t) 的卷积的傅里叶变换（表示为 x * y）是各个傅里叶变换的积：DFT (x * y)=X (·) Y，其中 * 表示卷积运算，(·) 表示逐点乘法。

通常，x(t) 是由一组 N 个数据点 x_j, j = 0, …, N -1 以离散方式表示的连续和周期性信号，这些数据点是在有限持续时间内取样的，通常是以等间隔在 x(t) 的一个周期内取样。y(t) 通常是以零开始的响应，峰值达到最大值，然后返回到零。在相等间隔对 y(t) 进行离散化处理将生成一组 N 个数据点，y_k, k = 0, …, N -1。如果 y_k 中的实际取样数量小于 N，则使用零填充数据。然后可将离散卷积定义为

image:(x * y) sub {j} equiv sum to {N over 2} from {k={-N over 2} + 1} x sub{j-k} y sub{k}, j = 0, ..., N-1

image:y sub {k},k=-N over 2 + 1, ..., N over 2

的值与 k = 0, …, N -1 的值相同，但使用环绕顺序。

Oracle Developer Studio 性能库例程允许使用上述定义通过 k = 0, …, N -1 或使用 FFT 来计算卷积。如果使用 FFT 计算两个序列的卷积，请执行以下步骤：

计算 X = x 的正向 FFT
计算 Y = y 的正向 FFT
计算 Z = X (·) Y <=> DFT (x * y)
计算 z = Z 的逆向 FFT；z = (x * y)

卷积的一个有趣的特性是两个多项式的积实际上是一个卷积。m 项多项式

a(x) = a₀ + a₁x + ... + a_m-1x^m-1

和 n 项多项式

b(x) = b₀ + b₁x + ... + b_n-1x^n-1

的积具有 m+n-1 系数，可通过以下方程获得该系数：

image:c sub {k} ~= sum to {

其中 k = 0, …, m + n - 2。

Oracle Developer Studio 性能库卷积和相关例程

Oracle Developer Studio 性能库包含如表 19所示的卷积例程。

表 19 卷积和相关例程

例程	参数	功能
`SCNVCOR`、`DCNVCOR`、`CCNVCOR`、`ZCNVCOR`	`CNVCOR,FOUR,NX,X,IFX, INCX,NY,NPRE,M,Y,IFY, INC1Y,INC2Y,NZ,K,Z, IFZ,INC1Z,INC2Z,WORK, LWORK`	包含一个或多个向量的过滤器的卷积或相关
`SCNVCOR2`、`DCNVCOR2`、`CCNVCOR2`、`ZCNVCOR2`	`CNVCOR,METHOD,TRANSX, SCRATCHX,TRANSY, SCRATCHY,MX,NX,X,LDX, MY,NY,MPRE,NPRE,Y,LDY, MZ,NZ,Z,LDZ,WORKIN, LWORK`	两个矩阵的二维卷积或相关
`SWIENER`、`DWIENER`	`N_POINTS,ACOR,XCOR, FLTR,EROP,ISW,IERR`	两个信号的维纳解卷积

[S,D,C,Z]CNVCOR 例程用于计算包含一个或多个输入向量的过滤器的卷积或相关。[S,D,C,Z]CNVCOR2 例程用于计算两个矩阵的二维卷积或相关。

卷积和相关例程的参数

一维卷积和相关例程使用表 20中显示的参数。

表 20 一维卷积和相关例程 SCNVCOR、DCNVCOR、CCNVCOR 和 ZCNVCOR 的参数

参数	定义
`CNVCOR`	`'V'` 或 `'v'` 指定要计算卷积。`'R'` 或 `'r'` 指定要计算相关。
`FOUR`	`'T'` 或 `'t'` 指定要使用傅里叶变换方法。`'D'` 或 `'d'` 指定要使用直接方法，即根据卷积和相关的定义计算卷积或相关。如果要进行卷积的两个序列的长度相同，FFT 方法要比直接方法快。但是，如果一个序列比另一个大得多，例如，在将大时序信号与小过滤器进行卷积时，直接方法的执行速度比基于 FFT 的方法快。
`NX`	过滤器向量的长度，其中 `NX`≥ 0。
`X`	过滤器向量
`IFX`	`X` 的第一个元素的索引，其中 `NX` ≥ `IFX` ≥ 1
`INCX`	`X` 中向量的元素之间的跨距，其中 `INCX` > 0。
`NY`	输入向量的长度，其中 `NY`≥ 0。
`NPRE`	放在 `Y` 向量前面的隐式零的数量，其中 `NPRE` ≥ 0。
`M`	输入向量的数量，其中 `M`≥ 0。
`Y`	输入向量。
`IFY`	`Y` 的第一个元素的索引，其中 `NY` ≥ `IFY` ≥ 1
`INC1Y`	`Y` 中输入向量的元素之间的跨距，其中 `INC1Y` > 0。
`INC2Y`	`Y` 中输入向量之间的跨距，其中 `INC2Y` > 0。
`NZ`	输出向量的长度，其中 `NZ`≥ 0。
`K`	`Z` 向量的数量，其中 `K` ≥ 0。如果 `K` < `M`，则仅处理前 `K` 个向量。如果 `K` > `M`，则处理所有输入向量，并在退出时将最后 `M`-`K` 个输出向量设置为零。
`Z`	结果向量
`IFZ`	`Z` 的第一个元素的索引，其中 `NZ` ≥ `IFZ` ≥ 1
`INC1Z`	`Z` 中输出向量的元素之间的跨距，其中 `INCYZ` > 0。
`INC2Z`	`Z` 中输出向量之间的跨距，其中 `INC2Z` > 0。
`WORK`	工作数组
`LWORK`	工作数组的长度。

二维卷积和相关例程使用表 21中显示的参数。

表 21 二维卷积和相关例程 SCNVCOR2、DCNVCOR2、CCNVCOR2 和 ZCNVCOR2 的参数

参数	定义
`CNVCOR`	`'V'` 或 `'v'` 指定要计算卷积。`'R'` 或 `'r'` 指定要计算相关。
`METHOD`	`'T'` 或 `'t'` 指定要使用傅里叶变换方法。`'D'` 或 `'d'` 指定要使用直接方法，即根据卷积和相关的定义计算卷积或相关。如果要进行卷积的两个矩阵的大小相同，FFT 方法要比直接方法快。但是，如果一个序列比另一个大得多，例如，在将大数据集与小过滤器进行卷积时，直接方法的执行速度比基于 FFT 的方法快。
`TRANSX`	`'N'` 或 `'n'` 指定 `X` 是过滤器矩阵，`'T'` 或 `'t'` 指定 `X` 的转置是过滤器矩阵
`SCRATCHX`	`'N'` 或 `'n'` 指定必须保留 `X`，`'S'` 或 `'s'` 指定可将 `X` 用于暂存空间。在从将 `X` 用于暂存空间的调用返回后，`X` 的内容为未定义。
`TRANSY`	`'N'` 或 `'n'` 指定 `Y` 是输入矩阵，`'T'` 或 `'t'` 指定 `Y` 的转置是输入矩阵
`SCRATCHY`	`'N'` 或 `'n'` 指定必须保留 `Y`，`'S'` 或 `'s'` 指定可将 `Y` 用于暂存空间。在从调用（其中 `X` 用于暂存空间）返回后，`Y` 的内容为未定义。
`MX`	过滤器矩阵 `X` 的行数，其中 `MX` ≥ 0
`NX`	过滤器矩阵 `X` 的列数，其中 `NX` ≥ 0
`X`	过滤器矩阵。当 `SCRATCHX` 是 `'N'` 或 `'n'` 时，`X` 在退出时保持不变；当 `SCRATCHX` 是 `'S'` 或 `'s'` 时，它在退出时将为未定义。
`LDX`	包含过滤器矩阵 `X` 的数组的主维。
`MY`	输入矩阵 `Y` 的行数，其中 `MY` ≥ 0。
`NY`	输入矩阵 `Y` 的列数，其中 `NY` ≥ 0
`MPRE`	放在输入矩阵 `Y` 向量每行前面的隐式零的数量，其中 `MPRE` ≥ 0。
`NPRE`	放在输入矩阵 `Y` 每列前面的隐式零的数量，其中 `NPRE` ≥ 0。
`Y`	输入矩阵。当 `SCRATCHY` 是 `'N'` 或 `'n'` 时，`Y` 在退出时保持不变；当 `SCRATCHY` 是 `'S'` 或 `'s'` 时，它在退出时将为未定义。
`LDY`	包含输入矩阵 `Y` 的数组的主维。
`MZ`	输出向量的数量，其中 `MZ` ≥ 0。
`NZ`	输出向量的长度，其中 `NZ` ≥ 0。
`Z`	结果向量
`LDZ`	包含结果矩阵 `Z` 的数组的主维，其中 `LDZ` ≥ `MAX(1,MZ)`。
`WORKIN`	工作数组
`LWORK`	工作数组的长度。

卷积和相关例程的工作数组 `WORK`

表 24中显示了与一维和二维卷积和相关例程一起使用的 WORK 工作数组的最小维数。一维卷积和相关例程的最小维数取决于参数 NPRE、NX、NY 和 NZ 的值。

二维卷积和相关例程的最小维数取决于表 22中显示的参数的值。

表 22 影响二维例程的最小工作数组大小的参数：SCNVCOR2、DCNVCOR2、CCNVCOR2 和 ZCNVCOR2

参数	定义
`MX`	过滤器矩阵中的行数
`MY`	输入矩阵中的行数
`MZ`	输出向量数
`NX`	过滤器矩阵中的列数
`NY`	输入矩阵中的列数
`NZ`	输出向量长度
`MPRE`	放在输入矩阵每行前面的隐式零的数量
`NPRE`	放在输入矩阵每列前面的隐式零的数量
`MPOST`	`MAX(0,MZ-MYC)`
`NPOST`	`MAX(0,NZ-NYC)`
`MYC`	`MPRE` + `MPOST` + `MYC_INIT`，其中 `MYC_INIT` 取决于过滤器和输入矩阵，如表 23中所示
`NYC`	`NPRE` + `NPOST` + `NYC_INIT`，其中 `NYC_INIT` 取决于过滤器和输入矩阵，如表 23中所示

MYC_INIT 和 NYC_INIT 取决于以下矩阵，其中 X 是过滤器矩阵，Y 是输入矩阵。

表 23 MYC_INIT 和 NYC_INIT 依赖项

例程	Y		Transpose(Y)
	X	Transpose(X)	X	Transpose(X)
`MYC_INIT`	`MAX(MX,MY)`	`MAX(NX,MY)`	`MAX(MX,NY)`	`MAX(NX,NY)`
`NYC_INIT`	`MAX(NX,NY)`	`MAX(MX,NY)`	`MAX(NX,MY)`	`MAX(MX,MY)`

表 24中显示了分配给最小工作数组大小的值。

表 24 与卷积和相关例程一起使用的 WORK 工作数组的最小维数和数据类型

例程	最小工作数组大小 (`WORK`)	类型
`SCNVCOR、DCNVCOR`	`4*(MAX(NX,NPRE+NY)` `+ MAX(0,NZ-NY))`	`REAL、REAL*8`
`CCNVCOR`、`ZCNVCOR`	`2*(MAX(NX,NPRE+NY)` `+ MAX(0,NZ‐NY)))`	`COMPLEX、COMPLEX*16`
`SCNVCOR2`^†、`DCNVCOR2`^†	`MY` `+` `NY` `+` `30`	`COMPLEX、COMPLEX*16`
`CCNVCOR2`^†、`ZCNVCOR2`^†	如果 `MY` = `NY`：`MYC + 8` 如果 `MY` ≥ `NY: MYC + NYC + 16`	`COMPLEX、COMPLEX*16`

如果由 LWORK 指定的工作区大小不够大，则在例程中分配内存。

示例程序：卷积

以下示例使用 CCNVCOR 执行两个复数向量的 FFT 卷积。

示例 19 使用傅里叶变换方法和复数数据的一维卷积

my_system% cat con_ex20.f
      PROGRAM TEST
C
      INTEGER           LWORK
      INTEGER           N
      PARAMETER        (N = 3)
      PARAMETER        (LWORK = 4 * N + 15)
      COMPLEX           P1(N), P2(N), P3(2*N-1), WORK(LWORK)
      DATA P1 / 1, 2, 3 /,  P2 / 4, 5, 6 /
C
      EXTERNAL          CCNVCOR
C
      PRINT *, 'P1:'
      PRINT 1000, P1
      PRINT *, 'P2:'
      PRINT 1000, P2
 
      CALL CCNVCOR ('V', 'T', N, P1, 1, 1, N, 0, 1, P2, 1, 1, 1,
     $              2 * N - 1, 1, P3, 1, 1, 1, WORK, LWORK)
C
      PRINT *, 'P3:'
      PRINT 1000, P3
C
 1000 FORMAT (1X, 100(F4.1,' +',F4.1,'i  '))
C
      END
my_system% f95 -dalign con_ex20.f -xlibrary=sunperf
my_system% a.out
 P1:
  1.0 + 0.0i   2.0 + 0.0i   3.0 + 0.0i  
 P2:
  4.0 + 0.0i   5.0 + 0.0i   6.0 + 0.0i  
 P3:
  4.0 + 0.0i  13.0 + 0.0i  28.0 + 0.0i  27.0 + 0.0i  18.0 + 0.0i

如果任何向量与可写向量发生重叠（由于使用了参数别名，或者由于为各个 INC 参数选择了不恰当的值），则结果都将为未定义，并且每次运行时都不同。

最常用的计算形式以及执行最快的情况是将过滤器向量 X 应用于存储在 Y 的列中的一系列向量，并将结果放在 Z 的列中。在这种情况下，INCX = 1，INC1Y = 1，INC2Y ≥ NY，INC1Z = 1，INC2Z ≥ NZ。另一种常用形式是将过滤器向量 X 应用于存储在 Y 的行中的一系列向量，并将结果存储在 Z 的行中，在这种情况下，INCX = 1，INC1Y ≥ NY，INC2Y = 1，INC1Z ≥ NZ，以及 INC2Z = 1。

可使用卷积来计算多项式的积。以下示例示例 20使用 SCNVCOR 计算 1 + 2x + 3x² 和 4 + 5x + 6x² 的积。

示例 20 使用傅里叶变换方法和实数数据的一维卷积

my_system% cat con_ex21.f
      PROGRAM TEST
      INTEGER     LWORK, NX, NY, NZ
      PARAMETER  (NX = 3)
      PARAMETER  (NY = NX)
      PARAMETER  (NZ = 2*NY-1)
      PARAMETER  (LWORK = 4*NZ+32)
      REAL        X(NX), Y(NY), Z(NZ), WORK(LWORK)
C
      DATA X / 1, 2, 3 /,  Y / 4, 5, 6 /, WORK / LWORK*0 /
C
      PRINT 1000, 'X'
      PRINT 1010, X
      PRINT 1000, 'Y'
      PRINT 1010, Y
      CALL SCNVCOR ('V', 'T', NX, X, 1, 1,   
     $NY, 0, 1, Y, 1, 1, 1,   NZ, 1, Z, 1, 1, 1, WORK, LWORK)
      PRINT 1020, 'Z'
      PRINT 1010, Z
 1000 FORMAT (1X, 'Input vector ', A1)
 1010 FORMAT (1X, 300F5.0)
 1020 FORMAT (1X, 'Output vector ', A1)
      END
my_system% f95 -dalign con_ex21.f -library=sunperf
my_system% a.out
 Input vector X
    1.   2.   3.
 Input vector Y
    4.   5.   6.
 Output vector Z
    4.  13.  28.  27.  18.

使输出向量比输入向量长，如上例所示，可将零隐式添加到输入末尾。实际上，零不是在任何向量中都是必需的，在示例中也未使用零，但使用隐式零进行填充会产生输入向量的舍入移位的效果，而不是循环移位效果。

以下示例示例 21将计算向量 [ 1, 2, 3 ] 和初始列向量 [ 4, 5, 6 ] 定义的循环矩阵之间的积。

示例 21 用于计算向量和循环矩阵的积的卷积

my_system% cat con_ex22.f 
      PROGRAM TEST
C
      INTEGER     LWORK, NX, NY, NZ
      PARAMETER  (NX = 3)
      PARAMETER  (NY = NX)
      PARAMETER  (NZ = NY)
      PARAMETER  (LWORK = 4*NZ+32)
      REAL        X(NX), Y(NY), Z(NZ), WORK(LWORK)
C
      DATA X / 1, 2, 3 /,  Y / 4, 5, 6 /, WORK / LWORK*0 /
C
      PRINT 1000, 'X'
      PRINT 1010, X
      PRINT 1000, 'Y'
      PRINT 1010, Y
      CALL SCNVCOR ('V', 'T', NX, X, 1, 1,   
     $NY, 0, 1, Y, 1, 1, 1,   NZ, 1, Z, 1, 1, 1, 
     $WORK, LWORK)
      PRINT 1020, 'Z'
      PRINT 1010, Z
C
 1000 FORMAT (1X, 'Input vector ', A1)
 1010 FORMAT (1X, 300F5.0)
 1020 FORMAT (1X, 'Output vector ', A1)
      END
my_system% f95 -dalign con_ex22.f -library=sunperf
my_system% a.out
 Input vector X
    1.   2.   3.
 Input vector Y
    4.   5.   6.
 Output vector Z
   31.  31.  28.

下面的示例和前面的示例之间的区别在于，输出向量的长度与输入向量的长度相同，所以在输入向量的末尾没有隐式零。由于没有可以移位到其中的隐式零，因此不会产生前面示例中的舍入移位的效果，而循环移位将产生循环矩阵积。

示例 22 使用直接方法的二维卷积

my_system% cat con_ex23.f
      PROGRAM TEST
C
      INTEGER           M, N
      PARAMETER        (M = 2)
      PARAMETER        (N = 3)
C
      INTEGER           I, J
      COMPLEX           P1(M,N), P2(M,N), P3(M,N)
      DATA P1 / 1, -2, 3, -4, 5, -6 /,  P2 / -1, 2, -3, 4, -5, 6 /
      EXTERNAL          CCNVCOR2
C
      PRINT *, 'P1:'
      PRINT 1000, ((P1(I,J), J = 1, N), I = 1, M)
      PRINT *, 'P2:'
      PRINT 1000, ((P2(I,J), J = 1, N), I = 1, M)
C
      CALL CCNVCOR2 ('V', 'Direct', 'No Transpose X', 'No Overwrite X',
     $   'No Transpose Y', 'No Overwrite Y', M, N, P1, M,
     $   M, N, 0, 0, P2, M, M, N, P3, M, 0, 0)
C
      PRINT *, 'P3:'
      PRINT 1000, ((P3(I,J), J = 1, N), I = 1, M)
C
 1000 FORMAT (3(F5.1,' +',F5.1,'i  '))
C
      END
my_system% f95 -dalign con_ex23.f -library=sunperf
my_system% a.out
 P1:
  1.0 +  0.0i    3.0 +  0.0i    5.0 +  0.0i  
 -2.0 +  0.0i   -4.0 +  0.0i   -6.0 +  0.0i  
 P2:
 -1.0 +  0.0i   -3.0 +  0.0i   -5.0 +  0.0i  
  2.0 +  0.0i    4.0 +  0.0i    6.0 +  0.0i  
 P3:
-83.0 +  0.0i  -83.0 +  0.0i  -59.0 +  0.0i  
 80.0 +  0.0i   80.0 +  0.0i   56.0 +  0.0i