25. Fourier Transforms¶
PDF pages 1133–1148
Chapter 25 Fourier Transforms Fourier transforms are the basis for a number of powerful analytical tools for analyzing linear (and sometimes nonlinear) systems. But they are also important numerical tools, in part because they can be performed accurately and very efficiently. Knowing how to use them both analytically and numerically is a good way to get deep intuition into physical systems, in particular dynamical systems, and allows for the construction for sophisticated numerical analysis techniques. There are numerous conventions for the Fourier transform, and we will discuss a few of them. For the sake of concreteness, we will consider the usual time-frequency convention for the Fourier transform F of a function f(t)
Z \infty −\infty f(t) ei\omegatdt (25.1) (Fourier-transform definition) and the inverse Fourier transform F −1 h ˜f i
2\pi Z \infty −\infty
(inverse-Fourier-transform definition) (25.2) We will then discuss adaptations to other normalization conventions. 25.1 Sampling Theorem A critical aspect of the numerical computation of a Fourier transform is adapting the integral transform to a finite sample of the data. The sampling theorem provides the basis for doing this, especially for understanding the errors involved in doing so. To understand the sampling theorem, suppose that the spectrum of f(t) has compact support—that is, suppose that ˜f(\omega) vanishes for |\omega| > \omegac, where \omegac is some ‘‘cut-off frequency.’’ Then defining the unit-rectangular-pulse function by
if |t| < 1/2 1/2 if
if |t| > 1/2, (25.3) we can write the compact-support condition for the spectrum as
\omega 2\omegac . (25.4) Recall (Section 17.1.2) that the convolution theorem for functions f(t) and g(t) reads F[f ∗g] = F[f]F[g], (25.5)
25.1.1 Critical Sampling¶
Chapter 25. Fourier Transforms where ‘‘∗’’ denotes the convolution operation
Z \infty −\infty f(t′)g(t −t′) dt′. (25.6) We will also later need the frequency-domain version, which reads F −1 h ˜f ∗˜g i
˜f i F −1 [˜g] (25.7) for frequency-domain functions ˜f(\omega) and ˜g(\omega), where the extra factor of 2\pi is due to the same factor in the inverse transform (25.2). Applying the original form (25.5) of the convolution theorem to the compact- support condition (25.4), we may write
\pi sinc \omegact ∗f(t), (25.8) (prelude to sampling theorem) where sinc x := sin x/x, and we have used the inverse Fourier transform
(25.9) Thus, we see that f(t) is an invariant under convolution with the sinc function. 25.1.1 Critical Sampling Now suppose we discretely sample the function f(t) at uniform time intervals ∆t. That is, we represent f(t) by the countable set of values
(25.10)
weighted comb function
\infty X j=−\infty fj ∆t \delta(t −j∆t). (25.11) This is a function that (1) is determined only by the samples fj, and (2) has the same coarse-grained area as f(t), at least to O(∆t2). To compute the Fourier transform of this function, note that we may rewrite it as the product of f(t) and the usual comb function:
\infty X j=−\infty \delta(t −j∆t). (25.12) Now using the fact that the Fourier transform of a comb function is a comb function, F \infty X j=−\infty \delta(t −j∆t)
∆t \infty X j=−\infty \delta
∆t , (25.13) we can use the convolution theorem in the form (25.7) to write the Fourier transform of f (∆t)(t) as the convolution of ˜f(\omega) with a comb function:
\infty X j=−\infty \delta
∆t . (25.14) This spectrum is periodic in frequency due to convolution with the comb function. Now suppose we choose the sampling interval ∆t such that it is determined by the cut-off frequency by
∆t. (25.15)
25.1.2 Reconstruction¶
25.1 Sampling Theorem Then¶
\infty X j=−\infty
, (25.16) and we see that now the spacing between the ‘‘teeth’’ of the comb function is 2\omegac, the same as the width of the (two-sided) spectrum ˜f(\omega). This means that, provided we are interested in the range |\omega| < \omegac, only one
h ˜f ∗\delta i
(25.17) Another way to state this is that for any frequency, \omega, the spectra satisfy ˜f (∆t)(\omega) rect \omega 2\omegac
\omega 2\omegac
(equivalence of sample-reconstructed and original spectra) (25.18) Thus, we see that the Fourier transforms of f(t) and f (∆t)(t) are the same. In other words, the samples fj, where the sampling interval satisfies the condition (25.15) for critical sampling of the signal, are sufficient to completely reconstruct the spectrum ˜f(\omega), and thus the original function f(t). This is essentially the content of the sampling theorem, without going into the precise conditions for its validity (as is usual in a physicist’s treatment of the subject). Of course, what this means is that a function whose spectrum has compact support is in some sense a very special function. However, so long as the weight of an arbitrary spectrum is very small outside the range (−\omegac, \omegac), the error in the reconstruction of the spectrum is correspondingly very small. The other situation in which this works is if the spectrum ˜f(\omega) is periodic in \omega, since we can see from Eq. (25.16) that the reconstructed spectrum is periodic with period 2\omegac. This makes sense, as we know from Fourier theory that a function that is either defined on a bounded interval or is periodic can be represented by a Fourier series, rather than a Fourier transform. Thus, in writing down Eq. (25.11), we have essentially written down the Fourier series for ˜f(\omega), but in the language of a Fourier transform. 25.1.2 Reconstruction The sampling theorem also provides a direct way to ‘‘reconstruct’’ the original function from its samples. We assume the spectrum has compact support—otherwise, a periodic spectrum implies that the modulated comb (25.11) is in fact the true form of f(t)—and then combine the compact-support condition (25.4) with the reconstructed-spectrum condition (25.18) to write
\omega 2\omegac . (25.19) Using the convolution theorem (25.5), we can thus write the inverse transform of this equation as
h\omegac
i , (25.20) where we have again used the inverse Fourier transform (25.9). Writing this out in terms of the samples,
\infty X j=−\infty fj∆t\omegac
(25.21)
Chapter 25. Fourier Transforms Again using the critical-sampling condition (25.15) and carrying out the convolution, this relation becomes the Whittaker–Shannon sampling interpolation formula:1
\infty X j=−\infty
\infty X j=−\infty fj sinc [\pi(t −tj)/∆t] (Whittaker–Shannon sampling interpolation formula) (25.22) The sinc functions are used to interpolate the function between the samples, and again, if the signal is bandwidth-limited and critically sampled, this reconstruction formula is exact. Note that the function
the sample times. What is less obvious is the exactness of the formula between the sample times. 25.2 Discrete Fourier Transform Any real-world numerical calculation of a Fourier transform will obviously require operations on finite data sets. By sampling the temporal signal, we have reduced the information from a function f(t) on an un- countable set to a function fj on a countable (discrete) set, but we must further reduce the information to a finite set. Of course, since sampling implied a spectrum of finite width, we can obtain a finite set simply by
frequency intervals ∆\omega. This is equivalent to a truncation of the time samples, so that the time samples only occur within some frequency interval. Again, this will amount to the assumption that the temporal signal is either a pulse with compact support, or that it is periodic. When time and frequency are both sampled, we can use the arguments above to impose constraints on the sample intervals and ranges. For example, as above, when the signal is temporally sampled with N points with interval ∆t, we may regard the signal as extending from t = 0 to t = 2tmax, where tmax = N 2 ∆t. (25.23) Note the factor of 2 here, since the sampling in both time and frequency imply the assumption that the
t = tmax, which is why we have set up our notation this way. The sampling interval ∆t, from our arguments in the last section, implies a maximum frequency
∆t, (25.24) which is called the Nyquist frequency, or the largest frequency that is critically sampled. Thus, the spectrum extends in frequency from \omega = −\omegamax to \omega = \omegamax. To have the same information, there will also be N samples in frequency, N/2 of which correspond to the range from \omega = 0 to \omegamax, so that the frequency-sampling interval is
N = 2\pi N∆t = \pi tmax . (25.25) Thus, given the two free parameters, the total time 2tmax for the sample and the total number of samples N, the above three formulae give the rest of the discretization parameters ∆t, ∆\omega, and \omegamax. To avoid all these details of time and frequency scales, the discrete Fourier transform (DFT) is defined to have the simple dimensionless form Fk = N−1 X j=0 fje2\piijk/N, (25.26) (discrete Fourier transform) 1The history of the sampling theorem is somewhat involved and unclear (convoluted?). For a nice overview see the Wikipedia entry, ‘‘Nyquist–Shannon sampling theorem,’’ http://en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem.
25.2 Discrete Fourier Transform and essentially amounts to multiplying the fj vector by the matrix exp(2\piijk/N). The inverse of this transform will obviously have the opposite sign in the exponential: ˆfj = N−1 X k=0 Fke−2\piijk/N. (discrete inverse Fourier transform, unscaled) (25.27) However, this transformation does not exactly invert the DFT, because the inverted numbers are a factor of N larger than the originals, ˆfj = Nfj. (25.28) This follows from the identity N−1 X k=0
(25.29) which is obvious for j = j′, and in the case j̸ = j′ can be seen to vanish because the left-hand side amounts to a sum over all Nth roots of unity. Thus, the inverse DFT is often defined with this factor of N: fj = 1 N N−1 X k=0 Fke−2\piijk/N. (discrete inverse Fourier transform, scaled) (25.30) Unfortunately, there is no standard convention over which inverse, (25.27) or (25.30) is implemented as the ‘‘inverse Fourier transform.’’ The unnormalized form (25.27) is commonly implented in low-level languages such as Fortran and C (a sensible choice since the user usually needs to multiply by extra normalization factors anyway), while the ‘‘normalized’’ inverse transform (25.30) is implemented in MATLAB and Octave. For any particular implementation of the DFT, you should consult the documentation, or just test the DFT, followed by an inverse DFT, on an array of ones to see if you get back the array with an extra factor of N, or if you get back the original array. To map our desired integral Fourier transforms onto the above DFTs, we start by writing down the truncated form for the approximation (25.12) for f(t), corresponding to N samples:
N−1 X j=0 \delta(t −j∆t). (25.31) The sum could equally well run over positive and negative times, since the function is effectively assumed to be periodic anyway. Now putting this into the Fourier-transform formula (25.1), we find
Z \infty −\infty f (∆t,N)(t)ei\omegatdt = Z \infty −\infty f(t) ∆t N−1 X j=0 \delta(t −tj)ei\omegatdt = N−1 X j=0 f(tj)ei\omegatj∆t. (25.32) The last expression is the usual discrete approximation for an integral, applied to the Fourier integral. The approximation formula usually has O(∆t2) error, but in view of the sampling theorem, the error is really due just to the truncation, and can thus be much smaller than O(∆t2), as we will discuss below. In particular, at the N frequencies \omegak = k∆\omega,
N−1 X j=0 fjei\omegaktj∆t. (25.33)
25.2.1 Periodicity and Transform Ordering¶
Chapter 25. Fourier Transforms
N−1 X j=0 fje2\piijk/N. (25.34) Thus, apart from the factor of ∆t, we have exactly that the DFT formula (25.26) approximates the Fourier integral, in the sense that ˜fk = Fk ∆t. (25.35) (DFT as approximation to Fourier integral) so that we must multiply by ∆t after applying the DFT to obtain the approximation to the Fourier integral. Similarly, if we adapt the inverse transform (without the factor of N) to the spectrum samples ˜fk, we find
2\pi N−1 X k=0 ˜fke−2\piijk/N. (25.36) That is, putting in the ˜fk for the Fk in the inverse DFT formula (25.27), the actual inverse transform samples are given by multiplying the results of the inverse DFT by ∆\omega: fj = ˆfj ∆\omega 2\pi . (inverse DFT as approximation to Fourier integral) (25.37) Of course, the factor of 2\pi here is the same factor in the inverse Fourier integral (25.2). 25.2.1 Periodicity and Transform Ordering One thing that should be clear from the above expressions is that the zero-frequency and zero-time com- ponents ( ˜f0 and f0, respectively) occur at the beginning of their respective arrays. This may seem odd, as the frequency spectrum has both positive and negative frequencies. Of course, the frequency spectrum is assumed by the DFT to be periodic with period 2\omegamax, which corresponds to a shift by N in the index. Thus, rather than thinking of the ordering ˜f0, ˜f1, . . . , ˜fN−1 (25.38)
we may regard the elements in the same order as
(25.39) Here, we have implicitly assumed that N is even, which is almost always the case in practice. It is often convenient to have the zero frequency actually in the center of the array, in which case we simply swap the left and right halves of the array to obtain
(25.40) Note the asymmetry here, since f−N/2 appears at the left boundary, while fN/2−1 appears at the right. However, the N-periodic nature of the discrete spectrum means that ˜fN/2 = ˜f−N/2, so in principle we can copy the left boundary onto the right boundary. The same ordering comments apply to the temporal signal. A straight signal out of your digital sampling oscilloscope would have the form f0, f1, . . . , fN−1 (25.41)
a spectrum of the form (25.38). On the other hand, you might have something like a correlation function
25.3 Aliasing¶
the negative-time values. This is also easy to do, since the temporal array is also N-periodic (with period
(25.42) For the correlation function, then we would impose the constraint f−N = f ∗ N. You can also swap the left and right halves to get a properly time-ordered array after all the DFT operations are done, and again it is still true that fN/2 = f−N/2. 25.3 Aliasing As we can see from the sampling theorem, the only error that we are introducing here is the fact that we are assuming f(t) and ˜f(\omega) both have compact support (or equivalently, are periodic). If both had compact support, the DFT would be exact (at least up to rounding errors, if the DFT is performed on a computer). However, very few functions have compact support in both the time and frequency domain: remember that having compact support in one domain is equivalent to convolution with a sinc function in the other, and sinc x only falls off like 1/x. Thus, the set of all functions for which the DFT can be exact is the set of all functions on a compact domain whose Fourier series are eventually zero (i.e., truncated). This set is obviously of zero measure in the space of all functions that we would want to transform. So, then, what is the error associated with the DFT? For a general time signal f(t) sampled with time interval ∆t, in general the true spectrum ˜f(\omega) will have some contributions beyond the Nyquist frequency \omegamax. Since the frequency spectrum is assumed by the DFT to be periodic with period 2\omegamax, the parts beyond \omegamax are spuriously ‘‘folded’’ into the computed spectrum. This effect is illustrated here for a Gaussian spectrum centered at \omega = 0. max -ow w max w Gaussian aliased spectrum Another way to see why this must be is that the DFT (25.26) is a unitary transformation on N points (except for the factor N −1/2), and thus the total power of the sampled time signal must be equivalent to the total power of the sampled frequency spectrum—the discrete form of Parseval’s theorem—whether or not the spectrum ‘‘fits’’ within the range \pm\omegamax. Yet another good way to visualize this is to look at a harmonic wave that is not critically sampled. Recall that a frequency is critically sampled if it is sampled at least twice per period. However, if the sampling is below the critical rate, the samples will be indistinguishable from those from a lower-frequency harmonic wave (not to mention an infinity of yet higher-frequency harmonic waves). This problem even occurs in the time domain in the laboratory on digital sampling oscilloscopes: if you’re measuring a fast sin wave on such a scope, and you’re seeing a wave with a frequency way lower than you expect, try cranking the ‘‘time/div’’ knob to see if you get a more appropriate (non-aliased) signal on a faster time scale. (Incidentally, this is the same effect that makes vibrating objects appear to be stationary or oscillate slowly with a strobe light, or that makes spinning wheels, propellers, and helicopter blades appear to precess slowly or even backwards on film or television.2) Of course, the same comments apply to starting 2For a good example, see http://www.youtube.com/watch?v=eJ6vadFVjYg.
Chapter 25. Fourier Transforms off with a frequency spectrum, and then using the inverse DFT to obtain a time signal (say, a correlation function): aliasing can still occur in the time domain. The point of all this is, to get accurate results with the DFT, say, when you are computing the spectrum of a pulse, you must pick your time grid such that the pulse and the spectrum fit well within 2tmax and 2\omegamax, respectively. That is, near the outer boundaries of the \omega and t grids, the signals in both domains should have fallen very close to zero, so you can neglect the power that falls outside the boundaries. Sometimes, you want to compute the DFT of a signal that doesn’t fall off, such as a stationary, fluctuating signal. In this case you should still choose an appropriate sampling rate to avoid spectral aliasing, and realize that there may be some artifacts if the signal is not exactly periodic (as when sampling a sin wave, when the length of the sample does not match the period of the wave). The artifacts should be minor, however, so long as the sampling time is long compared to the correlation time of the signal. 25.4 Fast Fourier Transform A main reason why the DFT is such an important computational tool is that it can be done so efficiently. From the DFT formulae (25.26) and (25.27), the DFT can be viewed as a multiplication of a matrix and a vector, and so for an array of length N, the computational effort (number of multiplications) should be O(N 2). However, using a class of algorithms called Fast Fourier Transforms (FFTs) can do the same calculation in only O(N log N) operations. By far the most common algorithms are specific to the case where N is a power of 2, in which case the operation count is O(N log2 N). Note that it is generally best to stick with these power-of-2 algorithms, and just live with this constraint on N: your data arrays can generally be interpolated or padded with zeros to get an appropriate array length. This savings in computational effort
points takes about 0.3 s. Using the scalings above, a direct implementation of the DFT formula would take on the order of 2.4 hours for the same calculation! You might think that a couple of hours isn’t so bad, but it is common for many FFTs to be needed in a single calculation. Also, because of the reduced operations count, round-off error does not accumulate nearly as much as in the straight DFT, and FFT algorithms are generally highly accurate. The FFT algorithms are so important that ‘‘FFT’’ is commonly used in place of ‘‘DFT.’’ The first FFT algorithm dates back to Gauss in 1805, but was not widely known until it was redis- covered by Cooley and Tukey3 in 1965 (the Cooley-Tukey algorithm is a recursive method and remains a popular algorithm).4 Many other algorithms beyond the Cooley-Tukey method are possible. We will not go into the details of various FFT algorithms,5 since different implementations have different advantages, especially regarding accuracy and execution time on different architectures.6 In this case, it is generally best to stick to algorithms written by specialists for reliability and good performance. Compiler vendors often provide highly optimized versions, but otherwise FFTPACK7 and FFTW8 are well-known libraries. 3James W. Cooley, John W. Tukey, ‘‘An Algorithm for the Machine Calculation of Complex Fourier Series,’’ Mathematics of Computation 19, 297 (1965). 4For a nice overview of FFT algorithms, accuracy, and history, see the Wikipedia entry ‘‘Fast Fourier transform,’’ http://en.wikipedia.org/wiki/Fast_Fourier_transform. Another good source, particularly for history and algorithm de- tails is the Wikipedia entry ‘‘Cooley-Tukey FFT algorithm,’’ http://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm. 5For the basic idea, see William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed. (Cambridge, 1992), Section 12.2, p. 498 (doi: 10.2307/2003354). 6For an example of how much a careful choice of method can make a big difference in execution time, see David H. Bailey, ‘‘A High-Performance FFT Algorithm for Vector Supercomputers,’’ International Journal of Supercomputer Applications 2, 82 (1988) (doi: 10.1177/109434208800200106). An algorithm that performs reasonably well with large arrays on modern cache- based computers is the Stockham FFT (implemented in FFTPACK), detailed in Paul N. Swarztrauber, ‘‘FFT algorithms for vector computers,’’ Parallel Computers 1, 45 (1984). For comparisons among many algorithms, see the benchFFT home page, http://www.fftw.org/benchfft/. 7http://www.netlib.org/fftpack/ 8http://www.fftw.org/
25.5.1 Temporal Signals¶
25.5 Conventions 25.5 Conventions We have seen above how to use the general DFT formulae to construct numerical approximations to Fourier integral transforms. But keeping track of all the different sign, normalization, and ordering conventions can be a pain, so we go through a few examples here, considering bits of pseudocode appropriate for MAT- LAB/Octave and Fortran 90/95. For MATLAB/Octave we will consider the built-in functions fft and ifft, which implement the formulae (25.26) and (25.30), respectively. In Fortran, we will use a fictitious (but typ- ical) subroutine fft(isign, data), where data is the data array (both input and output), and isign
since the generalization to higher dimensions is reasonably obvious. 25.5.1 Temporal Signals The simplest case we will consider is computing the energy spectrum | ˜f(\omega)|2 of a temporal signal f(t), or the power spectrum | ˜f(\omega)|2/T, where T is the duration of the sample. The DFT formula (25.26) is easily adapted to this case for computing ˜f(\omega), since we simply need to multiply by ∆t to obtain the correct scaling:
N−1 X j=0 fje2\piijk/N ∆t. (25.43) The resulting spectrum array will have its zero frequency component first, with the negative frequencies in the last half of the array. The negative frequencies can be discarded in the power spectrum, since they will just be a mirror reflection of the positive frequencies. In the MATLAB/Octave code below, we simply compute the Fourier transform of the signal array sig of length N, and multiply by dt. The result is stored in the temporary array spec (also of length N), and after the negative frequencies are discarded, the result is squared and scaled, and is then stored in energyspec or powerspec, which is of length N/2. The scaling factor is N*dt, which is the total time T of the sample, and we cast the vector scaling as a multiplication to avoid many slower divisions. A frequency grid w is also generated for illustration, using Eqs. (25.23)-(25.25), so that the resulting array could be plotted with plot(w, pwrspectrunc, '-'). % dt is the time sampling interval
powerspec
% w is the vector of frequencies wmax = pi/dt;
Below is the equivalent code in Fortran 90/95. Note that the variable declarations are not included, but some things to note: wp is a parameter declaring the floating-point precision (e.g., declare as
! dt is the time sampling interval
call fft( 1, spec)
powerspec(1:N/2)
! w is the vector of frequencies
25.5.2 Temporal Correlation Functions¶
Chapter 25. Fourier Transforms wmax = pi/dt
Typically then the values w and pwrspectrunc would then be dumped to a file or standard output for further processing. In either case, to be accurate and avoid artifacts, the length of the sample of f(t) should be well beyond the correlation time. This could be checked, for example, by computing the correlation function from the power spectrum (see below) and verifying that it decays to zero before the array ends. 25.5.2 Temporal Correlation Functions
want to compute the corresponding power spectrum. Your data array is of length N. Before starting, to minimize artifacts it is useful to ‘‘enforce’’ the periodicity by pasting the complex-conjugated mirror image of the correlation function onto the end of the array, to obtain a ‘‘periodic’’ array of length 2N. (There should be no discontinuity at the pasting point in the middle of the new array, because the samples should
sample, since we don’t want to duplicate it, and this requires pasting in an extra zero. That’s because the
Again, we can adapt the DFT formula (25.26) for this purpose,
2N−1 X j=0 fje2\piijk/2N ∆\tau, (25.44) so that again after computing the DFT we simply multiply by ∆\tau. Note that we have 2N in place of the usual N, since we have doubled the array length before the transform. Also, to center the dc component in the middle of the array, we swap the two halves of the transform with the fftshift function. To recover the correlation function from the power spectrum samples ˜gk, we adapt the inverse DFT formula (25.30) gj = " 2N 2N−1 X k=0 ˜gke−2\piijk/2N # 2N∆\omega 2\pi . (25.45) The difference here is that we need to multiply by ∆\omega instead of ∆\tau for the frequency integral, undo the factor of 2N in the inverse DFT, and divide by 2\pi for the time-frequency normalization convention of Eq. (25.2). This is the proper approach in MATLAB/Octave, where the ifft function implements the bracketed transform in Eq. (25.45). % dtau is the time sampling interval
% w is the vector of frequencies wmax = pi/dtau;
% to recover g from pwrspec
In Fortran 90/95, though, typical implementations leave out the factor of 2N, so we are performing an inverse DFT as in Eq. (25.27): gj = "2N−1 X k=0 ˜gke−2\piijk/2N # ∆\omega 2\pi . (25.46)
25.5.4 Wave Functions¶
25.5 Conventions Thus, we need not multiply by 2N after the DFT. Also, note that swapping the two halves of the frequency array is conveniently performed with the Fortran intrinsic cshift (cyclic array shift, specifically for half the array length). ! dtau is the time sampling interval pwrspec = 0
call fft( 1, pwrspec)
! w is the vector of frequencies
wmax = pi/dtau
! to recover g from pwrspec (destroys pwrspec in the process)
call fft(-1, pwrspec)
routine especially adapted for real inputs and outputs. 25.5.3 Standard Frequency Often in DFT applications, we will want to deal with the standard frequency \nu instead of the angular frequency \omega = 2\pi\nu, in which case the alternate Fourier transform becomes
Z \infty −\infty f(t)ei2\pi\nutdt, (25.47) and the inverse Fourier transform becomes
Z \infty −\infty
(25.48) That is, there is no longer the factor of 1/2\pi scaling the inverse transform, and now there are explicit factors of 2\pi in the exponents. Everything is the same as in the \omega–t convention, except the Nyquist frequency is now
2\pi = 2∆t, (25.49) and the frequency-sampling interval is
\piN = N∆t = 2tmax . (25.50) Other than these changes in the frequency grid, the only other difference in the above code snippets is that the division by 2\pi after the inverse DFT should be omitted. 25.5.4 Wave Functions To transform a wave function between position and momentum space—as is useful, for example, in impl- menting split-operator methods for evolving the Schrödinger equation, as in Chapter 26—the conventions are a bit different than for time and frequency. To compute the Fourier transform of \psi(x) to obtain the momentum-space version \phi(p), the integral is
\sqrt 2\pi¯h Z \infty −\infty \psi(x)e−ipx/¯hdx, (25.51)
Chapter 25. Fourier Transforms while the inverse Fourier transform is
\sqrt 2\pi¯h Z \infty −\infty \phi(p)eipx/¯hdp. (25.52) Note the symmetric normalization and the presence of the extra ¯h. Furthermore, the x = 0 and p = 0 components are generally at the center of the wave-function arrays, so the array halves must be swapped both before and after the transform. We can adapt the inverse DFT (25.30) for the Fourier transform here, because of the opposite sign convention of the phase factors,
1 N N−1 X j=0
N∆x \sqrt 2\pi¯h , (25.53) as well as the DFT (25.26) for the inverse Fourier transform,
"N−1 X k=0
∆p \sqrt 2\pi¯h . (25.54) Also, since p/¯h plays the role of frequency, the relations between the increments are as follows, if we take the number N of grid points and the extent of the grid from −xmax to xmax to be fixed: ∆x = 2xmax N¶
∆x ∆p = 2pmax N
N∆x. (25.55) Thus, the MATLAB/Octave code would read as follows: % specify N and xmax for the x grid % x is the x grid
% p is the vector of momenta pmax = hbar*pi/dx;
% to recover psi from phi
In Fortran 90/95, we again simply specify the sign of the exponent for the transform, and forget about extra factors of N. ! specify N and xmax for the x grid ! x is the x grid
call fft(-1, phi)
! p is the vector of momenta
pmax = hbar*pi/dx;
25.6 Discrete Wigner Transform¶
! to recover psi from phi
call fft( 1, psi)
In this code, psi and phi are obviously complex arrays of length N. 25.6 Discrete Wigner Transform Recall from Section 4.3 that the Wigner transform of a wave function \psi(x) is
2\pi¯h Z \infty −\infty
(25.56) We can simply regard W(x, p) as a (quantum) Fourier transform of the form
\sqrt 2\pi¯h Z \infty −\infty
(25.57) where the function to be transformed in
\sqrt 2\pi¯h
(25.58) However, because of the appearance of x′/2 in the argument of \psi, if \psi(x) is sampled with N samples with an interval of ∆x, the appropriate sampling interval to use for the Fourier transform is ∆\xi = n∆x, where
always pad \psi with zeros such that f(\xi) is always defined from −nxmax to nxmax, and thus that f(\xi) is still represented by N samples. Then we have the following modified parameters for the grid for W(x, p): ∆x = 2xmax N
¯h\pi n∆x ∆p = 2pmax N
2\pi¯h Nn∆x. (25.59) The Fourier transform must be repeated for each point in the x grid from −xmax to xmax −∆x in steps of ∆x. Of course, you can skip some of these position values, as when making a three-dimensional plot, it is generally best to keep the density of points the same in both the position and momentum directions. In any case, we can write the explicit formula for the discrete Wigner transform as
2\pi¯h N/2−1 X
(25.60)
the even integer n \ge 2. Note that with this ordering of l, the zero-momentum component is in the center of the array, so array-swapping operations as in Section 25.5.4 are necessary to use the DFT/FFT formulae to evaluate the summations here. For the discrete Wigner transform, the MATLAB/Octave code would read as follows (NN corresponds to N above): n = 4; % must be even and >= 2; controls aspect ratio of Wigner transform
Chapter 25. Fourier Transforms
for j=1:NN, % order of indices is p, x for faster access
end %for j
% transpose to x,p order, if desired, and scale
In Fortran 90/95, we again refer to the fictitious routine fft(isign, data). The wave-function input is in the array psi of length NN, and the output array W is a real array of dimension(NN,NN). The intermediate-storage array Wtmp is complex and of dimension(NN,NN). The variables lbd, ubd, lbdp, ubdp, rstep, and fstep are all of type integer, and pi is of type real.
if ( iand(NN, NN-1) .ne. 0 ) then write(0,) "Error: array length not power of 2" stop end if if ( size(W,1) .ne. NN .or. size(W,2) .ne. NN ) then write(0,) "Error: input and output array sizes do not match" stop end if if ( iand(n, 1) .ne. 0 .or. n .lt. 2 ) then write(0,*) "Error: n not even and positive" stop end if Wtmp = 0; do j = 1, NN ! order of indices is p, x for faster access ! do shifted copy
! do FT
call fft(-1, Wtmp(:,j))
end do ! transpose to x,p order and scale