# How to Filter 2D Fourier Coefficients

**URL:** <https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108>\
**Category:** New to Julia\
**Tags:** fftw\
**Created:** [March 25, 2024, 9:49pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108 "2024-03-25T21:49:05Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![balancedflow](https://avatars.discourse-cdn.com/v4/letter/b/22d042/32.png) [@balancedflow](https://discourse.julialang.org/u/balancedflow)\
**Post date:** [March 25, 2024, 9:49pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/1 "2024-03-25T21:49:05Z")

</div>

I have a real 2D Array (Nx,Ny) of some variable, let’s say absolute vorticity. I perform a 2d fft on the matrix by:

```julia
fc_2d = FFTW.r2r(arr,FFTW.R2HC)

```

which returns a 2d matrix of the same size of Fourier coefficients.

I want to do a 2d low pass filter by zero-ing out the scales greater than wavenumber 3 for example. I read the [documentation](https://www.fftw.org/fftw3_doc/The-Halfcomplex_002dformat-DFT.html) and I’m still a little confused about the format of ` fc_2d`. Are both dimensions in half-complex format? If so, I think I would just zero out the Fourier coefficients corresponding to wavenumbers 4 and up in both dimensions?

```julia
ind = 4
# x filter
fc_2d[ind+1:Nx-ind+1,:] .= 0 
# y filter
fc_2d[:,ind+1:Nx-ind+1] .= 0

```

Sorry if this is really basic / does not make sense. I am new to Julia and spectral analysis. Thanks.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 25, 2024, 9:57pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/2 "2024-03-25T21:57:03Z")

</div>

> [@balancedflow](#):
>
> `fc_2d = FFTW.r2r(arr,FFTW.R2HC)`

The `R2HC` format in 2d is pretty tricky to comprehend, I would recommend starting with a simple complex `fft` and then moving to the `r2c` format once you understand that in order to take advantage of the real input. (The `R2HC` transforms are typically slower than `r2c`, anyway.)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 25, 2024, 10:30pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/3 "2024-03-25T22:30:59Z")

</div>

That seems a job for GMT.jl `grdfft` but unfortunately we don’t yet have a manual tailored for the Julia wrapper (though the `grdfft` is wrapped). So far, mos complete manual is still the [CLI version](https://docs.generic-mapping-tools.org/dev/grdfft.html)

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [March 26, 2024, 9:19am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/4 "2024-03-26T09:19:43Z")

</div>

Using a `r2c` transform as suggested by @stevengj, you could do something like:

```julia
using FFTW

# Domain dimensions
Lx = 2π
Ly = 2π

# Grid size
Nx = 32
Ny = 32

# Wavenumbers (note: r2c transform is along the first dimension, so we use rfftfreq in x)
kx = rfftfreq(Nx, 2π * Nx / Lx) # only contains non-negative wavenumbers (kx ≥ 0)
ky = fftfreq(Ny, 2π * Ny / Ly)

arr = rand(Nx, Ny)
fc_2d = FFTW.rfft(arr) # perform real-to-complex (r2c) transform
summary(fc_2d) # 17×32 Matrix{ComplexF64}

# Zero-out Fourier coefficients associated to wavenumbers |k_i| > kmax in each direction.
kmax = 3
for i ∈ eachindex(kx)
    if abs(kx[i]) > kmax
        fc_2d[i, :] .= 0
    end
end
for j ∈ eachindex(ky)
    if abs(ky[j]) > kmax
        fc_2d[:, j] .= 0
    end
end

# Alternative (assuming Lx = Ly = 2π)
ind = 4
@assert kx[ind] == ky[ind] ≈ kmax
@assert kx[end] > kmax
@assert ky[Ny - ind + 2] ≈ -kmax
fc_2d[(ind + 1):end, :] .= 0 # zero-out modes kx > kmax
fc_2d[:, (ind + 1):(Ny - ind + 1)] .= 0 # zero-out modes |ky| > kmax

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 26, 2024, 9:24am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/5 "2024-03-26T09:24:20Z")

</div>

It seems like a gap in the API that there is not `rfftfreq` method for the R2HC transform. (Also that there is no convenient `rfftfreq` for 2D transforms.)

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [March 26, 2024, 9:32am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/6 "2024-03-26T09:32:20Z")

</div>

> Also that there is no convenient `rfftfreq` for 2D transforms.

I agree, it would be convenient to have a multidimensional `rfftfreq` in AbstractFFTs.jl.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [March 26, 2024, 9:38am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/7 "2024-03-26T09:38:40Z")

</div>

Wouldn’t this be the tensor product of an `rfftfreq` along the first dimension, and `fftfreq` along the others?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 26, 2024, 9:40am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/8 "2024-03-26T09:40:05Z")

</div>

Sidenote: This code essentially uses something like a Manhattan metric in Fourier space for the lowpass filter. I would intuitively do something like

```julia
fc_2d[abs2.(kx) .+ abs2.(ky') .> abs2(kmax)] .= 0

```

to zero out the higher frequencies. (Untested, double-check the transposition operator `'` which I applied to `ky` here to get a broadcast.)

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [March 26, 2024, 9:57am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/9 "2024-03-26T09:57:21Z")

</div>

That kind of Manhattan metric was what @balancedflow seemed to be after, if I understood correctly their initial post. It actually be useful, for example for dealiasing in pseudo-spectral methods.

For truncating modes |\boldsymbol{k}| \> k\_{\text{max}}, one could do something like:

```julia
kmax² = kmax^2
for j ∈ axes(us, 2), i ∈ axes(us, 1)
    k² = kx[i]^2 + ky[j]^2
    if k² > kmax²
        us[i, j] = 0
    end
end

```

Or in a more generic way (also valid for dimensions other than 2):

```julia
kmax² = kmax^2
ks = (kx, ky)
for I ∈ CartesianIndices(us)
    k⃗ = map(getindex, ks, Tuple(I))
    k² = sum(abs2, k⃗)
    if k² > kmax²
        us[I] = 0
    end
end

```

---

<div class="post-metadata">

**Author:** ![jipolanco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jipolanco/32/12129_2.png) [@jipolanco](https://discourse.julialang.org/u/jipolanco)\
**Post date:** [March 26, 2024, 9:59am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/10 "2024-03-26T09:59:11Z")

</div>

Exactly!

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 26, 2024, 11:09am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/11 "2024-03-26T11:09:55Z")

</div>

> [@balancedflow](#):
>
> I want to do a 2d low pass filter by zero-ing out the scales greater than wavenumber 3 for example.

Doing that is applying a boxcar filter in the frequency domain. A boxcar filter in frequency corresponds to a sinc filtering in space domain and that will cause ringing. The way to go is to apply a cosine-taper band pass, and that is option **F** (or `filter` in GMT.jl) of `grdfft` does. See [grdfft — GMT 6.6.0 documentation](https://docs.generic-mapping-tools.org/dev/grdfft.html#f)

---

<div class="post-metadata">

**Author:** ![balancedflow](https://avatars.discourse-cdn.com/v4/letter/b/22d042/32.png) [@balancedflow](https://discourse.julialang.org/u/balancedflow)\
**Post date:** [March 26, 2024, 10:50pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/12 "2024-03-26T22:50:21Z")

</div>

Hi everyone,

Thank you for all the help and reccomendations. I have a lot to look up and will respond back here if I have any questions!

---

<div class="post-metadata">

**Author:** ![brainandforce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brainandforce/32/211054_2.png) [@brainandforce](https://discourse.julialang.org/u/brainandforce)\
**Post date:** [March 27, 2024, 5:57am UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/13 "2024-03-27T05:57:42Z")

</div>

Perhaps you can use the package I wrote, FFTIndexing.jl, which allows you to reindex arrays based on Fourier bins.

At the moment I’m a bit too occupied to look at this problem, but I think it should be pretty easy to implement a boxcar filter with it.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 27, 2024, 1:00pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/14 "2024-03-27T13:00:32Z")

</div>

An example with GMT.jl

```julia
using GMT

# Download a Geoid grid
G = grdcut("@earth_geoid_02m", region=(-40,-20,32,45));

#Filter with Boxcar (low pass wavelength > 500 km). ´f=:g´ is to infor that grid is in geogs, `N=:a` is to select best FFT dims(
 Gf_box = grdfft(G, F="-/-/500000/500000", N=:a, f=:g);

# Filter with cosne taper. Pass wavelengths > 500 km and rejecting wavelengths < 100 km
Gf_ct = grdfft(G, F="-/-/500000/100000", N=:a, f=:g);

```

The results

```julia
viz(G, shade=true, coast=true, title="Geoid", colorbar=true, name="Goid.jpg")
viz(Gf_box, shade=true, coast=true, title="Geoid Boxcar", colorbar=true, name="Goid_box.jpg")
viz(Gf_ct, shade=true, coast=true, title="Geoid cos-taper", colorbar=true, name="Goid_ct.jpg")

```

 ![Goid](https://global.discourse-cdn.com/julialang/original/3X/a/0/a0dc32d51bbe92105806688bc9cb73f681959eee.jpeg)

 ![Goid_box](https://global.discourse-cdn.com/julialang/original/3X/1/b/1bd8575f19bc73bed0204ad453d6f2b15a9bdb29.jpeg)

 ![Goid_ct](https://global.discourse-cdn.com/julialang/original/3X/a/5/a5f30d6467b2243de897c0bc5b9f8a2953ab583a.jpeg)

---

<div class="post-metadata">

**Author:** ![balancedflow](https://avatars.discourse-cdn.com/v4/letter/b/22d042/32.png) [@balancedflow](https://discourse.julialang.org/u/balancedflow)\
**Post date:** [March 27, 2024, 9:26pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/15 "2024-03-27T21:26:00Z")

</div>

Thanks for the example code. This really helps. I was going through the documentation and trying grdfft on my data, but am having a couple of issues.

```julia
# get data
arr = nc["avo"][:,:,1] # 354x354 Matrix{Float32}
arr = mat2grid(arr); # Is it necessary to convert from a matrix to a GMT.GMTgrid here?

# do low-pass
# pass 100 km and cut out 90 km
filtered = grdfft(arr, F="-/-/100000/90000",f=:g,N=:a); 

```

This code should pass all scales \> 100 km and completely cut off scales \< 90 km, and the scales between these cutoffs are slowly filtered out, correct? Though the output of this is not what I expect. When I do something like `filtered = grdfft(arr, F="-/-/10000000/9000000",f=:g,N=:a);` I get something that looks more physically correct. These correspond to scales 10,000 km and 9,000 km (off by a factor of 2 of what my data resolution is actually in). My data resolution is 1 km so the grid domain is 354 km x 354 km.

Am I missing something obvious here? (I guess the function doesn’t know what the resolution of my data is here, since its just a 2d array of scalars)…

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [March 27, 2024, 10:57pm UTC](https://discourse.julialang.org/t/how-to-filter-2d-fourier-coefficients/112108/16 "2024-03-27T22:57:45Z")

</div>

> [@balancedflow](#):
>
> `arr = mat2grid(arr); # Is it necessary to convert from a matrix to a GMT.GMTgrid here?`

Absolutely, and what you did was not enough. You need also to pass at least the referencing information (either in the `x,y` vectors or the the `hdr` vector/matrix). The `mat2grid` help explain those things. And that is need because the `grdfft` manual also explains that.

> Description

> grdfft will take the 2-D forward Fast Fourier Transform and perform one or more mathematical operations in the frequency domain before transforming back to the space domain. An option is provided to scale the data before writing the new values to an output file. The horizontal dimensions of the grid are assumed to be in meters. Geographical grids may be used by specifying the -fflags option that scales degrees to meters.

So, your data must be in meters or in geographical coordinates. It’s in later case that one passes the `f=:g` option to inform the program that it has convert from degrees to meters. But since your original post mentioned _absolute vorticity_ I assumed that it was an oceanographic or atmospheric data. Another parameter to look at is the memory layout (row-major, column-major, top-down or botom-up). I recomend that you visualize the output of the `mat2grid` function to see if it makes sense.

EDIT

> [@balancedflow](#):
>
> `arr = nc["avo"][:,:,1]`

This makes me think that you are reading from a netCDF grid. You may save yourself the `mat2grid` hassle if you read the nc grid directly with `gmtread` or `gdalread`.