Hi everyone.
I have some 2D data in time and space. I am trying to do some Frequency-wavenumber analysis. Anyone has a code to compute the proper Frequency-wavenumber spectrum?
Let me know.
Kind regards.
M.
You mean you want to compute a FFT? Have a look at FFTW.jl 
I want a 2D FFT and then I need to do something like fftshift in 2D. I am unsure this is enough to achieve:
Then you use FFTW.jl and the functions fft (or one of its variants) and fftshift.
Using FFTW.jl and Plots.jl you can get a simple FK plot like this for two waves travelling at different speeds:
Code
# 1 - Define an oscillatory signal:
f(t,t0) = (t>t0) ? sin(2π*30*(t - t0)) * exp(-20*(t - t0)) : 0
# 2- Create input wavefield with two propagating waves:
Nt, Nx = 1001, 201
t = LinRange(0., 2., Nt) # time (s)
x = LinRange(0., 1000., Nx) # distance (m)
Vp, Vs = 2000., 1100. # velocity (m/s)
Wt = f.(t', x/Vp) + f.(t', x/Vs)
# 3 - Compute FK spectrum
using FFTW, Printf
Wf = abs.(fft(Wt)')
FK = 20*log10.(Wf .+ 1e-3*maximum(Wf))
FK0 = reverse(fftshift(FK[1:(Nt÷2+1),:], 2), dims=2) # physicists' FFT2D has opposite signs for time and space (outward travelling waves)
fi = rfftfreq(Nt, 1/(t[2]-t[1]))
ki = fftshift(fftfreq(Nx, 1/(x[2]-x[1])))
# 4 - Plot data
using Plots, Measures; gr(dpi=600, size=(1000,600))
p1 = heatmap(x, t, Wt', yflip=true, xlabel="x [m]", ylabel="t [s]")
p2 = heatmap(ki, fi, FK0, xlabel="k [1/m]", ylabel="f [Hz]", c=:jet)
plot(p1, p2, margins=10mm)
The aliasing is clearly seen.
Some numerical artifacts, due to the input not being periodic, would require additional work to minimize.
3 Likes
Hi @rafael.guerra can you send me the lines for this example?
The code was posted below the figure.
Sorry. I did not see! Thanks a lot!