`ffht` does not handle matrix form of the input

[sgkang]

Hi @prisae, I really like the idea of performing hankel transformation at once rather than looping over for frequency. So, for the time to frequency transform I want to do the same procedure, particularly when I am computing sensitivity. Hence, the input frequency response can have size (n_frequency x n_layer).

Following snippet of the code emulates the situation that I want to perform, but that outputs errors. Have you tried this?

from empymod import filters, transform, utils
import numpy as np

fftfilt = filters.key_81_CosSin_2009()
time = np.logspace(-5, -2, 21)
n_layer = 19
time_out, freq, ft, ftarg = utils.check_time(time, 0, 'sin', {'pts_per_dec':2, 'fftfilt': fftfilt}, 0)
a =1e2
fEM =1j*freq*2*np.pi/(a + 1j*freq*2*np.pi )
signal = fEM/(2j*np.pi*freq)
resp, _ = transform.ffht(np.tile(signal.reshape([-1, 1]), (1, n_layer)), time, freq, ftarg) 
resp *= 2/np.pi

[prisae]

@sgkang, transform.ffht/transform.dlf take care of the tiling/reshaping. Just change resp, _ = transform.ffht(np.tile(signal.reshape([-1, 1]), (1, n_layer)), time, freq, ftarg) to resp, _ = transform.ffht(signal, time, freq, ftarg) and you should be good to go.

Re-open it if it doesn't resolve the issue.

[prisae]

By the way, how reshaping is done depends on pts_per_dec, if it is smaller than, equal to, or greater than zero. Might be worth checking if pts_per_dec:-1 (Lagged Convolution DLF) is not faster than pts_per_dec:2 (splined DLF). It will most likely be more precise, because two points per decade for spline is pretty low.

[sgkang]

Hmm... not sure we are in the same page. Let's suppose I have multipe receiver measuring signal in time domain. For this case my frequency domain reponse first need to be computed: response (nfreq, nrx). Currently, I can input a column of the response (nfreq,1) to ffht, but I was wondering that I can put all of my response (nfreq, nrx). Is t possible? Or I just need to loop over for frequency?

Then now I want to convert them into time.

[prisae]

Ah. You need to loop over offsets. In empymod.model.tem (lines 1551 & 1552):

    for i in range(off.size):
        out = getattr(transform, ft)(fEM[:, i]*fact, time, freq, ftarg)

where getattr(transform, ft) in your case translates to transform.ffht.

I don't think there is a way to speed that up. The reason you can do it for multiple offsets in the Hankel-transform is that you go from wavenumber-to-offsets, so offset is part of your transform and you can do the tricks there. However, in the Fourier transform you go from frequency-to-time, so offset does not play a part in that, and each offset needs its own transform. So you would basically need a 2D DLF, but I doubt you would save time over a loop.

[sgkang]

I see. Got it now. Thanks @prisae!

[prisae]

@sgkang, just a follow-up comment: If you have a LOT of offsets, you could potentially do some interpolation along the offsets: only calculate certain offsets, and interpolate the responses in-between. I think the time-domain responses do not only vary smoothly along the time-axis, but also along the offset-axis.