Fourier Transform of Signals: Theory, Computation and Physical Interpretation

1. From a Measured Signal to Its Frequency Components

A signal is normally recorded as a quantity that changes with time. Examples include an electrical voltage $V(t)$, acoustic pressure $P(t)$, optical intensity $I(t)$, or a terahertz electric field $E(t)$. In a THz time-domain experiment, the detector directly provides a sequence of electric-field values measured at successive time delays. The resulting waveform may contain a main pulse, oscillations, reflections, noise, delayed echoes, or several pulses. The purpose of the Fourier transform is to determine which frequencies construct this measured temporal waveform and what phase each frequency possesses. Thus, the time-domain representation answers the question “what happens as a function of time?”, while the frequency-domain representation answers “which frequencies are present, with what amplitude and phase?”

The continuous Fourier transform of a signal $x(t)$ is

$X(f)=\int_{-\infty}^{+\infty}x(t)e^{-i2\pi ft},dt$

and the inverse Fourier transform is

$x(t)=\int_{-\infty}^{+\infty}X(f)e^{i2\pi ft},df$

The function $x(t)$ is the time-domain signal and $X(f)$ is its complex frequency-domain representation. The exponential term can be understood from Euler’s relation

$e^{-i2\pi ft}=\cos(2\pi ft)-i\sin(2\pi ft)$

so the Fourier transform effectively compares the measured signal with sine and cosine waves of every possible frequency. If a certain frequency is strongly present in the data, contributions at that frequency add coherently and produce a large Fourier coefficient. If that frequency is absent, the positive and negative contributions largely cancel.

Suppose the measured waveform contains two frequencies,

$x(t)=A_1\cos(2\pi f_1t+\phi_1)+A_2\cos(2\pi f_2t+\phi_2)$

The temporal waveform may look complicated because the two oscillations interfere. In the Fourier domain, however, the components separate into peaks around $f_1$ and $f_2$. The Fourier transform therefore performs a mathematical decomposition of the signal into individual harmonic components.

The Fourier result is generally complex and may be written as

$X(f)=\operatorname{Re}[X(f)]+i\operatorname{Im}[X(f)]$

or, more usefully,

$X(f)=|X(f)|e^{i\phi(f)}$

where the spectral amplitude is

$A(f)=|X(f)|=\sqrt{\operatorname{Re}[X(f)]^2+\operatorname{Im}[X(f)]^2}$

and the spectral phase is

$\phi(f)=\arg[X(f)]=\operatorname{atan2}\left(\operatorname{Im}[X(f)],\operatorname{Re}[X(f)]\right)$

The amplitude $|X(f)|$ tells us how strongly a particular frequency contributes to the signal. The phase $\phi(f)$ describes how that frequency is positioned relative to the other spectral components. Both amplitude and phase are required to reconstruct the original waveform correctly.

This point is extremely important. Two signals may have exactly the same spectral amplitude but different spectral phases and therefore completely different temporal shapes. A Fourier spectrum is therefore not fully described by amplitude alone. The complete information is contained in the complex Fourier spectrum.

A simple temporal shift demonstrates the importance of phase. If the original signal is delayed by a time $\tau$,

$x(t)\rightarrow x(t-\tau)$

then its Fourier transform becomes

$X(f)\rightarrow X(f)e^{-i2\pi f\tau}$

and therefore the spectral phase changes according to

$\phi(f)\rightarrow\phi(f)-2\pi f\tau$

The amplitude spectrum remains unchanged. Therefore, moving a THz pulse in time changes the slope of its spectral phase but does not change its Fourier amplitude. The delay can be determined from the phase slope as

$\tau=-\frac{1}{2\pi}\frac{d\phi}{df}$

or, using angular frequency $\omega=2\pi f$,

$\tau=-\frac{d\phi}{d\omega}$

This relation is particularly important in THz spectroscopy because spectral phase contains propagation delay, dispersion and timing information. Removing a linear phase slope is mathematically equivalent to shifting the pulse in time. Therefore, one should not automatically flatten spectral phase unless the goal is deliberately to remove temporal delay.

Short signals are generally broadband because many frequencies are required to construct a localized waveform. A long sinusoidal signal can contain a very narrow frequency range, whereas an ultrashort pulse requires a broad spectrum. An ideal infinitely short temporal impulse contains all frequencies. This time-frequency relationship is one of the fundamental concepts behind ultrafast optics and broadband THz science.

