convolutioncomplex¶

Computes the convolution, and optionally, the correlation of two complex vectors.

Synopsis¶

convolutionComplex (x, y)

Required Arguments¶

complex x[] (Input): Real vector of length nx.
complex y[] (Input): Real vector of length ny.

Return Value¶

An array of length nz containing the convolution of x and y. If no zeros are computed, then None is returned.

Optional Arguments¶

periodic: The input is periodic.
correlation: Return the correlation of x and y.
firstCall: If the function is called multiple times with the same nx and ny, select this option on the first call.
continueCall: If the function is called multiple times with the same nx and ny, select this option on intermediate calls.
lastCall: If the function is called multiple times with the same nx and ny, select this option on the final call.

Description¶

The function convolutionComplex, by default, computes the discrete convolution of two sequences x and y. More precisely, let \(n_x\) be the length of x, and \(n_y\)denote the length of y. If a circular convolution is desired, the optional argument periodic must be selected. We set

\[n_z = max \{n_y, n_x\}\]

and we pad out the shorter vector with zeros. Then, we compute

\[z_i = \sum_{j=1}^{n_z} x_{i-j+1} y_j\]

where the index on x is interpreted as a positive number between 1 and \(n_z\), modulo \(n_z\).

The technique used to compute the \(z_i\)’s is based on the fact that the (complex discrete) Fourier transform maps convolution into multiplication. Thus, the Fourier transform of z is given by

\[\hat{z}(n) = \hat{x}(n) \hat{y}(n)\]

where the following equation is true.

\[\hat{z}(n) = \sum_{m=1}^{n_z} z_m e^{-2 \pi i(m-1)(n-1)/n_z}\]

The technique used here to compute the convolution is to take the discrete Fourier transform of x and y, multiply the results together component-wise, and then take the inverse transform of this product. It is very important to make sure that \(n_z\) is the product of small primes if option periodic is selected. If \(n_z\) is a product of small primes, then the computational effort will be proportional to \(n_z \log(n_z)\). If option periodic is not selected, then a good value is chosen for \(n_z\) so that the Fourier transforms are efficient and \(n_z \geq n_x + n_y − 1\). This will mean that both vectors will be padded with zeros.

Optionally, the function convolutionComplex computes the discrete correlation of two sequences x and y. More precisely, let n be the length of x and y. If a circular correlation is desired, then option periodic must be selected. We set (on output)

\(n_z = n\)	if `periodic` is chosen
\(n_z = 2^\alpha 3^\beta 5^\gamma \geq 2n - 1\)	if `periodic` is not chosen

where α, β, and γ are nonnegative integers yielding the smallest number of the type \(2^\alpha 3^\beta 5^\gamma\) satisfying the inequality. Once \(n_z\) is determined, we pad out the vectors with zeros. Then, we compute

\[z_i = \sum_{j=1}^{n_z} x_{i+j-1} y_j\]

where the index on x is interpreted as a positive number between one and \(n_z\), modulo \(n_z\). Note that this means that

\[z_{n_z-k}\]

contains the correlation of \(x(k − 1)\) with y as \(k = 0, 1, ..., n_z∕2\). Thus, if \(x(k − 1) = y(k)\) for all k, then we would expect

\[\Re z_{n_z}\]

to be the largest component of \(\Re z\). The technique used to compute the \(z_i\)’s is based on the fact that the (complex discrete) Fourier transform maps correlation into multiplication.

Thus, the Fourier transform of z is given by

\[\hat{z}_j = \hat{x}_j \overline{y}_j\]

where the following equation is true.

\[\hat{z}_j = \sum_{m=1}^{n_z} z_m e^{-2 \pi i(m-1)(j-1)/n_z}\]

Thus, the technique used here to compute the correlation is to take the discrete Fourier transform of x and the conjugate of the discrete Fourier transform of y, multiply the results together component-wise, and then take the inverse transform of this product. It is very important to make sure that \(n_z\) is the product of small primes if periodic is selected. If \(n_z\) is the product of small primes, then the computational effort will be proportional to \(n_z \log(n_z)\). If periodic is not chosen, then a good value is chosen for \(n_z\) so that the Fourier transforms are efficient and \(n_z \geq 2n − 1\). This will mean that both vectors will be padded with zeros.

No complex transforms of x or y are taken since both sequences are real, and real transforms can simulate the complex transform above. Such a strategy is six times faster and requires less space than when using the complex transform.

Examples¶

Example 1¶

This example computes a nonperiodic convolution. The purpose is to compute a moving average of the type found in digital filtering. The averaging operator in this case is especially simple and is given by averaging five consecutive points in the sequence. We try to recover the values of an exponential function contaminated by noise. The large error for the last value has to do with the fact that the convolution is averaging the zeros in the “pad” rather than the function values. Notice that the signal size is 100, but only report the errors at ten points.

