![]() |
Home · Pages · Index · Overviews |
Fast Fourier Transform Module ReferenceFast fourier transforms, convolution, correlation, and normalized correlation. More... #include <fft.h> Routines
Detailed DescriptionThe fast fourier transform (FFT) provides a fast calculation of convolution and correlations (among other things). Given a vector of n values, A = [ a0, a1, ... an-1 ], consider the polynomial PA(x) = an-1xn-1 + an-2xn-2 + ... + a0. The discrete fourier transform of A is the vector FFT(A) = [ f0, f1, ... fn-1 ], where for each i, fi = PA(φi), where φ is the nth root of unity. In other words, the discrete fourier transform of A is the values of the polynomial PA(x) evaluated at the n roots of unity φi for i ∈ [0,n-1]. The key observation is that the pairwise product of the elements of FFT(A) and FFT(B) for two n-element vectors are then the values of the product of the polynomials PA(x) and PB(x) at the n roots of unity, and the coefficients of the product of the two polynomials are the convolutions of the coefficients of said polynomials. Ergo, inverting the transform gives the convolution for each offset of A with respect to b, where the inversion is simply a matter of evaluating PFFT(A)(x) at each of the inverse roots of unity φ-i for i ∈ [0,n-1]. The kicker is that the FFT and its inverse can be computed in O(nlogn) time as opposed to O(n2) time because of the structure of the roots of unity in the field of complex numbers. Correlations are closely related to convolutions and its not hard to show that they amount to taking the product of FFT(A) and the conjugate values of FFT(B). Moreover, the transform can be extended to arrays of 2 or more dimensions. A more traditional view of the FFT is to consider the values A as being in the spatial domain, and then if one thinks of roots of unities as frequencies (a thing physicist find rather useful) FFT(A) is the same object described in the frequency domain. The (forward) FFT maps from the spatial to the frequency domain, and the inverse FFT maps from the frequency domain back to the spatial domain. The mappings are rightly called transforms because A = FFT -1(FFT(A)), i.e. the inverse transform of the transform restores the object. The module provides a single routine FFT that maps an n-dimensional Complex_Array in either the forward or reverse direction depending on the setting of a boolean argument invert. The FFT is performed in-place. It executes optimized codes in the event the argument is real-valued and/or 1-dimensional and/or 32-bit floating point versus 64-bit floating point. For example, FFT(FFT(data,0),1) leaves data unalterred save for any loss of precision due to rounding error (which is typically on the order of 1/2 the precision of the numbers). Unfortunately, Mylib's implementation while matching the best in class for FFT's where dimensions are powers of 2, only works on arrays that have dimensions that are powers of 2. Power_Of_2_Pad returns the next largest power of 2 as a convenience, and the the routines Pad_Array or Pad_Array_Inplace in the Array class will pad the arrays for you as needed. There is a routine FFT_Convolution that computes the pairwise complex product of two FFTs, and another FFT_Correlation that computes the pairwise complex product with the conjugate of the second FFT. Then as expected, given two arrays a and b in the spatial domain, their convolution c and correlation d in the spatial domain can be computed as follows: c = FFT(FFT_Convolution(FFT(Copy_Array(a),0),FFT(Copy_Array(b),0)),1); d = FFT(FFT_Correlation(FFT(Copy_Array(a),0),FFT(Copy_Array(b),0)),1);
Note carefully, that the computed convolution and correlation "wrap"
around the boundary. For example, in the case of 1-dimensional arrays, if a, b, c,
and d are N-element vectors, then:
Because of this wrapping behavior, one often needs to zero-pad the subject arrays in order
to eliminate the effect of wrapping. That is, if a is of length A, and b is of length
B, then zero-padding a by B-1 zeros on the right, and b by A-1 zeros on the left,
will give a convolution or correlation vector of length A+B-1 for which the terms are
the convolution of correlations of a versus b at every possible offset/lag. For example,
for the case of the correlation d, with the given padding it follows that:
Numeric_Array *Zero_Pad(Numeric_Array *a) { Coordinate *shp = AForm_Shape(a); Coordinate *anc = Copy_Array(shp); int i; for (i = 0; i < a->ndims; i++) { ADIMN(shp)[i] *= 2; ADIMN(anc)[i] = 0; } Use_Zero_Boundary(); return (Pad_Array(a,anc,shp)); } Numeric_Array *NCC(Numeric_Array *a, Numeric_Array *b, boolean wrap) { Numeric_Array *ap, *bp; if (wrap) { ap = Copy_Array(a); bp = Copy_Array(b); } else { ap = Zero_Pad(a); bp = Zero_Pad(b); } FFT(FFT_Correlation(FFT(ap,0),FFT(bp,0)),1); Free_Array(bp); return (Normalize_FFT_Correlation(a,b,ap)); } Routine Documentation
Returns the smallest power of 2 not less than m.
FFT of a complex- or real-valued float or double array (see the definition of Complex_Array). If invert is true then the inverse FFT is performed. The operation is performed in-place within the data vector of data, and a pointer to said is returned as a convenience. The code works for arrays of any dimensionality. The operations only work on arrays whose dimensions are powers of 2! This routine is mostly a stub that calls various routines in the MY_FFT sub-library, that have been highly optimzed for various scenarios, i.e. real versus complex, 1-dimensional versus multi-dimensional. It is imperative to note that taking the FFT of a real-valued array, i.e. PLAIN_KIND and either FLOAT32_TYPE or FLOAT64_TYPE, results in a PLAIN_KIND array that is an encoding of a complex, conjugate symmetric array that results when one transforms a real array. This setting of the array's kind, indicates to the routines of the module that are expecting an array in the frequency domain, that the array is so encoded. The reason for the encoding, is that without it, a real valued fft could not be performed in place. See the main section above for the details on how the symmetries of a conjugate symmetric matrix are used to pack its value into an encoding that occupies half as much space.
Given two ffts, fft1 and fft2, of the same shape in the frequency domain, produce their convolution (in the frequency domain), in-place within fft1 and return a pointer to it as a convenience. The routine works just fine if fft1 and fft2 point to the same array. However, both must be complex or both must be an encoding of the fft of a real array. In general, given two arrays a and b in the spatial domain, i.e. untransformed, the way to compute their convolution in the spatial domain is via the template: c = FFT(FFT_Convolution(FFT(Copy_Array(a),0),FFT(Copy_Array(b),0)),1);
If you don't mind a or b being overwritten, then you can remove the respective
calls to
Copy_Array.
Note carefully, that the computed convolution "wraps"
around the boundary. For example, in the case of 1-dimensional arrays, if c[0..N] is the
result of convoluting
a[0..N] and b[0..N] in the spatial domain, then:
Given two ffts, fft1 and fft2, of the same shape in the frequency domain, produce their correlation (in the frequency domain), in-place within fft1 and return a pointer to it as a convenience. The routine works just fine if fft1 and fft2 point to the same array. However, both must be complex or both must be an encoding of the fft of a real array. In general, given two arrays a and b in the spatial domain, i.e. untransformed, the way to compute their correlation in the spatial domain is via the template: c = FFT(FFT_Correlation(FFT(Copy_Array(a),0),FFT(Copy_Array(b),0)),1);
If you don't mind a or b being overwritten, then you can remove the respective
calls to
Copy_Array.
Note carefully, that the computed correlation "wraps"
around the boundary. For example, in the case of 1-dimensional arrays, if c[0..N] is the
result of correlating
a[0..N] and b[0..N] in the spatial domain, then:
Normalize the spatial correlation array cor that is presumed to be the correlation matrix of versions of ref1 and ref2 that were zero-padded to the right to give the dimensions of cor. Normally, each dimension of cor is the sum of the corresponding dimensions in ref1 and ref2 less one in order that the correlation not wrap values about the domain, but this need not be the case.
To explain this further, consider the 1-dimensional case, and suppose cor is of length
C, ref1 of length R1, and ref2 of length R2, and that C = R1+R2-1.
Then, if cor is the correlation of ref1' which is padded with R2-1 zeros on the right,
and ref2' that is padded with R1-1 zeros on the right, then
In full generality, assume cor[i] is the correlation of ref1 and ref2 over the possibly
cyclic domains I1(i) and I2(i) depending on i and the amount of zero-padding.
The normalization is exactly that discussed in the Correlation and Regression module, namely:
|