2. How the Fourier Transform Is Actually Performed Computationally

A computer does not receive a continuous function $x(t)$. It receives a finite sequence of measured numbers. Suppose an experiment records $N$ samples,

$x[0],x[1],x[2],\ldots,x[N-1]$

at equally spaced times

$t_n=n\Delta t$

where $\Delta t$ is the temporal sampling interval. The computer therefore works with the Discrete Fourier Transform, or DFT,

$X[k]=\sum_{n=0}^{N-1}x[n]e^{-i2\pi kn/N}$

where $k$ is the frequency-bin index,

$k=0,1,2,\ldots,N-1$

The inverse DFT is

$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{i2\pi kn/N}$

In practice, we normally calculate the DFT using the Fast Fourier Transform, or FFT. The FFT is not a different transform. It is simply an efficient numerical algorithm for calculating exactly the same DFT.

The computational procedure should be understood as a sequence of clearly defined steps. Consider that your experimental data contain two columns: one column containing time and one containing measured electric field. The first step is to load these arrays,

$t=[t_0,t_1,t_2,\ldots,t_{N-1}]$

and

$E=[E_0,E_1,E_2,\ldots,E_{N-1}]$

The time points should normally be equally spaced. The temporal sampling interval is calculated from

$\Delta t=t_1-t_0$

or more robustly from the mean,

$\Delta t=\operatorname{mean}(t_{n+1}-t_n)$

If the time axis is measured in picoseconds, it should normally be converted to seconds before obtaining frequencies in hertz,

$t_{\mathrm{s}}=t_{\mathrm{ps}}\times10^{-12}$

The sampling frequency is then

$f_s=\frac{1}{\Delta t}$

The total measured temporal duration is approximately

$T=N\Delta t$

and the Fourier frequency spacing is

$\Delta f=\frac{f_s}{N}=\frac{1}{N\Delta t}\approx\frac{1}{T}$

The maximum frequency that can be uniquely represented is the Nyquist frequency,

$f_{\mathrm{Nyq}}=\frac{f_s}{2}=\frac{1}{2\Delta t}$

These relations should be calculated before performing the FFT, because they determine the physical frequency axis.

Now consider a numerical example. Suppose the measured time spacing is

$\Delta t=0.05~\mathrm{ps}=5.0\times10^{-14}~\mathrm{s}$

The sampling frequency becomes

$f_s=\frac{1}{5.0\times10^{-14}}=2.0\times10^{13}~\mathrm{Hz}=20~\mathrm{THz}$

and therefore

$f_{\mathrm{Nyq}}=10~\mathrm{THz}$

If the measurement contains $N=2000$ temporal samples, then

$T=N\Delta t=2000\times0.05~\mathrm{ps}=100~\mathrm{ps}$

and the frequency spacing becomes

$\Delta f=\frac{1}{100~\mathrm{ps}}=0.01~\mathrm{THz}$

Thus, the computational Fourier spectrum will contain discrete frequencies separated by approximately $0.01~\mathrm{THz}$.

The actual Fourier transform is then calculated computationally. In Python this may be written as

E_FD = np.fft.fft(E_TD)

where E_TD contains the measured time-domain data and E_FD contains the complex Fourier coefficients.

However, the FFT output alone does not automatically provide the physical frequency axis. The corresponding frequency values must also be generated,

freq_Hz = np.fft.fftfreq(N, d=dt_s)

If frequencies in terahertz are desired,

freq_THz = freq_Hz / 1e12

The raw FFT contains both positive and negative frequencies. For a real-valued experimental waveform, these two parts contain redundant information. Therefore, it is often more convenient to use the real FFT,

E_FD = np.fft.rfft(E_TD)
freq_Hz = np.fft.rfftfreq(N, d=dt_s)
freq_THz = freq_Hz / 1e12

The rfft output contains frequencies from zero to the Nyquist frequency.

The amplitude spectrum is calculated from the absolute value of the complex FFT,

amplitude = np.abs(E_FD)

and the spectral phase is calculated using

phase = np.angle(E_FD)

Because np.angle returns phase wrapped between approximately $-\pi$ and $+\pi$, a continuous phase can be obtained using

phase_unwrapped = np.unwrap(phase)

The fundamental computational chain is therefore

$E(t)\rightarrow\mathrm{FFT}\rightarrow E(f)\rightarrow|E(f)|,\phi(f)$

