Extracting Wave Parameters from Noisy Data: A Fourier Projection Approach

Using temporal Fourier projection to resolve waves in noisy data

The Problem

You’ve got noisy spatio-temporal data from simulation or observation—atmospheric waves, turbulent flows, or another field. A 2D FFT can identify the wave in these tests, but its discrete spatial frequency bins limit wavelength resolution. If the temporal frequency is known, projection onto that frequency gives a more precise estimate of the spatial structure.

In the experiments below, projection’s typical wavelength error ranges from 0.006% to 0.102% across the tested noise levels. The full-plane 2D FFT also finds the wave, but returns 50.5 km, a fixed 1.0% error in this setup.

The Key Insight: Orthogonality as a Filter

The foundation of this method is one of the most beautiful properties in mathematics: Fourier orthogonality. You’ve probably seen this before, but let’s refresh on why it’s so powerful.

Complex exponentials form an orthogonal basis:

where $\delta_{mn}$ is the Kronecker delta (1 if $m = n$, zero otherwise).

What does this mean physically? Imagine your signal is a superposition of different frequency components. To isolate one specific frequency $\omega$, take the inner product of the signal with $e^{i\omega t}$. For Fourier basis frequencies over a matching interval, orthogonality cancels the other components. A finite observation window or an imperfect frequency estimate can still cause leakage, so frequency selection matters in practice.

This is the matched filter principle from signal processing: when you know the signal you are looking for, you can design a detector specifically for it. Under suitable assumptions about the noise, that detector can be optimal for finding the target.

And here’s the kicker for this problem: this isn’t just frequency discrimination. When you apply this to spatio-temporal data, you can use a known temporal frequency to isolate spatial structure. That’s the superpower I’m about to exploit.

The Setup

Let’s formalize the problem. You have a field $\rho(z, t)$ that varies with altitude $z$ and time $t$. This could be atmospheric density from DSMC, temperature from climate models, or any other field carrying wave information. Embedded in this noisy data is a propagating wave:

The parameters are:

• $A$: wave amplitude
• $k = 2\pi/\lambda$: vertical wavenumber (this is what I want!)
• $\omega = 2\pi/T$: angular frequency
• $\phi$: phase offset

Here’s the asymmetry in what you typically know:

Estimated first: The temporal frequency $\omega$. In these experiments I estimate it with a 1D FFT at one altitude, which works for most draws at the tested noise levels. It is not guaranteed: the estimate is wrong in 4% of draws at three times the signal and 27% at four times the signal. That failure can then corrupt the spatial phase estimate.

Unknown: The vertical wavenumber $k$. This is harder because spatial sampling is typically coarser than temporal sampling (fewer altitude points than time points). But once you have $k$, you immediately get the wavelength: $\lambda = 2\pi/k$. This tells you how the wave propagates vertically, its energy distribution, whether it satisfies dispersion relations, and much more.

So the central question is: How do you extract $k$ robustly from noisy spatio-temporal data when you already know $\omega$?

The Method: Temporal Projection

Here’s the elegant solution. For each altitude $z$, I project the time series at that altitude onto the known temporal frequency $\omega$:

This looks like computing a single Fourier coefficient—and that’s exactly what it is! But notice the clever structure: instead of transforming in both space and time simultaneously (which is what 2D FFT does), I’m:

1. Keeping altitude $z$ as a parameter (not transforming it)
2. Transforming only in time
3. Evaluating that transform at the single frequency $\omega$ I care about

The result $I(z)$ is a complex-valued function of altitude. It encodes both the amplitude and phase of the wave component at frequency $\omega$ as a function of height.

Why is this powerful? Think about what happens during that integration:

• The wave component at frequency $\omega$ oscillates in phase with $e^{i\omega t}$ → they multiply constructively → the integral accumulates coherently, growing linearly with $T$ (the total observation time)
• Noise has components at all frequencies → they oscillate randomly relative to $e^{i\omega t}$ → uncorrelated contributions add like a random walk → the integral only grows like $\sqrt{T}$ (standard deviation of $N$ uncorrelated samples scales as $\sqrt{N}$)

