Understanding FFT Scaling with Matlab (or Python, or…)

The Fourier Transform is one of the most frequently used computational tools in earthquake seismology. Using an FFT requires some understanding of the way the information is encoded (frequency ordering, complex values, real values, etc) and these are generally well documented in the various software packages used in the field. I’ll assume that you have a good handle on these issues and that you know the basic idea about using forward and inverse transforms. I want to focus on an aspect of FFT use that often confuses students, the scaling provided by the transforms. Scaling can be confusing because to make the functions more efficient and flexible, the scaling is often omitted, or it is implicitly included by assuming that you plan to forward and inverse transform the same signal.

In seismology we usually start with a signal that has a specified sample rate, \(\delta t\), and a length, \(N\). I’ll assume that we are dealing with a signal for which \(N\) is a power of two, which used to be the case whenever we used FFT’s but is less required now that computers are much faster. FFT routines don’t use a specific value of sample rate, \(\delta t\), at all, but they don’t always completely ignore it – they sometimes include \(\delta t\) implicitly in their scaling. For most cases, when you compute a forward transform, you have to scale the result by the sample rate to get the correct values (this \(\delta t\) arises from the \(\delta t\) in the Fourier Transform that the DFT is approximating).

Consider the case in Matlab. We’ll construct a unit area “spike” in Matlab by creating a signal with a single non-zero value. So we can keep track of the physical units, which are also ignored in the computation but are essential to keep track of when doing science, I’ll assume that our signal is a displacement seismogram with units of meters and that the units of \(\delta t\) are seconds. To create a signal with unit area, we must set the amplitude equal to \(1 / \delta t ~\mbox{meters}\). That is the spike amplitude has the same numerical value as \(1/\delta t\), but units of length. For a spike with an amplitude \(1/\delta t\), the area under the curve in the time domain is (using a the formula for the area of a triangle with width 2 dt):

\[Area = \frac{1}{2} \cdot 2 \times \delta t~\mbox{[s]} \cdot \frac{1}{\delta t}~\mbox{[m]}= 1 \mbox{ [m-s] .} \]

The Fourier Transform of a unit area spike should have a value of unity. Let’s see what we get in Matlab.

 >> dt = 0.1;
 >> s = [1/dt,0,0,0,0,0,0,0];
 >> shat = fft(s) % amplitudes are wrong without the correct scaling
 shat = 10 10 10 10 10 10 10 10

The spectrum is a factor of \(1/\delta t\) too large. We get the correct result  by multiplying by \(\delta t\).

>> shat = fft(s) * dt % you must scale by dt to get the correct values shat = 1 1 1 1 1 1 1 1

When we inverse transform, we have to account for the \(\delta f\) in the Fourier Transform Integral. That means we multiply by \(\delta f\). For a signal of length \(N\) and sample rate \(\delta t\), the frequency sampling rate, \(\delta f\) is

\[\delta f = \frac{1}{N \cdot \delta t}~. \]

Here’s the tricky part, Matlab multiplies the inverse transform by \(1/N = \delta t \cdot \delta f\). So we have to correct for that by either dividing by \(\delta t\) or simply remove the \(1/N\) factor and apply the \(\delta f\), which makes things easier to understand. Here’s a summary, I explicitly multiply the results of the inverse transform (ifft) by \(N\) to remove Matlab’s scaling, then I scale by \(\delta f\).

dt = 0.1;
N = 8;
df = 1/(N*dt);
% define a unit area signal
s = [1/dt,0,0,0,0,0,0,0]
% forward transform - scale by dt
shat = fft(s) * dt
% inverse transform - remove 1/N factor, then scale by df
ifs = (N * ifft(shat)) * df


The first thing that I do when I start using a new FFT, is to test the amplitudes using a unit-area delta-function approximation (which is really a triangle). The unit area time domain signal must have a value of unity at the zero frequency, and if it’s a spike, the spectrum should be flat. If you think this is not that critical, then you need to think about deconvolutions, earthquake-source spectra calculations, etc. To get the correct amplitudes, you have to apply your scaling (unless you do a forward and then inverse transform without applying any operations using physical quantities to the original signal).

There’s nothing special about Matlab here – you can use the same tests to understand the scaling in Mathematica, Python, Fortran, C, JavaScript, …

Seismic Reflection/Transmission Coefficients with Matlab

A common computation in seismology is the calculation of reflection and transmission coefficients that describe the partitioning of energy when a seismic wave strikes a boundary between elastic materials. The coefficients depend on the seismic wave speeds and densities on either side of the boundary and on the incident wave horizontal slowness, or ray parameter. The formulas for the values can be found in Table 3.1 of Lay and Wallace (1995) (note that an error in the sign of the second term in the computation of b) or in equations 5.38 and 5.39 of Aki and Richards (1980). For a welded elastic boundary, an incident P or SV wave results in four waves, two reflected and two transmitted or refracted waves. Thus for each case we need four coefficients – notation for each is illustrated in the figure below.

Notation used to identify reflection and transmission coefficients. We must consider two cases, that for an incident P wave and that for an incident SV wave. R identifies “Reflected” waves, T represents “Transmitted” or “refracted” waves.

Continue reading