or explicitly,

$\text{measured samples}\rightarrow\text{complex FFT}\rightarrow\text{frequency axis}\rightarrow\text{amplitude and phase}$

A compact Python example is

import numpy as np
import matplotlib.pyplot as plt

time_ps = np.asarray(time_ps)
E_TD = np.asarray(E_TD)

dt_s = np.mean(np.diff(time_ps)) * 1e-12
N = len(E_TD)

E_FD = np.fft.rfft(E_TD)
freq_THz = np.fft.rfftfreq(N, d=dt_s) / 1e12

amplitude = np.abs(E_FD)
phase = np.angle(E_FD)
phase_unwrapped = np.unwrap(phase)

fig, ax = plt.subplots(1, 2, figsize=(9, 3.5))

ax[0].plot(freq_THz, amplitude)
ax[0].set_xlabel("Frequency (THz)")
ax[0].set_ylabel("Amplitude")

ax[1].plot(freq_THz, phase_unwrapped)
ax[1].set_xlabel("Frequency (THz)")
ax[1].set_ylabel("Unwrapped phase (rad)")

plt.tight_layout()
plt.show()

Computationally, every Fourier-bin index $K$ corresponds to a physical frequency. For an evenly sampled signal,

$f_K=K\Delta f$

approximately for the positive-frequency region. Therefore, if

$\Delta f=0.0746~\mathrm{THz}$

then

$K=13\rightarrow f\approx13\times0.0746=0.970~\mathrm{THz}$

$K=14\rightarrow f\approx14\times0.0746=1.044~\mathrm{THz}$

$K=15\rightarrow f\approx15\times0.0746=1.119~\mathrm{THz}$

This explains exactly what the variable $K$ means in numerical FFT analysis. $K$ is not itself frequency; it is the index identifying one discrete Fourier bin. The actual physical frequency is obtained from the FFT frequency array.

For experimental signals it is often useful to remove a constant baseline before performing the FFT,

E_TD = E_TD - np.mean(E_TD)

This suppresses an artificial DC component caused by a nonzero detector offset. However, one should only remove the baseline when it does not represent meaningful low-frequency physics.

FFT normalization also requires attention. The direct output from np.fft.rfft scales with the number of temporal samples. A basic normalized amplitude may be calculated as

amplitude = np.abs(E_FD) / N

For a one-sided sinusoidal amplitude spectrum, positive-frequency components other than DC and the Nyquist component are often multiplied by two,

amplitude = np.abs(E_FD) / N
amplitude[1:-1] *= 2

The exact normalization depends on whether the goal is to display relative spectral amplitude, recover sinusoidal peak amplitude, calculate signal energy, or calculate a power spectral density. For many THz spatial-field plots, spectra are simply normalized by their maximum value,

amplitude_norm = amplitude / np.max(amplitude)

The inverse FFT provides an important verification. If the complex spectrum has not been modified, the original temporal waveform can be reconstructed using

E_reconstructed = np.fft.irfft(E_FD, n=N)

Mathematically,

$E(t)=\mathcal{F}^{-1}{E(f)}$

The reconstructed waveform should agree with the original sampled waveform within numerical precision. This demonstrates an essential point: the Fourier transform does not destroy temporal information. The complete complex spectrum contains everything required to recover the original signal.

3. Sampling, Frequency Resolution, Phase, Windowing and Common Computational Problems

The accuracy of a Fourier transform depends strongly on how the original signal was sampled. The two most important experimental parameters are the temporal sampling interval $\Delta t$ and total acquisition window $T$. They control different aspects of the spectrum. The sampling interval controls the highest measurable frequency,

$f_{\max}\approx\frac{1}{2\Delta t}$

while the total recorded time controls frequency-bin spacing,

$\Delta f\approx\frac{1}{T}$

Therefore, if you want to measure higher THz frequencies, you require smaller temporal steps. If you want finer frequency resolution, you require a longer temporal measurement window.

Aliasing occurs when the experimental signal contains frequencies higher than the Nyquist frequency. Consider a sampling rate $f_s$. Any signal component above

$f_s/2$

cannot be represented uniquely. Such frequencies fold back into the lower-frequency spectrum and appear as false components. Therefore,

$f_{\mathrm{signal,max}}<\frac{f_s}{2}$

should be satisfied.