from __future__ import print_function
from numpy import *
from pyimsl.math.convolutionComplex import convolutionComplex
from pyimsl.math.constant import constant
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform


def F1(a):
    return complex(exp(a), 0.0) * complex(cos(a), sin(a))


ny = 100
nfltr = 5
pi = constant("pi")
fltr = empty(nfltr, dtype=complex)
y = empty(ny, dtype=complex)

# Set up the filter
for i in range(0, nfltr):
    fltr[i] = complex(0.2, 0.0)

# Set up the y-vector for the periodic case
twopi = 2.0 * pi
randomSeedSet(1234579)
noise = randomUniform(2 * ny)
for i in range(0, ny):
    x = float(i) / (ny - 1)
    temp = complex(0.5 * noise[i] - 0.25, 0.5 * noise[ny + i] - 0.25)
    y[i] = F1(x) + temp

# Call the convolution routine for the periodic case
z = convolutionComplex(fltr, y)

# Print results
print("Periodic Case")
print("        x          F1(x)        Original Error   Filtered Error")
total1 = 0.0
total2 = 0.0
for i in range(0, ny):
    x = float(i) / (ny - 1)
    origer = abs(y[i] - F1(x))
    fltrer = abs(z[i + 2] - F1(x))
    if ((i % 11) == 0):
        print("%10.4f   (%6.4f,%6.4f) %12.4f %15.4f"
              % (x, (F1(x)).real, (F1(x)).imag, origer, fltrer))
    total1 = total1 + origer
    total2 = total2 + fltrer

print(" Average absolute error before filter:%10.5f" % (total1 / (ny)))
print(" Average absolute error after filter:%11.5f" % (total2 / (ny)))

Output¶

Periodic Case
        x          F1(x)        Original Error   Filtered Error
0000   (1.0000,0.0000)       0.1684          0.3524
1111   (1.1106,0.1239)       0.0582          0.0822
2222   (1.2181,0.2752)       0.1991          0.1054
3333   (1.3188,0.4566)       0.1487          0.1001
4444   (1.4081,0.6706)       0.2381          0.1004
5556   (1.4808,0.9192)       0.1037          0.0708
6667   (1.5307,1.2044)       0.1312          0.0904
7778   (1.5508,1.5273)       0.1695          0.0856
8889   (1.5331,1.8885)       0.1851          0.0698
0000   (1.4687,2.2874)       0.2130          1.0760
 Average absolute error before filter:   0.19057
 Average absolute error after filter:    0.10024

Example 2¶

This example computes both a periodic correlation between two distinct signals x and y. There are 100 equally spaced points on the interval \([0, 2\pi]\) and \(f_1(x) = \cos (x) + i \sin (x)\). Define x and y as follows:

\[\begin{split}\begin{array}{l} x_i = f_1 \left(\frac{2\pi(i-1)}{n-1}\right) \phantom{...} i = 1, \ldots, n \\ y_i = f_1 \left(\frac{2\pi(i-1)}{n-1} + \frac{\pi}{2}\right) \phantom{...} i = 1, \ldots, n \end{array}\end{split}\]

Note that the maximum value of z (the correlation of x with) occurs at \(i = 25\), which corresponds to the offset.

from __future__ import print_function
from numpy import *
from pyimsl.math.convolutionComplex import convolutionComplex
from pyimsl.math.constant import constant

n = 100
pi = constant("pi")
x = empty(n, dtype=complex)
y = empty(n, dtype=complex)
zreal = empty(n, dtype=double)

# Set up y-vector for the nonperiodic case
for i in range(0, n):
    x[i] = complex(cos(2.0 * pi * i / (n - 1)),
                   sin(2.0 * pi * i / (n - 1)))
    y[i] = complex(cos(2.0 * pi * i / (n - 1) + pi / 2.0),
                   sin(2.0 * pi * i / (n - 1) + pi / 2.0))

# Call the correlation function for the nonperiodic case
z = convolutionComplex(x, y, correlation=True, periodic=True)

sumx = sumy = 0.0
for i in range(0, n):
    sumx = sumx + abs(x[i] * x[i])
    sumy = sumy + abs(y[i] * y[i])
xnorm = sqrt(sumx)
ynorm = sqrt(sumy)
for i in range(0, n):
    zreal[i] = z[i].real / (xnorm * ynorm)

max = zreal[0]
k = 0
for i in range(1, n):
    if (max < zreal[i]):
        max = zreal[i]
        k = i

print("The element of z with the largest normalized ")
print("value is Z(%2d)." % (k))
print("The normalized value of Z(%2d) is %6.3f" % (k, zreal[k]))

Output¶

The element of z with the largest normalized 
value is Z(25).
The normalized value of Z(25) is  1.000