# Implementation of dispersion compensation of lamb waves

I am trying to implement the method from Paul Wilcox paper "A rapid signal processing technique to remove the effect of dispersion from guided wave signals" on data from a lamb wave propagating in a plate (I have the dispersion curves of the plate, and can assume that I only look at A0 mode).

The signal I am working with has 501 samples, with a sample rate of 500kHz; the pulse has a center frequency of 15kHz. So it is very oversampled.

He has split the method up into steps.

• Pad the received time-trace, g(t), with zeros to yield a time domain signal containing a number of points that is an integral power of two and at least eight times as many as in the original signal.
• Perform an FFT on the zero-padded time-domain signal to obtain G(ω).
• Calculate the number of points, n, in the wavenumber domain, the size of the wavenumber step, ∆k, and the size of the distance step, ∆x, in the final distance-trace using the rules given in (12), (10), and (11), respectively.
• Interpolate G(ω) to find its value, G(k), at points equally spaced in k for the desired mode.
• Calculate the group velocity of the guided wave mode at the same wavenumber points, vgr(k).
• Compute H(k) = G(k)vgr(k).
• Apply an inverse FFT to H(k) to obtain the dispersion compensated distance-trace h(x).

I struggle to understand how to calculate ∆k, ∆x, and the new k array with equally spaced points.

I have managed to get some results by just inserting them into the formulas from the paper: $$(5): dω = vgr(ω)dk, ω = vph(ω)k$$ $$(9): n∆x>m∆t v_max.$$ $$(11): k_{Nyq} \geq k(f_{Nuq})$$ $$k_{Nyq} \geq k(1/(2\Delta t))$$ $$\frac{1}{2\Delta x} \geq k(1/(2\Delta t))$$ $$\Delta x \leq \frac{1}{2k(1/2\Delta t)}$$ $$(10): \Delta k = \frac{1}{n\Delta x} < \frac{1}{m \Delta t v_{max}}$$ $$(12): n > 2\frac{k_{Nyq}}{\Delta k}$$ Using this gives a reasonably compensated signal, but I don't know how to get the correct distance axis. What I would like my end result to be is a dispersion-compensated signal with a spatial resolution of 1mm and an axis from 0 to 1m (as the maximum wave travel length is around 0.6m). Also using a lower Nyquist frequency seems reasonable, but I am unsure of the effects of just lowering f_nyq. For example, this is how the new k array looks compared to the calculated k array, if I use the original fs and f_nyq. If I change f_nyq to 60kHz I would assume that it would be similar but I get this: I assume that there is something wrong with my implementation, or is it okay that the new k array does not match up with the calculated one?

To summarize:

• I would like some explanation on the wavenumber domain and its effect, and how the calculations of dk, dx, and new k are done.

• What I have to do if I would like to specify the max distance and the spatial resolution.

• What I have to be careful of when lowering the f_nyq.

Here is the code I am using for now if that is interesting.

def dispersion_compensation_Wilcox(file_n=2, postion=5, fs=500000):
"""
Performs dispersion compensation on the input signal.

Args:
signal (ndarray): The input time-domain signal.
postion (int): The postion we want to use.
fs (int): Samplerate

Returns:
ndarray: The dispersion compensated distance-trace.

"""

def get_k_value(freq):
v_freq = get_velocity_at_freq(freq)['A0']['phase_velocity'] #fetches the velocity at the upper frequency
k_freq = (2*np.pi*freq)/v_freq
return k_freq

wave_data_top, x_pos_top, y_pos_top, z_pos_top, time_axis_top = get_comsol_data(9) #fetches data from comsol files
wave_data_bottom, x_pos_bottom, y_pos_bottom, z_pos_bottom, time_axis_bottom = get_comsol_data(10) #fetches data from comsol files
signal = (wave_data_top[postion]+wave_data_bottom[postion])/2 #A0 mode
n_times = 16 #number of points that is an integral power of two and at least eight times as many as in the original signal.
m = len(signal)
n_fft = 2 ** int(np.ceil(np.log2(n_times * m)))
dt = 1/fs # 1/500000
freq_vel = np.fft.fftfreq(G_w.size, dt)
G_w = G_w[freq_vel>0]
freq_vel = freq_vel[freq_vel>0]
f_nyq = fs/2
v_gr, v_ph = wp.theoretical_group_phase_vel(freq_vel, material='LDPE_tonni20mm', plot=True) #group and phase velocity with the same length as freq_vel
k = (2*np.pi*freq_vel)/v_ph #same length as freq_vel. Wavenumber domain
v_nyq = get_velocity_at_freq(f_nyq)['A0']['phase_velocity'] #fetches the velocity at the nyquist frequency
v_max = np.max(v_gr)
k_nyq = get_k_value(f_nyq)
dx = 1/(2*k_nyq)
print(f'dx: {dx}')
dk = 1/(n_fft*dx)
k_max = k[-1] #doesnt matter if i use this or this 2*np.pi*upper_freq/v_max since both are equal or 2 times k_nyq
print(f'k_nyq: {k_nyq}, kmax: {k_max}')
w = 2*np.pi*freq_vel
print(f'shape of w: {w.shape}')
n = len(k)
print(f'length of k: {n}')

k_new = np.arange(0, k_max + dk, dk)

x = np.arange(0, n_fft * dx, dx)

# Interpolate G(w) to find G(k)
G_interp = interpolate.interp1d(k, G_w, kind='linear', bounds_error=False, fill_value=0)(k_new)

v_gr_interp = interpolate.interp1d(k, v_gr, kind='linear', bounds_error=False, fill_value=0)(k_new)

# Compute H(k) = G(k) * vgr(k)
H_k = G_interp * v_gr_interp

# Apply inverse FFT to H(k) to obtain the dispersion compensated distance-trace
h_x = np.fft.ifft(H_k)
print(f'shape of h_x: {h_x.shape} before removing zero-padding')
# Remove zero-padding from the compensated signal
h_x = h_x[:m]
x = np.arange(len(h_x))*dx
#normalize h_x and signal
h_x = h_x/np.max(h_x)

return h_x.real

• Are you solve the problem finally? Can you share the final code? Thanks you so much. Nov 23, 2023 at 15:44

• get_velocity_at_freq(),
• wp.theoretical_group_phase_vel()