Frequency resolution must also be understood correctly. Suppose two spectral resonances occur at $1.00~\mathrm{THz}$ and $1.01~\mathrm{THz}$. If the measurement window is only $20~\mathrm{ps}$,

$\Delta f=\frac{1}{20~\mathrm{ps}}=0.05~\mathrm{THz}$

so these two frequencies cannot be clearly resolved. Increasing the sampling density without increasing the total window does not solve this problem. To achieve approximately $0.01~\mathrm{THz}$ bin spacing, the observation window needs to be roughly

$T\approx\frac{1}{0.01~\mathrm{THz}}=100~\mathrm{ps}$

A related concept is zero padding. If a measured signal contains $N$ samples, zeros may be appended before performing the FFT,

N_fft = 4 * N
E_FD = np.fft.rfft(E_TD, n=N_fft)

This produces more plotted frequency points, giving a smoother-looking spectrum. However, it does not add experimental information and does not improve the true physical spectral resolution. The true resolution is still primarily determined by the original acquisition duration $T$.

Another practical problem is spectral leakage. The FFT mathematically assumes that the measured data segment repeats periodically. If the beginning and end of the measured signal do not connect smoothly, the repeated waveform contains an artificial discontinuity. This discontinuity generates additional Fourier frequencies. A pure sinusoid can therefore spread into neighboring frequency bins instead of appearing as a perfectly narrow spectral peak.

Window functions reduce this boundary discontinuity. A windowed signal is

$x_w[n]=x[n]w[n]$

where $w[n]$ may be a Hann, Hamming, Blackman, or another window. For example,

window = np.hanning(N)
E_windowed = E_TD * window
E_FD = np.fft.rfft(E_windowed)

Windowing reduces spectral side lobes but broadens spectral features, meaning that there is a trade-off between leakage suppression and spectral resolution.

For a short THz pulse that naturally returns close to zero well before the boundaries of the time window, strong windowing may not be necessary. However, echoes or long oscillations extending toward the beginning or end of the acquisition can produce leakage. A window should therefore be applied because the signal requires it, not simply because FFT analysis is being performed.

Spectral phase must also be treated carefully. The computational phase is obtained from

$\phi(f)=\operatorname{atan2}(\operatorname{Im}[E(f)],\operatorname{Re}[E(f)])$

which is wrapped between $-\pi$ and $\pi$. The wrapped phase may contain sudden jumps from approximately $+\pi$ to $-\pi$, even when the physical phase varies continuously. Unwrapping removes these artificial $2\pi$ jumps,

phase_unwrapped = np.unwrap(np.angle(E_FD))

Phase should generally only be interpreted where the spectral amplitude is sufficiently large. If

$|E(f)|\rightarrow0$

then both the real and imaginary Fourier components are dominated by noise, making the calculated phase essentially arbitrary. A useful phase mask can therefore be defined using a relative amplitude threshold,

amp_norm = amplitude / np.max(amplitude)
phase_plot = np.ma.array(phase_unwrapped, mask=amp_norm < 0.05)

The threshold should be selected from the signal-to-noise characteristics of the experiment rather than simply to make the phase appear smoother.

Another important computational issue is choosing the temporal zero. If a pulse is centered at a different time, its Fourier amplitude is unchanged but its spectral phase receives a linear slope. Therefore, if two measurements are being compared for material dispersion or intrinsic phase differences, they may first need to be temporally aligned. If the purpose is instead to determine propagation delay, that linear phase slope should be retained because it contains the desired physical information.

A reflection or echo in the time domain produces characteristic modulation in the frequency domain. Suppose

$E(t)=E_0(t)+aE_0(t-\tau)$

The Fourier transform becomes

$E(f)=E_0(f)\left[1+ae^{-i2\pi f\tau}\right]$

The interference between the original pulse and delayed copy produces spectral fringes. Their approximate spacing is

$\Delta f_{\mathrm{fringe}}\approx\frac{1}{\tau}$

Thus, a reflection separated by a large temporal delay produces closely spaced spectral fringes, while a short-delay reflection produces widely spaced fringes. This is an excellent example of the complementary relationship between time and frequency.

4. Fourier Transform of THz Spatio-Temporal Data

The same computational procedure extends directly from a one-dimensional waveform $E(t)$ to a spatially resolved THz dataset $E(x,y,t)$. Suppose an experiment records $N_t$ temporal frames, and every frame contains an $N_y\times N_x$ image. The data array therefore has dimensions

