# Extracting wave numbers using FFT and 2D data

**URL:** <https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758>\
**Category:** General Usage\
**Tags:** fftw\
**Created:** [April 6, 2025, 2:17pm UTC](https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758 "2025-04-06T14:17:02Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [April 6, 2025, 2:17pm UTC](https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758/1 "2025-04-06T14:17:02Z")

</div>

Dear all,

I am trying to write a basic script that extract sthe wave numbers of a 2D signal (defined in ‘spatial coordinates’) using FFTW.  
To test the script, I want to define a signal where the wave numbers in x and y are predefined. Then, I would like to perform a spectral analysis and I expect that this analysis returns the a priori defined wave numbers (for each spatial dimension).  
So far my script fails achieving this. Does anyone see obvious error or misuses of FFTW?

Thanks and cheers!

```julia
# Define the 2D signal with known wavenumbers
Lx, Ly = 10, 10 # Domain size in x and y directions
Nx, Ny = 256, 256 # Number of samples in x and y directions
dx, dy = Lx / Nx, Ly / Ny # Spatial step sizes

# Coordinate arrays
x, y = LinRange(0, Lx, Nx), LinRange(0, Ly, Ny)

# Define the wavenumbers in the x and y directions
kx, ky = 2 * π * 10 / Lx, 2 * π * 2 / Ly   

# Create a signal with two sinusoidal components in the 2D space
signal = cos.(kx * x .+ ky * y')

# Compute the 2D FFT of the signal 
fft_signal = fft(signal)

# Frequency bins in x and y directions
fx = fftfreq(Nx, dx)  
fy = fftfreq(Ny, dy)  

# Shift the FFT (move zero frequency to the center)
fft_shifted = fftshift(fft_signal)

# Magnitude 
magnitude = abs.(fft_shifted)

# Find dominant frequencies
max_idx = argmax(magnitude)

# Convert Cartesian index to row and column
max_row, max_col = Tuple(max_idx)

# Get the corresponding (fx, fy) of the peak (frequencies)
peak_fx = fx[max_col]
peak_fy = fy[max_row]

# Wavenumber in x and y
detected_kx = 2 * π * peak_fx
detected_ky = 2 * π * peak_fy  

# Print detected and input wavenumbers
println("Detected wavenumber in x: $detected_kx cycles/unit")
println("Detected wavenumber in y: $detected_ky cycles/unit")

println("Input wavenumber in x: $kx cycles/unit")
println("Input wavenumber in y: $ky cycles/unit")

p1 = heatmap(x, y, signal', xlabel="x", ylabel="y", title="Original 2D Signal", aspect_ratio=1)
p2 = heatmap(sort(fx), sort(fy), magnitude, xlabel="Frequency in x", ylabel="Frequency in y", title="2D Magnitude Spectrum of FFT", aspect_ratio=1)
plot(p1, p2)

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 6, 2025, 5:18pm UTC](https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758/2 "2025-04-06T17:18:00Z")

</div>

The issue is in the way the bins are computed, I believe there are two mistakes. The following should fix it:

```julia
# Frequency bins in x and y directions
fx = fftshift(fftfreq(Nx, 1/dx))  
fy = fftshift(fftfreq(Ny, 1/dy))

```

You need also to multiply by 2\*pi the detected wavenumbers:

```julia
# Print detected and input wavenumbers
println("Detected wavenumber in x: $(2*pi*detected_kx) cycles/unit")
println("Detected wavenumber in y: $(2*pi*detected_ky) cycles/unit")

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [April 6, 2025, 6:12pm UTC](https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758/4 "2025-04-06T18:12:48Z")

</div>

Also, your spatial sampling `dx` and `dy` are not well computed, as you can confirm by doing `diff(x)`.

I recommend you do:

```julia
x, y = LinRange(0, Lx, Nx+1)[1:Nx], LinRange(0, Ly, Ny+1)[1:Ny]

```

The resulting dx, dy should be as in your:

```julia
dx, dy = Lx / Nx, Ly / Ny

```

and it should have better endpoint properties.

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [April 6, 2025, 7:13pm UTC](https://discourse.julialang.org/t/extracting-wave-numbers-using-fft-and-2d-data/127758/5 "2025-04-06T19:13:36Z")

</div>

Thanks a lot! This is all working now. Here’s the updated code:

```julia

# Define the 2D signal with known wavenumbers
Lx, Ly = 10, 10 # Domain size in x and y directions (in arbitrary units)
Nx, Ny = 256, 256 # Number of samples in x and y directions
dx, dy = Lx / Nx, Ly / Ny # Spatial step sizes

# Generate the 1D grid for x and y, ensure they are sorted in ascending order
x, y = LinRange(0+dx/2, Lx-dx/2, Nx), LinRange(0+dy/2, Ly-dy/2, Ny)

# Define the wavenumbers in the x and y directions
kx, ky = 2 * π * 10 / Lx, 2 * π * 2 / Ly   

# Create a signal with two sinusoidal components in the 2D space
signal = cos.(kx * x .+ ky * y')

# Compute the 2D FFT of the signal using fft
fft_signal = fft(signal)

# Frequency bins in x and y directions using fftfreq
fx = fftshift(fftfreq(Nx, 1/dx))  
fy = fftshift(fftfreq(Ny, 1/dy))

# Shift the FFT output for better visualization (move zero frequency to the center)
fft_shifted = fftshift(fft_signal)

# Compute the magnitude of the 2D FFT
magnitude = abs.(fft_shifted)

# Find dominant frequencies
max_idx = argmax(magnitude)

# Convert Cartesian index to row and column
max_row, max_col = Tuple(max_idx)

# Get the corresponding (fx, fy) of the peak (frequencies)
peak_fx = fx[max_row]
peak_fy = fy[max_col]

# Wavenumber k = 2 * π * frequency
detected_kx = 2 * π * peak_fx # Wavenumber in the x-direction (physical space)
detected_ky = 2 * π * peak_fy # Wavenumber in the y-direction (physical space)

# Print detected and input wavenumbers
println("Detected wavenumber in x: $(detected_kx) cycles/unit")
println("Detected wavenumber in y: $(detected_ky) cycles/unit")

println("Input wavenumber in x: $kx cycles/unit")
println("Input wavenumber in y: $ky cycles/unit")

p1 = heatmap(x, y, signal', xlabel="x", ylabel="y", title="Original 2D Signal", aspect_ratio=1)
p2 = heatmap(sort(fx), sort(fy), magnitude, xlabel="Frequency in x", ylabel="Frequency in y", title="2D Magnitude Spectrum of FFT", aspect_ratio=1)
plot(p1, p2)

```