How to Filter 2D Fourier Coefficients

New to Julia
1

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:

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 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?

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.

2

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.)

1 Like
3

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

1 Like
4

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

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
1 Like
5

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.)

6

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

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

7

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

8

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

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.)

1 Like
9

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:

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):

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
2 Likes
10

Exactly!

11

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

3 Likes
12

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!

13

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.

1 Like
14

An example with GMT.jl

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

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
Goid2029×1384 390 KB

Goid_box
Goid_box2029×1384 223 KB

Goid_ct
Goid_ct2029×1384 201 KB

2 Likes
15

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.

# 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)…

16

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

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.

1 Like