So the signal-to-noise ratio improves during integration by a factor of $\sqrt{T}$. The longer you observe, the better your extraction. It’s a matched filter, tuned to the temporal signature of your wave. The experiments below test this prediction across 500 noise draws at each of four noise levels. Levels 0.1, 0.5, 1, and 2 correspond to +17, +3, −3, and −9 dB, respectively. At level 2, the typical wavelength error is 0.102%.

The Math

Replace $\rho(z, t)$ in Eq. (2) with the wave form from Eq. (1), and expand the cosine using Euler’s formula: $\cos\theta = \frac{1}{2}(e^{i\theta} + e^{-i\theta})$

The first term has time phases cancel: $e^{-i\omega t} e^{i\omega t} = 1$, so: $\text{Term 1} = \frac{AT}{2} e^{i(kz + \phi)}$

The second term oscillates at $2\omega$ → integrates to zero for $n$ integer cycles.

Result:

The phase is: $\arg[I(z)] = kz + \phi$

Phase varies linearly with altitude! The slope gives the wavenumber:

Procedure:

1. Compute $\Phi(z) = \arg[I(z)]$ at each altitude
2. Unwrap phase (handle $2\pi$ jumps)
3. Linear fit → slope = $k$
4. Wavelength: $\lambda = 2\pi/k$

The orthogonality does the heavy lifting. Integration amplifies signal (grows like $T$) while noise averages down (grows like $\sqrt{T}$).

Does It Work? Synthetic Test

Theory is nice, but does this actually work? Let me test it rigorously on synthetic data where I know the ground truth.

I generated a synthetic atmospheric gravity wave with Mars-like parameters:

• Spatial domain: 100–300 km altitude, 100 grid points
• Temporal domain: 0–3600 s (1 hour), 500 time points
• True wavelength: λ = 50 km
• True period: T = 600 s (10 minutes)
• Amplitude: 10% perturbation

Then I added Gaussian noise at increasing levels to see when the method breaks.

Results

At the highest tabulated noise level, 2.0, the typical wavelength error across 500 noise draws is 0.102%. The method begins to fail as noise rises beyond about 2.5 times the signal; the frequency estimate from a single altitude is wrong in 4% of draws at three times the signal and 27% at four times the signal.

Why It Works, and Where It Breaks

Remember the matched-filter argument: during integration, the signal component at $\omega$ accumulates coherently while uncorrelated noise contributions add like a random walk. The error curve in Figure 3 compares this prediction with 500 noise draws at each level. The tested results support the prediction through noise level 2; beyond about 2.5 times the signal, single-altitude frequency selection becomes a limiting step. These experiments do not establish a general tolerance to frequency error or a minimum observation-window length.

vs. Standard 2D FFT

The comparison depends on searching the correct part of the 2D spectrum. For a traveling wave written as

the power lies in the opposite-sign quadrants, at $(k,-\omega)$ and $(-k,\omega)$. The old comparison searched only the positive-$k$, positive-$\omega$ quadrant, which contains noise for this sign convention. Its 14.43 km result was a random bin: $202/14$ km.

Searching the full frequency plane finds the wave at every tested noise level. The 2D FFT returns 50.5 km, however, a fixed 1.0% error for the 50 km wave. Its spatial wavelengths are quantized by the domain length divided by an integer; here, $202/4 = 50.5$ km. Projection gives finer wavelength resolution in these tests:

Across these four noise levels, projection’s typical errors are about 10–167 times smaller than the 2D FFT’s fixed 1.0% error. This is a resolution advantage in this synthetic setup, not evidence that a full-plane 2D FFT cannot detect the wave.

Why Projection Has Finer Resolution Here

