Attenuation and Dispersion correction

Hello,

I am currently working on correcting the attenuation and dispersion in my measurement signal, but I am not getting the desired result.

Briefly about the measurement signal: These are pressure signals, some of which pass through a material and are measured at a piezo sensor at time t1 = 1 us and t2 = 1.45 us. The pressure signal (signal1) at time t1 can be assumed to be undamped and non-dispersive, as this is formed directly in front of the piezo sensor and is measured without passing through the material. For the pressure signal (signal2), which arrives at the sensor at time t2, an attenuation and dispersion correction must be carried out, as this signal2 has passed through the material, but originally also had the same absolute value and pulse width as signal1 (just with a different sign).

I performed the calculation using the usual formulae, but did not get the expected and desired correction.

I have the feeling that it's not such a big problem, but I can't find the solution. Perhaps one of you will spot my mistake(s) or have a new or different approach.

I am grateful for any help.

%% Correction of Dispersion and Attenuation
clear;
clc;
close all;

dt = 5e-10;
fs = 1/dt;
time = 0:dt:2.000500000000000e-06;

%Variables for gaussian-functions

GRD_Ampl = -1.02251282295196;
HS_Ampl = 0.497967501946497;
a1 = 6.28704925677956e+15;
a2 = 3.38523088091301e+14;
GRD_time = 2e-07; 
HS_time = 1.453e-06;


% Variables of the specimen

d = 0.002;
v = 1600;


% Building of the gaussian-functions

Ideal_GRD = GRD_Ampl*exp(-a1 * (time - GRD_time).^2);
Ideal_HS = HS_Ampl*exp(-a2 * (time - HS_time).^2);

Signal = Ideal_HS + Ideal_GRD;

% Fourier transformation of the pressure signals

freq = fs/length(Signal)*(-length(Signal)/2+1:length(Signal)/2);

F_Ideal_GRD = fft(Ideal_GRD);
F_Ideal_HS = fft(Ideal_HS);
F_Signal = fft(Signal);

% Calculation of the attenuation and dispersion factors

alpha = -(1/d)*log(abs(F_Ideal_HS./F_Ideal_GRD)); % attenuation factor
beta = -(1/d)*(phase(F_Ideal_HS)-phase(F_Ideal_GRD)); % dispersion factor

figure;
tiledlayout(2,1);
nexttile(1)
plot(freq, fftshift(alpha), 'LineWidth', 2);
title('attenuation factor alpha');
xlabel('freqency');
ylabel('Amplitude');
grid on;

nexttile(2)
plot(freq, fftshift(beta), 'LineWidth', 2);
title('dispersion factor beta');
xlabel('freqency');
ylabel('Amplitude');
grid on;

% Calculation of the corrected signal

G = exp((-alpha*d)+(-1i*beta*d));

F_Signal_corrected = F_Signal.*G;
Signal_corrected = ifft(F_Signal_corrected);

figure;
plot(time, Signal_corrected, 'LineWidth', 2);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

title('Ideal corrected pressure pulse signal');
xlabel('Time in s');
ylabel('Amplitude');
grid on

Answer 1

[Sanju]

Hi @Christian Kohlsaat,

I understand that you are trying to correct the attenuation and dispersion in your measurement signal. However, you mentioned that you are not getting the desired result. It would be helpful if you could provide more information about the expected and actual results, as well as any specific issues or error messages you encountered.

In the provided code, you are calculating the attenuation and dispersion factors using the Fourier transformation of the pressure signals. Then, you are applying these factors to correct the signal. One potential issue could be calculation of the dispersion factor, as the phase difference between the Fourier transforms might not be accurate. You can try using the “unwrap” function to improve the accuracy of the phase calculation.

Here is an example of how you can use the unwrap function to calculate the dispersion factor:

beta = -(1/d)*unwrap(phase(F_Ideal_HS) - phase(F_Ideal_GRD));

Additionally, there’s a warning about ignoring the imaginary parts of the complex arguments when plotting the corrected signal. This warning suggests that the corrected signal might have complex values. To visualize the real part of the corrected signal, you can modify the plotting code as follows:

plot(time, real(Signal_corrected), 'LineWidth', 2);  

You can refer to the below documentation on “unwrap” function if required,

https://www.mathworks.com/help/matlab/ref/unwrap.html

Hope this Helps!

Thanks.