$(N_t,N_y,N_x)$

Every spatial pixel contains a complete temporal THz waveform,

$E_{ij}(t)=E(x_i,y_j,t)$

The Fourier transform is performed independently along the temporal dimension for every pixel,

$E(x,y,f)=\mathcal{F}_t{E(x,y,t)}$

Computationally, if time is axis zero,

E_FD = np.fft.rfft(E_TD, axis=0)

The frequency array is obtained from

Nt = E_TD.shape[0]
freq_THz = np.fft.rfftfreq(Nt, d=dt_s) / 1e12

The resulting array now has dimensions approximately

$(N_f,N_y,N_x)$

where $N_f$ is the number of positive-frequency bins.

The spectral amplitude map is

$A(x,y,f)=|E(x,y,f)|$

and computationally,

E_FD_amp = np.abs(E_FD)

The spectral phase map is

$\Phi(x,y,f)=\arg[E(x,y,f)]$

and computationally,

E_FD_phase = np.angle(E_FD)

If a selected Fourier index is

$K=15$

then

amp_map = E_FD_amp[K]
phase_map = E_FD_phase[K]
frequency = freq_THz[K]

This extracts the complete two-dimensional amplitude and phase distributions corresponding to one discrete THz frequency.

The key point is that the spatial coordinates are not transformed in this operation. Only the temporal axis is Fourier transformed. Therefore,

$(x,y,t)\rightarrow(x,y,f)$

The result is a frequency-resolved spatial electric field.

For a THz vortex beam, a particular frequency might show a doughnut-like amplitude distribution and a helical phase distribution. For an $\ell=1$ vortex,

$\Phi(x,y)\approx\phi$

where

$\phi=\operatorname{atan2}(y-y_0,x-x_0)$

and the phase changes by approximately

$\Delta\Phi=2\pi$

around a complete circle enclosing the vortex center.

Performing the same extraction for every frequency allows the spatial mode to be studied as a function of frequency,

$E(x,y,f_1),E(x,y,f_2),E(x,y,f_3),\ldots$

This is appropriately described as frequency-resolved THz field reconstruction or spectral decomposition of the spatio-temporal THz field.

If only one Fourier bin is retained and transformed back to time, the reconstructed signal behaves approximately like a monochromatic oscillation. For example,

$E_K(x,y,t)=A_K(x,y)\cos[2\pi f_Kt+\Phi_K(x,y)]$

Such a signal is not temporally localized because a truly narrow frequency component must extend over a long time.

To reconstruct a temporally localized frequency range, select a finite bandwidth instead,

$E_{\mathrm{band}}(x,y,f)=E(x,y,f)W(f)$

and then calculate

$E_{\mathrm{band}}(x,y,t)=\mathcal{F}^{-1}{E_{\mathrm{band}}(x,y,f)}$

Computationally,

E_filtered = np.zeros_like(E_FD)

mask = (freq_THz >= 0.8) & (freq_THz <= 1.2)

E_filtered[mask] = E_FD[mask]

E_band_TD = np.fft.irfft(E_filtered, n=Nt, axis=0)

This reconstructs the spatial-temporal field contributed by frequencies between $0.8$ and $1.2~\mathrm{THz}$.

The full computational workflow for experimental THz data can therefore be summarized as

$E(x,y,t)\rightarrow\mathrm{FFT~along~}t\rightarrow E(x,y,f)$

followed by

$E(x,y,f)\rightarrow|E(x,y,f)|$

and

$E(x,y,f)\rightarrow\Phi(x,y,f)$

followed, when required, by frequency selection,

$E(x,y,f)W(f)$

and reconstruction,

$E(x,y,f)W(f)\rightarrow\mathrm{IFFT}\rightarrow E_{\mathrm{selected}}(x,y,t)$

The Fourier transform therefore does considerably more than produce a frequency spectrum. Computationally it converts a measured temporal dataset into a complex spectral representation containing amplitude and phase. In THz science this makes it possible to isolate selected spectral components, calculate wavefront phase, determine arrival times, investigate dispersion, analyze vortex structure, reconstruct frequency bands, compare spatial modes at different frequencies, and propagate measured fields numerically. The central concept is simple but powerful: the FFT converts the measured sequence of temporal electric-field samples into complex frequency coefficients, and those coefficients contain the complete spectral amplitude and phase information required to reconstruct and understand the original THz signal.