Two applications of the FFT that are frequently encountered especially in the signal processing area are the discrete convolution and discrete correlation operations.
Given two functions x(t) and y(t), the Fourier transform of the convolution of x(t) and y(t), denoted as x * y, is the product of their individual Fourier transforms: DFT (x * y)=X (•) Y where * denotes the convolution operation and (•) denotes pointwise multiplication.
Typically, x(t) is a continuous and periodic signal that is represented discretely by a set of N data points xj, j = 0, …, N -1, sampled over a finite duration, usually for one period of x(t) at equal intervals. y(t) is usually a response that starts out as zero, peaks to a maximum value, and then returns to zero. Discretizing y(t) at equal intervals produces a set of N data points, yk, k = 0, …, N -1. If the actual number of samplings in yk is less than N, the data can be padded with zeros. The discrete convolution can then be defined as
The values of
are the same as those of k = 0, …, N -1 but in the wrap-around order.
The Oracle Developer Studio Performance Library routines enable you to compute the convolution by using the definition above with k = 0, …, N -1, or by using the FFT. If the FFT is used to compute the convolution of two sequences, the following steps are performed:
Compute X = forward FFT of x
Compute Y = forward FFT of y
Compute Z = X (•) Y <=> DFT (x * y)
Compute z = inverse FFT of Z; z = (x * y)
One interesting characteristic of convolution is that the product of two polynomials is actually a convolution. A product of an m-term polynomial
and an n-term polynomial
has m+n-1 coefficients that can be obtained by:
where k = 0, …, m + n - 2.
Closely related to convolution is the correlation operation. It computes the correlation of two sequences directly superposed or when one is shifted relative to the other. As with convolution, we can compute the correlation of two sequences efficiently as follows using the FFT:
Compute the FFT of the two input sequences.
Compute the pointwise product of the resulting transform of one sequence and the complex conjugate of the transform of the other sequence.
Compute the inverse FFT of the product.
The routines in the Performance Library also allow correlation to be computed by the following definition:
There are various ways to interpret the sampled input data of the convolution and correlation operations. The argument list of the convolution and correlation routines contain parameters to handle the following cases:
The signal and/or response function can start at different sampling times.
You might want only part of the signal to contribute to the output.
The signal and/or response function can begin with one or more zeros that are not explicitly stored.
Oracle Developer Studio Performance Library contains the convolution routines shown in Table 19.
|
The [S,D,C,Z]CNVCOR routines are used to compute the convolution or correlation of a filter with one or more input vectors. The [S,D,C,Z]CNVCOR2 routines are used to compute the two-dimensional convolution or correlation of two matrices.
The one-dimensional convolution and correlation routines use the arguments shown in Table 20.
|
The two-dimensional convolution and correlation routines use the arguments shown in Table 21.
|
The minimum dimensions for the WORK work arrays used with the one-dimensional and two-dimensional convolution and correlation routines are shown in Table 24. The minimum dimensions for one-dimensional convolution and correlation routines depend upon the values of the arguments NPRE, NX, NY, and NZ.
The minimum dimensions for two-dimensional convolution and correlation routines depend upon the values of the arguments shown in the table Table 22.
|
MYC_INIT and NYC_INIT depend upon the following, where X is the filter matrix and Y is the input matrix.
|
The values assigned to the minimum work array size is shown in Table 24.
|
Memory will be allocated within the routine if the workspace size,
indicated by LWORK, is not large enough.
The following example uses CCNVCOR to perform FFT convolution of two complex vectors.
Example 19 One-Dimensional Convolution Using Fourier Transform Method and COMPLEX Datamy_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
If any vector overlaps a writable vector, either because of argument aliasing or ill-chosen values of the various INC arguments, the results are undefined and can vary from one run to the next.
The most common form of the computation, and the case that executes fastest, is applying a filter vector X to a series of vectors stored in the columns of Y with the result placed into the columns of Z. In that case, INCX = 1, INC1Y = 1, INC2Y ≥ NY, INC1Z = 1, INC2Z ≥ NZ. Another common form is applying a filter vector X to a series of vectors stored in the rows of Y and store the result in the row of Z, in which case INCX = 1, INC1Y ≥ NY, INC2Y = 1, INC1Z ≥ NZ, and INC2Z = 1.
Convolution can be used to compute the products of polynomials. The following example Example 20 uses SCNVCOR to compute the product of 1 + 2x + 3x2 and 4 + 5x + 6x2.
Example 20 One-Dimensional Convolution Using Fourier Transform Method and REAL Datamy_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.
Making the output vector longer than the input vectors, as in the example above, implicitly adds zeros to the end of the input. No zeros are actually required in any of the vectors, and none are used in the example, but the padding provided by the implied zeros has the effect of an end-off shift rather than an end-around shift of the input vectors.
The following example Example 21 computes the product between the vector [ 1, 2, 3 ] and the circulant matrix defined by the initial column vector [ 4, 5, 6 ].
Example 21 Convolution Used to Compute the Product of a Vector and Circulant Matrixmy_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.
The difference between the following example and the previous example is that the length of the output vector is the same as the length of the input vectors, so there are no implied zeros on the end of the input vectors. With no implied zeros to shift into, the effect of an end-off shift from the previous example does not occur and the end-around shift results in a circulant matrix product.
Example 22 Two-Dimensional Convolution Using Direct Methodmy_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