The 2D FFT places power on discrete spatial and temporal frequency bins. In this domain, the selected spatial bin corresponds to 50.5 km, so that estimate cannot land on the 50 km truth. Projection uses the known temporal frequency to build a phase profile across altitude, then estimates the wavenumber from its slope. That fitted slope is not restricted to the 2D FFT’s discrete spatial bins, which explains its smaller errors in this test.

What About More Complex Wave Fields?

So far I’ve assumed an idealized traveling wave with constant amplitude and wavelength. But real atmospheric waves are messier. What happens when:

• The wavelength changes with altitude due to stratification?
• The wave amplitude grows or decays due to dissipation or density changes?
• The wave structure is more complex than a simple sinusoid?

The projection can describe these variations, provided the wave remains coherent at a well-defined temporal frequency and that frequency is estimated accurately enough. The derivation below shows what the projected signal contains; it does not guarantee accurate recovery in every noisy or real-world field.

Varying Wavelength with Altitude

In a stratified atmosphere, the vertical wavenumber $k(z)$ can change with altitude as the background density, temperature, and buoyancy frequency vary. Instead of a constant $k$, the wave might have:

When you apply the projection method, you still get:

The phase is no longer linear in $z$—it’s the integral of the local wavenumber. But you can still extract the local wavenumber by taking the derivative:

So instead of fitting a straight line to the phase profile, you can:

• Compute the local slope at each altitude (using finite differences or smoothing splines)
• Track how $k(z)$ varies with height
• Identify where the wave refracts, reflects, or changes character

The method gives you the full spatial structure of $k(z)$, not just a single average value.

Amplitude Variation and Dissipation

Real waves don’t have constant amplitude. In decreasing atmospheric density, wave amplitude typically grows with altitude (energy conservation). Or dissipation might cause exponential decay. Your wave might look like:

The projection gives:

Now the magnitude varies: $|I(z)| = A(z)T/2$. You can extract:

• Phase: $\arg[I(z)] = kz + \phi$ → still gives wavenumber from the slope
• Amplitude envelope: $A(z) = 2|I(z)|/T$ → tracks how amplitude varies with altitude

So the same extraction gives you both the spatial wavelength structure and the amplitude envelope. With $A(z)$ and $k(z)$ in hand, you can:

• Identify dissipation rates from exponential fits to $A(z)$
• Check energy conservation (in decreasing density, $A \propto \rho^{-1/2}$)
• Locate wave breaking or saturation (where $A$ grows faster than expected)

What Parameters Handle This?

Looking back at the derivation in “The Math” section:

The beauty is that $A$, $k$, and $\phi$ were never assumed to be constant—they were always implicitly functions of $z$! The derivation holds locally at each altitude. The projection integral naturally captures:

• $A(z)$ in the magnitude $|I(z)|$
• $k(z)$ in the phase gradient $d\arg[I(z)]/dz$
• $\phi(z)$ in the absolute phase

The key assumption is that the temporal frequency $\omega$ is well-defined and consistent across altitudes. The synthetic tests in this post use a coherent wave; they do not establish performance for every real atmospheric field. As shown above, high noise can also make the single-altitude frequency estimate fail and corrupt the recovered spatial phase.

Quick Implementation

# 1. Extract omega from temporal FFT
z_mid = len(z) // 2
freqs = fftfreq(len(t), dt)
omega = 2 * np.pi * freqs[np.argmax(np.abs(fft(field[z_mid, :])))]

# 2. Project at each altitude
I_z = np.array([np.trapz(field[i, :] * np.exp(1j*omega*t), t)
                for i in range(len(z))])

# 3. Extract wavenumber from phase slope
phase = np.unwrap(np.angle(I_z))
k = np.polyfit(z, phase, 1)[0]
wavelength = 2 * np.pi / k

Part of my research notes on atmospheric wave analysis. For more on DSMC and wave dynamics, see other posts in this series.