1
$\begingroup$

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. enter image description here If I change f_nyq to 60kHz I would assume that it would be similar but I get this: enter image description here 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)))
    signal_padded = np.pad(signal, (0, n_fft - m), mode='constant')
    #computing fft of padded signal
    G_w = np.fft.fft(signal_padded)
    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
$\endgroup$
1
  • $\begingroup$ Are you solve the problem finally? Can you share the final code? Thanks you so much. $\endgroup$ Commented Nov 23, 2023 at 15:44

1 Answer 1

1
$\begingroup$

In order to help you I need you to expand the functions:

  • get_velocity_at_freq(),
  • wp.theoretical_group_phase_vel()

that are not detailed in your code.

$\endgroup$
1
  • $\begingroup$ Hello and welcome to Signal Processing SE. Please note that this is a Q&A site and not a forum, so please refrain from posting comments as answers. I understand that you don’t enough “reputation” yet to leave comments on other people’s posts, but if you get involved just a little bit (answer a couple of questions, post another, etc.) this “privilege” will be acquired quite fast. Until then, please do not post comments as answers as they will most probably deleted. $\endgroup$
    – ZaellixA
    Commented Nov 30, 2023 at 9:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.