HELP FFT David Young April 1996, revised Oct 2001 LIB * FFT is an implementation of the Fast Fourier Transform algorithm, and is part of the popvision package. The routines are based on that in "Numerical Recipes in C" by W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling (CUP 1988), p. 411. CONTENTS - (Use g to access required sections) 1 Introduction 1.1 The 1-D Discrete Fourier Transform 1.2 The 2-D Discrete Fourier Transform 2 The Fast Fourier Transform 2.1 The procedure fft_1d 2.2 Array subscripting and array bounds 2.3 Invertibility 2.4 The procedure fft_2d 3 Transforms for real arrays 3.1 Symmetries for real data 3.2 Representing symmetrical data 3.3 Routines for 1-D real arrays 3.4 Routines for 2-D real arrays 3.5 Routine for 2-D real symmetric arrays 4 Convenience routines ----------------------------------------------------------------------- 1 Introduction ----------------------------------------------------------------------- You are strongly advised to consult a textbook dealing with the Fourier Transform to gain familiarity with the mathematical ideas underlying these routines. The Discrete Fourier Transform is the name of the mathematical transform outlined in this section and described in mathematics books. It is the discrete version of the Fourier Transform itself, which is an integral transform for functions of a continuous variable. The Fast Fourier Transform (FFT) is an algorithm which implements the Discrete Fourier Transform. 1.1 The 1-D Discrete Fourier Transform --------------------------------------- Given some data D(n), n = 0..N-1, the 1-dimensional "forward" transform is defined by 1 N-1 F(k) = - * SIGMA D(n) * exp(2 * pi * i * n * k / N) N n=0 where i**2 = -1 and both D(n) and F(k) may be complex. SIGMA means summation over the range of indices given. Note that F is periodic with period N: F(k+j*N) = F(k) for any integer j. Thus if F is calculated for k = 0..N-1, it is defined for all k, using F(k) = F(k mod N). F(k) is therefore usually calculated only for k in the range 0..N-1. The "reverse" transform is N-1 D(n) = SIGMA F(k) * exp(-2 * pi * i * n * k / N) k=0 which will reconstruct D(n) for n in the range 0..N-1. If the reverse transform is used to find D(n) for other values of n, then the reconstructed D is also periodic, with D(n) = D(n mod N) for all n. It is therefore convenient to regard the original data as one cycle of a periodic series - the cycle on 0..N-1 - of which the other cycles are implicit. The normalising factor 1/N in the forward transform is an arbitrary choice, since the only mathematical requirement is that the product of this factor and the factor for the reverse transform (here equal to 1) is 1/N. Sometimes the transforms are made more symmetrical by using a factor of 1/sqrt(N) in each direction, and sometimes the factor of 1/N is applied to the reverse transform. The choice used here has the merit that the components of F can be viewed as the amplitudes of the constituent harmonic waves of the data D. This view breaks the symmetry, and identifies D as the time or space domain and F as the frequency domain. Given this, the component F(k) of the transform is the (complex) amplitude of a harmonic curve of frequency 2*pi*k/N and wavelength N/k in the data domain. In particular, F(0) is just the average of the data, F(1) is the amplitude of the component that has a single cycle in 0..N-1, and F(N/2) is the amplitude of the component that alternates between positive and negative on successive data points. F(N/2) thus describes the highest frequency that can be represented in the data, and F(k) for k > N/2 is best thought of as describing the negative frequency associated with F(k-N). In the implementation given here, N is always a power of 2, for efficiency. As described below, the Pop-11 interface is capable of concealing this constraint, which can be convenient for simple applications. However in applications in which the frequency components are to be interpreted, it is essential to use the true value of N. 1.2 The 2-D Discrete Fourier Transform --------------------------------------- The generalisation to 2 dimensions is straightforward. If D(n, m) is defined on n = 0..N-1, m = 0..M-1, the formula is 1 N-1 M-1 F(k,j) = -- * SIGMA SIGMA D(n, m) * exp(2 * pi * i * (n*k/N + m*j/M)) NM n=0 m=0 but it is important to realise that this can be written as two successive 1-D transforms, thus: 1 N-1 ( F(k,j) = - * SIGMA ( exp(2*pi*i*n*k/N) * N n=0 ( 1 M-1 ) - SIGMA exp(2*pi*i*m*j/M) * D(n,m) ) M m=0 ) With this in mind, all the statements about periodicity that apply to the 1-D transform apply equally to each dimension of the 2-D transform. Thus D and F are periodic in the plane, with the region n = 0..N-1, m = 0..M-1 being repeated over the rest of the plane in both cases. The inverse of the 2-D transform is the straightforward analogue of the 1-D inverse; that is, it is obtained by negating the exponent and omitting the factor 1/NM. The separability is also used to implement the 2-D transform: a 1-D transform is applied to each row of a 2-D array, then to each column of the resulting array (or vice versa). ----------------------------------------------------------------------- 2 The Fast Fourier Transform ----------------------------------------------------------------------- 2.1 The procedure fft_1d ------------------------- fft_1d(N, dir, inr, ini, outr, outi) -> (outr, outi) Carries out a 1-dimensional Fast Fourier Transform on complex data. (See Section 3 for transforms dealing with real data.) N may be an integer giving the transform length as defined above. It must be a power of 2. Alternatively, N may be , in which case the transform length is the number of input data points rounded up to the next power of 2. dir must be for a forward transform and for a backward transform. inr and ini are 1-D arrays with the same bounds, containg the real and imaginary parts respectively of the input data. There are no other restrictions on the bounds - see below. outr and outi may both be input as , in which case new arrays with boundslists [0 N-1], created with *newsfloatarray, are created and returned. Alternatively, outr and outi may be input as 1-D arrays with the same bounds as one another, in which case the results are stored in these arrays and the same arrays are returned. There are no other restrictions on the bounds when the arrays are supplied - see below. If they are arrays, outr can be the same array as inr, and outi can be the same as ini, in which case the original data are overwritten by the transform. Finally, outr and outi may be non-false and not arrays on entry. In this case, recyclable arrays with bounds [0 N-1] are obtained using oldsfloatarray (see *oldarray and *newsfloatarray) and returned. The supplied values are used as the tags in the calls to oldsfloatarray. The point of this is to reduce garbage creation when a procedure needs to carry out repeated transforms. On return, the results outr and outi contain the real and imaginary parts respectively of the transformed data. The computation is carried out in single-precision floating point arithmetic. For best efficiency, all arrays passed to the procedure should be packed arrays of single-precision floats, as created, for example, by *newsfloatarray. If this is not the case, data will be copied as necessary to and from a work array (created using *oldarray to reduce garbage). 2.2 Array subscripting and array bounds ---------------------------------------- To make sensible use of Pop-11 arrays, the following relations hold (apart from machine arithmetic inaccuracies) between array elements and the mathematical quantities for which the transforms are defined, regardless of the array bounds: if dir is , for forward transforms inr(n) = Re(D(n mod N)) outr(k) = Re(F(k mod N)) ini(n) = Im(D(n mod N)) outi(k) = Im(F(k mod N)) else for backward transforms inr(k) = Re(F(k mod N)) outi(n) = Im(D(n mod N)) ini(k) = Im(F(k mod N)) outr(n) = Re(D(n mod N)) Here D and F are assumed to be defined for the range 0..N-1. (Note that these relations are not true in many implementations of the FFT - in particular, a shift of 1 element is common.) In the simplest case, all the arrays have bounds [0 N-1]. The range of the input and output data then corresponds to the formulae in the introduction. Output arrays created by the procedure always have this boundslist. However, this may not always be convenient. The procedure deals with other cases using zero-padding or truncation on inputs, together with the periodicity of both inputs and outputs, as follows. (The procedure used is copy_modn, described in section 4.) Suppose the boundslist of the input arrays is [in0 in1] and that of the output arrays, if supplied by the caller, is [out0 out1]. The number of input data is nin=in1-in0+1 and of output data is nout=out1-out0+1. If nin is greater than or equal to N, N data pairs from the input arrays are used. If in1 is great enough, these will start from the least multiple of N that is not less than in0; otherwise they will start from in0. If nin is less than N, then the input data is padded on the right with zeros, as if the input array actually had bounds [in0 in0+N-1] but with zeros in the elements from in1+1 onwards. (Actually, the padding can equally well be regarded as being on the left, given the cyclic nature of the data.) For a reverse transform, padding will usually only make sense if the boundslist of the input array straddles 0 (or a multiple of N) symmetrically, and truncation will almost never make sense. For reverse transforms, therefore, you should normally expect to have nin=N, and in fact this is advisable where possible for forward transforms too. If in0 is not 0, the input data is then in effect mapped onto the range 0..N-1 using the assumed periodicity of the data, as described above. For example, the case N = 8, in0 = 5, in1 = 10 is handled like this: n: 0 1 2 3 4 5 6 7 8 9 10 11 12 Original data: a b c d e f Data for transform: d e f 0 0 a b c The outputs are actually calculated for k in the range 0..N-1, and for user-supplied output arrays are then mapped onto the range out0..out1 using periodicity as defined above. 2.3 Invertibility ------------------ If a forward transform is followed by a reverse transform, the original data will be reconstructed to within machine precision, provided that: N is the same in both calls; N is greater than or equal to the number of input data; and the output arrays have at least N elements. These conditions are satisfied when N and the output arrays are supplied as and allowed to default. 2.4 The procedure fft_2d ------------------------- fft_2d(N, M, dir, inr, ini, outr, outi) -> (outr, outi) Carries out a 2-dimensional Fast Fourier Transform on complex data. N may be an integer giving the transform length in the first dimension as defined in the introduction. It must be a power of 2. Alternatively, N may be , in which case the transform length is the size of the input array in the first dimension, rounded up to the next power of 2. M is like N, but for the second dimension. dir must be for a forward transform and for a backward transform. inr and ini are 2-D arrays with the same bounds, containg the real and imaginary parts respectively of the input data. There are no other restrictions on the bounds - see above. Both arrays must (at present) be stored "by row" - the normal Pop-11 default. outr and outi may both be input as , or as non-false non-array tags for oldsfloatarray (see fft_1d above), in which case arrays with boundslists [0 N-1 0 M-1] are created and returned. Alternatively, outr and outi may be input as 2-D arrays with the same bounds as one another, in which case the results are stored in these arrays and the same arrays are returned. There are no other restrictions on the bounds when the arrays are supplied - see above. If supplied, the arrays must be stored "by row". If supplied as arrays, outr can be the same array as inr, and outi can be the same as ini, in which case the original data are overwritten by the transform. (If outr shares inr's arrayvector and they overlap, but are not identical arrays, then the offsets and boundslists must be such that outr(k, j) and inr(n, m) refer to the same element of the arrayvector whenever k mod N = n mod N and j mod M = m mod M. A similar restriction applies to outi and ini. This is only likely to be a consideration in the most esoteric applications.) On return, the results outr and outi contain the real and imaginary parts respectively of the transformed data. The computation is carried out in single-precision floating point arithmetic. For best efficiency, all arrays passed to the procedure should be packed arrays of single-precision floats, as created, for example, by *newsfloatarray. If this is not the case, data will be copied as necessary to and from a work array (created using *oldarray to allow subsequent reuse). Note that the term "first dimension" means the dimension relating to the first index of the array (the first argument of the array procedure). In image processing, this usually refers to the column or x-coordinate of a pixel. In matrix algebra, it usually refers to the row of a matrix element. The remarks above about array subscripting, array bounds and invertibility for the 1-D case extend in the obvious way to the 2-D case. Zero-padding or truncation on the inputs, on both dimensions, and tesselation of the plane on the outputs, are used to map between the fundamental region [0 N-1 0 M-1] and whatever boundslists the arrays happen to have. (This is implemented in copy_modn2, described in section 4.) ----------------------------------------------------------------------- 3 Transforms for real arrays ----------------------------------------------------------------------- 3.1 Symmetries for real data ----------------------------- The routines above require both the real and imaginary parts of the data to be specified. This is wasteful when the input data are real, as the resulting symmetry properties of the transformed data can be used to reduce both memory and time requirements. When D is real, the following relation holds (see Section 1 above for details of the notation): F(k) = F*(-k) Here F* means the complex conjugate of F. This is in addition to the periodicity condition F(k) = F(k mod N). Thus in particular we have F(k) = F*(N-k). The corresponding relation in 2-D is F(j, k) = F*(-j, -k) The library has two routines for 1-D transforms of real data, one for each direction, and likewise two routines for 2-D transforms. There are separate routines for forward and backward transforms because the omission of the imaginary part affects the number of arguments and results. Finally there is a special routine for transforming 2-D data that is both symmetric to reflection in the origin and real: D(n, m) = D(-n, -m), D real. In this case the transform F has exactly the same symmetry properties: F(j, k) = F(-j, -k), F real. 3.2 Representing symmetrical data ---------------------------------- In 1-D this is fairly straightforward. For symmetrical data (such as the real part of F when D is real) there are N/2+1 independent values. However, not every span of this length is equally useful, because data on opposite sides of any point congruent to 0 are redundant (or inconsistent). A complete set of data is represented by a span of the form [j j+N/2] where j mod N/2 = 0. Where input data is expected to be symmetrical, the routines will use the leftmost such span in the input region, or if there is no such span, the longest available segment that falls into such a span, with missing data treated as zero. (This is implemented in copy_modn_ref, described in section 4.) Antisymmetric 1-D data, such as the imaginary part of F when D is real, have N/2-1 independent values, because for such data I(i)=0 if i mod N = 0. Values in such elements of input data are ignored. The routines do not use the trick, sometimes found in other packages, of packing the "extra" point R(N/2) into the "spare" point I(0), reducing the array lengths to N/2. In 2-D the situation is more complex. We assume that N and M are both even. The following figure indicating regions of redundancy may be helpful - regions marked with the same letters (ignoring case) determine each other. O is at (0, 0). Taking horizontal to be the first dimension (image convention), N is 10 and M is 12 (for the FFT they would actually be powers of 2). Multiples of N/2 and M/2 are marked in the margins. | | | a a h b b b b d a a a a h b b a a h b b b b d a a a a h b b - c c i c c c c o c c c c i c c - b b h a a a a d B B B B h a a b b h a a a a d B B B B h a a b b h a a a a d B B B B h a a b b h a a a a d B B B B h a a b b h a a a a d B B B B h a a - e e g e e e e F E E E E G e e - a a h b b b b D A A A A H b b a a h b b b b D A A A A H b b a a h b b b b D A A A A H b b a a h b b b b D A A A A H b b a a h b b b b D A A A A H b b - c c i c c c c O C C C C I c c - b b h a a a a d b b b b h a a b b h a a a a d b b b b h a a | | | For symmetrical data, say R, there are MN/2+2 independent values, and a complete set does not form a rectangle. One complete set (shown by capital letters above) is contained in the rectangular region [0 N/2 0 M-1], but within this there is some redundancy because R(0,i)=R(0,M-i) and R(N/2,i)=R(N/2,M-i), as marked with d and h above. Other regions that also contain a complete set are found by shifting the region [0 N/2 0 M-1] any distance in y, and any multiple of N/2 in x. An arbitrary distance in x is not allowed because duplicates will be included. (Here x and y are taken to refer to the first and second dimensions respectively.) Alternatively, it is possible to start from the region [0 N-1 0 M/2], which can shifted arbitrarily in x and by multiples of M/2 in y. Thus suitable input regions are of the form [x0 x0+N/2 y0 y0+M-1], x0 mod N/2 = 0 or [x0 x0+N-1 y0 y0+M/2], y0 mod M/2 = 0 Where input data is expected to be symmetrical, the routines will use such a region if one is contained in the input data. If both kinds of region are available, that of the first type is chosen. If possible the region used will also satisfy y0 mod M = 0 for the first type, or else x0 mod N = 0 for the second type; if this extra condition cannot be satisfied, x0 or y0 will be set to the lower bounds of the input region. If no region of this form is contained in the input, a rectangular region will be chosen to use as many data points as possible, subject to its being contained within one of the "complete data" regions. As many non-redundant points from this region as possible are then used. Missing data is treated as zero. (This is implemented in copy_modn_ref2 described in section 4.) Antisymmetric 2-D data, such as the imaginary part of F when D is real, for which I(j,k) = -I(-j,-k), have MN/2-2 independent values, because the points at I(j,k) for j mod N/2 = 0 and k mod M/2 = 0 are all 0. These are marked as o, i, f and g in the diagram above. 3.3 Routines for 1-D real arrays --------------------------------- fft_1d_real_fwd(N, inr, outr, outi) -> (outr, outi) This procedure carries out a forward transform for real data. The arguments are as for fft_1d, except that ini and dir are omitted. If outr and outi are not arrays, then the arrays returned have bounds [0 N/2]. The symmetry and periodicity conditions are used to fill supplied output arrays regardless of their boundslists. Since outi(0) = outi(N/2) = 0 there are actually only N values computed. fft_1d_real_bckwd(N, inr, ini, outr) -> outr This procedure carries out a backward transform, assuming that the symmetry condition for the transform of real data holds for the inputs. The arguments are as for fft_1d, except that outi and dir are omitted. If N is , it is determined from the number of data points input. N is the smallest power of 2 such that the number of data points is no greater than N/2+1. The comments in the previous section about symmetrical data apply to the inputs. If outr is not an array, the array returned has bounds [0 N-1]. 3.4 Routines for 2-D real arrays --------------------------------- fft_2d_real_fwd(N, M, inr, outr, outi) -> (outr, outi) This procedure carries out a forward transform for 2-D real data. The arguments are as for fft_2d, except that ini and dir are omitted. N and M are determined, if not given, by rounding up the dimensions of inr as for fft_2d. If outr and outi are not arrays, the arrays returned have bounds [0 N/2 0 M-1]. If output arrays are supplied, then they will be filled appropriately, using symmetry and periodicity, so that the correct relations hold for all their elements. fft_2d_real_bckwd(N, M, inr, ini, outr) -> outr This procedure carries out a backward transform for 2-D real data. The arguments are as for fft_2d, except that ini and dir are omitted. inr is assumed to be symmetric by reflection about the origin and ini to be antisymmetric. If N is it is set to be the smallest power of 2 such that the size of the first dimension of inr is no greater than N/2+1. If M is it is set to the size of the second dimension of inr, rounded up to the next power of 2. The comments in the previous section about the region from which data are actually taken apply. If outr is not an array, the array returned has bounds [0 N-1 0 M-1]. 3.5 Routine for 2-D real symmetric arrays ------------------------------------------ fft_2d_real_sym(N, M, dir, inr, outr) -> outr This procedure carries out transforms for 2-D real data that are assumed to be symmetric by reflection through the origin. The arguments are as for fft_2d except that ini and outi are omitted. The forward and backward transforms are essentially the same in this case. The dir argument's main effect is on normalisation: if it is , division by 1/(N*M) takes place, and not otherwise. For this routine, it is probably wise to set N and M explicitly. If they are allowed to default by giving them values, the rules that apply are as follows. Assuming that xsize and ysize are the dimensions of inr, and that N and M are powers of 2, and as small as possible: if dir is true (forward transform) N is no smaller than xsize M/2+1 is no smaller than ysize if dir is false (backward transform) N/2+1 is no smaller than xsize M is no smaller than ysize That is, for forward transforms, inr is expected to contain a region of the second type described in section 3.1 above, and for backward transforms to contain a region of the first type. The reason for this is that the internal structure of the routine makes it most efficient to swap representations on transforming. Note that if N and M are given explicitly, the value of dir has no effect on the region of the input array from which data are taken. inr is assumed to contain data obeying the symmetry and periodicity rules inr(j, k) = inr(-j, -k) and inr(j, k) = inr(j mod N, K mod M). Inconsistent data are ignored - a consistent set is taken from the largest available region, as described in section 3.1. For highest efficiency, the boundslist of inr should be [0 N/2 0 M-1]. However, a consistent approach is to use [0 N-1 0 M/2] for forward transforms and [0 N/2 0 M-1] for reverse transforms. On exit, outr receives data which satisfy the symmetry and periodicity requirements. If outr is supplied, highest efficiency is obtained if its boundslist is [0 N-1 0 M/2], but a consistent approach is to adopt the default output boundslists that depend on dir, as follows. If outr is not supplied, the array returned has boundslist [0 N/2 0 M-1] if dir is , and [0 N-1 0 M/2] if dir is . The arrays inr and outr can be the same, in which case the original data are overwritten. This routine makes considerable use of internal work arrays. These are created using oldanyarray (see HELP * OLDARRAY) and so are garbage-efficient. ----------------------------------------------------------------------- 4 Convenience routines ----------------------------------------------------------------------- A few routines which are mainly for internal use in the library are made global as they may be generally useful. nextpow2(int1) -> int2 Returns the smallest power of 2 which is not less than int1. copy_modn(s0, s1, soff, src, d0, d1, doff, dst, N) src and dst must be vectors (of any type); the other arguments must all be positive integers. src and dst may be the same vector. If soff and doff are zero, elements d0 to d1 of dst are filled from elements s0 to s1 of src, such that dst(i) = src(j) if i mod N = j mod N and j is one of the elements actually accessed. See the discussion in section 2.2 which explains what elements are actually used and how. The result is always periodic with period N. soff and doff are offsets, and all vector references are actually relative to src(soff) and dst(doff). That is, they give the indices in the real vector which correspond to the zero points of the mathematical vectors for which the relation above is defined. 1-D arrays can be passed to this routine as the arrayvector and the array offset in the arrayvector (that is min_sub result of *arrayvector_bounds minus the first element of the boundslist). copy_modn_ref(s0, s1, soff, src, d0, d1, doff, dst, N, revsgn) As copy_modn but incorporating reflection symmetry. See the discussion in section 3.1 which explains what elements are actually used and how. If revsgn is not , those segments which are reversed in direction on copying are negated, to give antisymmetry. copy_modn2(sreg, sarr, dreg, darr, N, M) 2-D equivalent of copy_modn. sreg and dreg must be 4-element lists giving regions of sarr and darr from which to take and into which to place data. If they default to the boundslist of the arrays. sarr and darr may be the same array. sarr and darr must be 2-D arrays, of any type. Data are copied from the region sreg of sarr to the region dreg of darr, respecting periodicity of N in the first dimension and M in the second. See sections 2.2 and 2.4 for what data are actually used, and how. copy_modn_ref2(sreg, sarr, dreg, darr, N, M, revsgn) 2-D equivalent of copy_modn_ref. Arguments are as for copy_modn2, with the addition of revsgn which if non- causes data in regions which have been reflected to be negated, giving antisymmetry. See section 3.1 for a full description of the fairly complex behaviour of this routine. ----------------------------------------------------------------------- 5 Examples ----------------------------------------------------------------------- Examples below should all work if executed in order. First load some libraries, uses popvision, fft, showarray, rc_graphplot, rci_show, arrayset; 5.1 Examples for 1-D data -------------------------- Define a handy graph plotting routine: define graph(arr); dlocal rc_window; rci_show([1 300 1 250]) -> rc_window; lvars (i0, i1) = explode(boundslist(arr)); rc_graphplot(i0, 1, i1, 'i/j', arr, 'D/F') -> enddefine; 10 -> rci_show_xshift; 35 -> rci_show_yshift; ;;; space out windows Create 1-D sinusoidal data, with 230 points, starting from D(1094): vars b = [1094 ^(1094+229)]; ;;; boundslist of data vars d1r = newsfloatarray(b, procedure(x); sin(10*360*x/200) endprocedure); vars d1i = newsfloatarray(b); ;;; zero complex part Create a window (click on this and any subsequent graphics windows if you want to get rid of them) and look at the input: graph(d1r); Take the FFT of the data using the general purpose routine and plot the results: vars (f1r, f1i) = fft_1d(false, true, d1r, d1i, false, false); graph(f1r); graph(f1i); Note the symmetries because the data were real. The length of the transform got rounded up to 256 (the next power of 2 above 230) so that is the length of the output arrays. The wiggles either side of the peaks are causes by the finite length of the sine wave input. Transform back: vars (d1r_n, d1i_n) = fft_1d(false, false, f1r, f1i, false, false); graph(d1r_n); graph(d1i_n); Note that the complex part is close to zero, but the graph shows numerical noise. The real part is now on the range 0 to 255, the default, and the zero padding which was inserted for the forward transform is revealed. You should be able to see why it appears where it does from the discussion in section 2.2, bearing in mind that (1094+230) mod 256 is 44, and 1094 mod 256 is 70. To transform back to the original span, one could do: vars (d1r_n, d1i_n) = fft_1d(false, false, f1r, f1i, copy(d1r), copy(d1r)); graph(d1r_n); By passing in a new array with the same bounds as d1r, we get back something that is almost identical with the original input, and the zero padding vanishes again. Since the data are real, all this can be done at lower cost using the transforms for real data: rci_show_destroy("all"); graph(d1r); vars (f1r, f1i) = fft_1d_real_fwd(false, d1r, false, false); graph(f1r); graph(f1i); vars d1r_n = fft_1d_real_bckwd(false, f1r, f1i, copy(d1r)); graph(d1r_n); Note that f1r and f1i now have half the length they did before - what is missing was redundant, and the original data are still recreated almost exactly. For a simple example of frequency domain filtering, we replace the input with a square wave: vars d1r = newsfloatarray(b, procedure(x); 2 * (x div 20) mod 2 - 1; endprocedure); graph(d1r); vars (f1r, f1i) = fft_1d_real_fwd(false, d1r, false, false); Now set the high frequencies to zero and transform back: vars r, i; for r, i in_array f1r, f1i updating_last 2 in_region [16 128] do 0 ->> r -> i endfor; vars d1r_n = fft_1d_real_bckwd(false, f1r, f1i, copy(d1r)); graph(d1r_n); Remove all the windows if necessary: rci_show_destroy("all"); 5.2 Examples for 2-D data -------------------------- Make a 2-D checkerboard test pattern, and a zero complex part: vars b = [30 190 20 170]; vars d2r = newsfloatarray(b, procedure(x, y); (sin(9*x) * sin(14*y)) ** 2 - 0.5 endprocedure); rci_show(d2r) -> ; vars d2i = newsfloatarray(b); Transform with the standard fft: vars (f2r, f2i) = fft_2d(false, false, true, d2r, d2i, false, false); rci_show(f2r) -> ; rci_show(f2i) -> ; The transforms may not appear to have much information, but if you look closely you should be able to see the bright pixels representing the spikes in the output. Transform back with all the defaults: vars (d2r_n, d2i_n) = fft_2d(false, false, false, f2r, f2i, false, false); rci_show(d2r_n) -> ; The output reconstructs the input except that we now see the zero-padding that was used to make the dimensions a power of 2 on the forward transform. We could lose this by passing in an array the same size as d2r to accept the result. The same results can be obtained by using the versions of the arrays for real data - the transform arrays will then each be half the size, and the transforms will be faster. We can do simple frequency-domain filtering in the 2-D case also. In this case we will lose the low frequencies. It will be simpler if the transform arrays are centred on zero, so we will redo the transform first, giving it centred arrays to work on. We supply N and M explicitly now as otherwise the reverse routine will get the wrong defaults. Note how the spikes are now all close to the origin when the arrays are displayed. vars tb = [-128 128 -128 128]; vars f2r = newsfloatarray(tb); vars f2i = newsfloatarray(tb); fft_2d(256, 256, true, d2r, d2i, f2r, f2i) -> (f2r, f2i); rci_show(f2r) -> ; rci_show(f2i) -> ; (We could also have recentred the previous outputs using copy_modn2. This would have been faster, but the example above starts from the original data to illustrate how fft_2d can produce outputs in any required region.) Now get rid of the low frequencies, transform back and display: arrayset(0, f2r, [-20 20 -20 20]); arrayset(0, f2i, [-20 20 -20 20]); vars (d2r_n, d2i_n) = fft_2d(256, 256, false, f2r, f2i, false, false); rci_show(d2r_n) -> ; Note how the boundary of the data region appears as a strong high-frequency edge. Again, this all works just the same, though about twice as fast, with the real array routines, because we have been using real data. fft_2d_real_fwd(256, 256, d2r, f2r, f2i) -> (f2r, f2i); arrayset(0, f2r, [-20 20 -20 20]); arrayset(0, f2i, [-20 20 -20 20]); fft_2d_real_bckwd(256, 256, f2r, f2i, false) -> d2r_n; rci_show(d2r_n) -> ; We did not have to alter the calls to arrayset because the routines take care of filling the transform arrays properly. However, we could have made them smaller had we wished to and achieved the same results by only dealing with values in one half-plane. This is left as an exercise. --- $popvision/help/fft --- Copyright University of Sussex 2001. All rights reserved.