Plate resonances spectra: Solution of Inverse Problem.
Solution of the inverse problem of thickness resonances in plates at normal incidence
Obtention of material properties from the measured spectra (magnitude and phase) of the ultrasonic transmission coefficient of a layer of the material measured at normal incidence and over a frequency range that encompases, at least, one order of the thickness resonances.
Two cases are considered:
Homogeneous plate (one layer). Gradient Descent method
Plate made of several layers (layered model). Stochastic Gradient Descent method
One layer plate
Layered plate
One Layer Model Fitting
IPython notebook (html) example
The goal of this code is to extract material properties of a plate (density, thickness, ultrasound velocity and attenuation) from the measured phase and magnitude spectra of the transmission coefficient measured at normal incidence, with the condition that one or several thickness resonances are observed. The experimental procedure and the theoretical basis that make possible to fully solve the inverse problem in this case are explained in the Ultrasonics 2010 paper (see publications). The nummerical optimization of the fitting procedure in the program below take the idea of the Ultrasonics 2010 paper to generate an initial guess and implement and optimization technique using a Gradient Descent technique.Data Repository: ultrasonic-thickness-resonance Code: download link Feedback: wiki link
# -*- coding: utf-8 -*-
"""
1 March 2014
Program to estimate material properties from the measured first order thickness resonance in the
ultrasonic transmission coefficient of a homogeneous plate at normal incidence
The following material properties are obtained:
1. Plate density
2. Plate thickness
3. Ultrasound velocity in the plate
4. Ultrasound attenuation coefficient at resonant frequency in the plate
5. Variation of the attenuation coefficient with the frequency assuming a power law.
DEPENDENCIES:
NUMPY
PYLAB (only if graphic display is enabled)
DEFINITIONS:
Measurement: FFT of the signal transmitted through the sample (magnitude -dB- and phase -rad-)
Calibration: FFT of the signal transmitted from transmitter to receiver without the sample between them
(amplitude -dB- and phase -rad-).
INPUT:
An ascii file with three columns (floats),
1st col.: Frequency (in Hz),
2nd col.: Amplitude spectrum of the measurement (in dB)
- the amplitude spectrum of the calibration (in dB).
3rd col.: Phase spectrum of the measurement (in rad).
- the phse spectrum of the calibration (in rad).
This file must be in: c:\Python27\Scripts\medida.dat
The default path is the one shown before: c:\Python27\Scripts\
Alternatively, it is possible to re-define the path to allocate input and output files in a user defined path.
This can be done by storing the user defined path in the file:
c:\Python27\Scripts\path_CT.txt
This file must contain one line of text with the desired path
(for example: c:\Users\JS\Dropbox\Python\ )
If the file: \medida.dat is no found, then he program ask for the file name, including extension, (raw_input)
ADVANCED INOUT DATA:
Other input data can be provided by creating the file: 'input_python14.txt' in c:\Python27\Scripts\
or in the user defined path.
The format of this file is one row of floats space separated. Example:
1.2000000e+000 3.4300000e+002 2.0000000e+001 4.0000000e+002 1.0000000e-002 9.9999e-001 1.0000000e+000
These data correspond to:
1. Density of the fluid where the plate is located (air at 20 C is the default case: 1.2 kg/m3)
2. Ultrasound velocity in this fluid (343m/s in the default case)
3. Loss in dB mesured from the resonance peak to limit the analysis to this frequency band (Default value: 10 dB)
4. Maximum number of steps allowed in the Gradient Descent (400 in the default case)
5. Step size in the Gradient Descent, constant (0.01 in the default case)
6. Minimum increment of the error allowed in the Gradient Descent, if no, the routine is interrupted and values returned
(0.99999 in the default case).
7. Option for graphical display 1: on, 0: off. If enabled, the the program shows plots with transmission
coefficient spectra, both, claculated and measured. This option requires PYLAB
OUTPUT FILES:
(path)...\output.dat
A file with theoretically calculated transmission coefficient spectra (with same format as the input data):
[frequency (Hz) Magnitude spectrum (dB) Phase spectrum (Hz)]
(path)...\output2.dat
A file with 3 rows and 6 columns:
[1:plate thickness (m); 2:ultrasound velocity(m/s); 3:attenuation coeff./Np/m); 4.plate density (kg/m3); 5.exponent in the attenuation vs freq power law; 6.resonant frequency (Hz)]
[1:IEE Magnitude spectrum; 2: IEE Phase Spectrum; 3. FEE Magnitude spectrum; 4. FEE Phase spectrum; 5. Number of steps in the GD routine; 6. Loss from peak for bandwidth (dB)]
[1:Q-factor of the resonance; 2:Maximum of the Magnitude spectrum (dB); 3:Running time (s) ; 4: Maximum number of steps allowed in GD routine; 5: Step size in GD routine; 6: Minimum error improvement for another step in GD]
IEE: Initial error estiamtion
FEE: Final error estimation
CALCULATION METHOD:
Initial values are obtained from the resonant frequency, the amplitude and the phase at resonance,
and the Q factor of the resonance. The layer parameters are modified to improve the fitting of the
theoretically predicted spectra into the experimentally measured one. The minimization of the error
is performed by using a Gradient Descent algorithm.
Author: Tomas Gómez Alvarez-Arenas
ITEFI-CSIC
www.us-biomat.com
t.gomez@csic.es
03/25/2014
Starting point is CT12_GD.py
"""
import time
start_time = time.time()
import numpy as np
class CoefTrans(object):
"""
Transmission Coefficient class: sound transmission through a layer at normal incidence
3 columns of floats
Format: [frequency (Hz); magnitud spectrum (dB); phase spectrum (rad)] """
def __init__(self,CT):
self.CT = CT
def __str__(self):
return self.CT
def freq(self):
return self.CT[:,0]
def mag(self):
return self.CT[:,1]
def phase(self):
return self.CT[:,2]
def errorM(self,other):
""" Difference (square minima) between measured and calculated magnitude spectrum """
return (sum(((self.mag() - other.mag()) / self.mag()) ** 2)) ** (1.0/2.0) / len(self.CT[:,1]) * 100
def errorPh2(self,other):
return (sum(((self.phase() - other.phase()) ) ** 2)) ** (1.0/2.0) / len(self.CT[:,2]) * 100
def fres(self):
""" Fres function applies to a Transmission Coefficient
and returns four parameters: Resonant frequency (Hz); Magnitud spectrum value at resonance (dB); Phase spectrum value at resonance (rad)
index of the resonant frequency (in the frequency vector """
n = 0
while not (self.mag()[n] == max(self.mag()) ):
n += 1
return [self.freq()[n], self.mag()[n], self.phase()[n], n]
def Q(self):
""" This function applies to a transmission coefficient and returns Q factor and the index of the
frequency vector that corresponds to the masimum value of the amplitude spectrum (resonant frequency) """
ind = self.fres()[3]
indm = ind
indM = ind
limit1 = 0
limit2 = 0
while self.mag()[ind] - self.mag()[indm] < 3:
indm -= 1
if indm == 0:
limit1 = 1
break
while self.mag()[ind] - self.mag()[indM] < 3:
indM += 1
if indM == len(self.mag()):
limit2 = 1
break
if limit1 + limit2 == 0:
return [(self.freq()[indM] - self.freq()[indm]) / self.fres()[0], (indm + indM) / 2]
if limit1 == 1:
return [2 * (self.freq()[indM] - self.fres()[0]) / self.fres()[0], self.fres()[3]]
if limit2 == 1:
return [2 * (self.fres()[0] - self.freq()[indm]) / self.fres()[0], self.fres()[3]]
if limit1 + limit2 == 2:
raise NameError('There is no resonance')
def ParamPunt(self,other):
""" Other es un medium
calcula los parametros de la capa a partir de la medida de coef Trans
en particular, emplea los valores obtenidos a fres """
tiempo = self.fres()[2] / (2 * 3.14159 * self.fres()[0])
espesor = other.v_med * (tiempo + 1 / (2 * self.fres()[0]))
velocidad = 2 * espesor * self.fres()[0]
z1 = other.v_med * other.rho_med
att = 3.1415926549 * self.fres()[0] / ( velocidad / self.Q()[0] )
c = attyRho(att, espesor, 10 ** (self.fres()[1] /20))
att = c[0]
rho = z1 / c[1] / velocidad
return [espesor, velocidad, att, rho, 1, self.fres()[0]]
def ResWindow(self, BOUND):
""" This function applies to a Transmission Coefficient and takes the argument BOUND which is a negative integer.
in (dB). It returns a transmission coefficient that is the one taken as argument but limited to the frequency band
that correspond to Maximum value of transmission coefficient magnitude -BOUND. """
ind = self.fres()[3]
indm = ind
indM = ind
limit1 = 0
limit2 = 0
while self.mag()[ind] - self.mag()[indm] < BOUND:
indm -= 1
if indm == 0:
limit1 = 1
break
while self.mag()[ind] - self.mag()[indM] < BOUND:
indM += 1
if indM == len(self.mag()):
limit2 = 1
break
teoarr=np.vstack((self.freq()[indm:indM], self.mag()[indm:indM], self.phase()[indm:indM]))
return CoefTrans(teoarr.T)
class medium(object):
def __init__(self,rho_med,v_med):
self.rho_med = rho_med
self.v_med = v_med
def __str__(self):
return 'Density (kg/m3): ' + str(self.rho_med) + '; Velocity (m/s): ' + str(self.v_med)
class layer(object):
def __init__(self,lay):
self.lay = lay
def thick_layer(self):
return self.lay[0]
def v_layer(self):
return self.lay[1]
def att_layer(self):
return self.lay[2]
def rho_layer(self):
return self.lay[3]
def efa_layer(self):
return self.lay[4]
def frec0_layer(self):
return self.lay[5]
def __str__(self):
return 'Thickness (um): ' + str(self.thick_layer()*1e6) + '; Velocity (m/s): ' + str(self.v_layer()) + '; Attenuation @ resonant frequency(Np/m): ' + str(self.att_layer()) + '; Density (kg/m3): ' + str(self.rho_layer()) + '; EFA: ' + str(self.efa_layer()) + '; Resonant frequency (Hz): ' + str(self.frec0_layer())
def attyRho(alf,esp,Tmax):
""" Calculates attenuation and impedance ratio from initial parameter estimation """
r = np.sin(alf * esp) / 2 * Tmax
R = (( 1 - r ) ** 2 ) / (( 1 + r ) ** 2 )
alf = alf - np.log(R) / (2 * esp)
return [alf, r]
def cT_teo(medium,cT_exp):
z1 = medium.rho_med * medium.v_med
freq = cT_exp.freq()
j = complex(0, 1)
k = 2 * 3.14159 * freq / cT_exp.ParamPunt(medium)[1] + j * cT_exp.ParamPunt(medium)[2] * ((freq / cT_exp.fres()[0]) ** (cT_exp.ParamPunt(medium)[4] ))
veloc = 2 * 3.141592 * freq / k
z2 = veloc * cT_exp.ParamPunt(medium)[3]
espesor = cT_exp.ParamPunt(medium)[0]
T = - 2.0 * z1 * z2 / (2 * z1 * z2 * np.cos(k * espesor) + j * (z1 ** 2 + z2 **2) * np.sin(k * espesor))
Tm = 20 * np.log10( abs (T))
Tphas = np.unwrap(2 * 3.1415926 * freq * espesor / medium.v_med - np.angle(T))
teoarr=np.vstack((freq, Tm, Tphas))
return CoefTrans(teoarr.T)
def cT_layer(medium,layer,freq):
z1 = medium.rho_med * medium.v_med
j = complex(0, 1)
k = 2 * 3.14159 * freq / layer.v_layer() + j * layer.att_layer() * (freq / layer.frec0_layer()) ** layer.efa_layer()
veloc = 2 * 3.141592 * freq / k
z2 = veloc * layer.rho_layer()
T = - 2.0 * z1 * z2 / (2 * z1 * z2 * np.cos(k * layer.thick_layer()) + j * (z1 ** 2 + z2 **2) * np.sin(k * layer.thick_layer()))
Tm = 20 * np.log10( abs (T))
Tphas = np.unwrap(2 * 3.1415926 * freq * layer.thick_layer() / medium.v_med - np.angle(T))
teoarr=np.vstack((freq, Tm, Tphas))
return CoefTrans(teoarr.T)
def fitt_EFA2(medium, cT_exp, layer0):
Flag = 0
cTlayer0 = cT_layer(medium, layer0, cT_exp.freq())
ErrorM0 = cTlayer0.errorM(cT_exp)
ErrorPh0 = cTlayer0.errorPh2(cT_exp)
for n in range(0,20,1):
layerT = layer([layer0.thick_layer(), layer0.v_layer(), layer0.att_layer(), layer0.rho_layer(), n/10.0, layer0.frec0_layer()])
cTlayerT = cT_layer(medium, layerT, cT_exp.freq())
ErrorM = cTlayerT.errorM(cT_exp)
ErrorPh = cTlayerT.errorPh2(cT_exp)
if ErrorM + ErrorPh < ErrorM0 + ErrorPh0: #and ErrorPh < ErrorPh0:
Flag = 1
layerMin = layerT
if Flag == 1:
return layerMin
else:
return layer0
def F_Error(medium,cT_exp,layer0):
"""" Computes the error of a given calculated transmission coefficient spectra and the measured one """
cTlayer0 = cT_layer(medium, layer0, cT_exp.freq())
Error = [cTlayer0.errorM(cT_exp), cTlayer0.errorPh2(cT_exp)]
return Error
def Grad_Error(medium,cT_exp,layer0,delta):
""" Computes the local gradient of the error function"""
Error0 = F_Error(medium,cT_exp,layer0)
grad = [1 - delta, 1.0, 1 + delta]
layer_list = [layer0.thick_layer(), layer0.v_layer(), layer0.att_layer(), layer0.rho_layer(), layer0.efa_layer(), layer0.frec0_layer()]
layer_fin = layer_list[:]
layer_test = layer_list[:]
for step1 in grad:
layer_test[0] = layer_list[0] * step1
for step2 in grad:
layer_test[1] = layer_list[1] * step2
for step3 in grad:
layer_test[2] = layer_list[2] * step3
for step4 in grad:
layer_test[3] = layer_list[3] * step4
for step5 in grad:
layer_test[4] = layer_list[4] * step5
layerT = layer(layer_test)
Error = F_Error(medium,cT_exp,layerT)
if Error[0] < Error0[0] and Error[1] < Error0[1]:
layer_fin[:] = layer_test[:]
Error0 = Error[:]
return layer(layer_fin)
def Walk_downGradient(medium,cT_exp,layer0,delta,MaxSteps,Prec):
""" This function takes one step in the fitting procedure douwn the gradient in the error hiperspace """
max = 0
layerF = Grad_Error(medium,cT_exp,layer0,delta)
Error0 = F_Error(medium,cT_exp,layer0)
Error = F_Error(medium,cT_exp,layerF)
while Error[0] < Error0[0] * Prec and Error[1] < Error0[1] * Prec: layer0 = layerF layerF = Grad_Error(medium,cT_exp,layer0,delta) #print layerF Error0 = F_Error(medium,cT_exp,layer0) Error = F_Error(medium,cT_exp,layerF) max += 1 if max > MaxSteps:
break
return [layer0.thick_layer(), layer0.v_layer(), layer0.att_layer(), layer0.rho_layer(), layer0.efa_layer(), layer0.frec0_layer(), max]
try:
fpath = open('c:\Python27\Scripts\path_CT.txt')
for line in fpath:
path = line
fpath.close()
except IOError:
path = 'c:\\Python27\\Scripts\\'
try:
data = np.loadtxt(path + 'input_python14.txt')
except IOError:
data = np.array([1.204, 343, 10, 400, 0.01, 0.99999, 1])
#default input data:
#surrounding fluid: air @ 20 C
#resonance bandwidth: -20 dB of the peak value
#Maximum 400 steps in GD
#Fixed step size: 0.01
#Minimum error improvement
#Graphic display 1: yes; 0: no
caso1 = medium(data[0], data[1])
LIMITE = data[2]
MaxSteps = data[3]
delta = data[4]
Prec = data[5]
GraphOutput = data[6]
try:
test0 = np.loadtxt(path + 'test23.dat')
except IOError:
FileName = raw_input('Name of the file with the transmission coefficient spectra, include extension, (default path is c:\Python27\Scripts\) ')
test0 = np.loadtxt(path + FileName)
test1 = test0[np.argsort(test0[:,0])]
cttest1 = CoefTrans(test1)
aa = layer(cttest1.ParamPunt(caso1))
fit = Walk_downGradient(caso1,cttest1.ResWindow(LIMITE),aa,delta,MaxSteps,Prec)
fitt = layer(fit[0:6])
lfinEFA = fitt_EFA2(caso1, cttest1, fitt )
cT_fin = cT_layer(caso1, lfinEFA, cttest1.freq())
f0 = open(path + 'output.dat','w')
f2 = open(path + 'output2.dat','w')
ind = 0
for FR in cT_fin.freq():
linea = str(cT_fin.freq()[ind].real) + ' ' + str(cT_fin.mag()[ind].real) + ' ' + str(cT_fin.phase()[ind].real) +'\n'
f0.write(linea)
ind += 1
ERROR_0 = F_Error(caso1,cttest1.ResWindow(LIMITE),aa)
ERROR_FIN = F_Error(caso1,cttest1.ResWindow(LIMITE),lfinEFA)
etime = time.time() - start_time
data0 = str(lfinEFA.thick_layer().real) + ' ' + str(lfinEFA.v_layer().real) + ' ' + str(lfinEFA.att_layer().real) + ' ' + str(lfinEFA.rho_layer().real) + ' ' + str(lfinEFA.efa_layer().real) + ' ' + str(lfinEFA.frec0_layer().real) +'\n'
data1 = str(ERROR_0[0]) + ' ' + str(ERROR_0[1]) + ' ' + str(ERROR_FIN[0]) + ' ' + str(ERROR_FIN[1]) + ' ' + str(fit[6]) + ' ' + str(LIMITE) +'\n'
data2 = str(1.0/cttest1.Q()[0]) + ' ' + str(max(cttest1.mag())) + ' ' + str(etime) + ' ' + str(MaxSteps) + ' ' + str(delta) + ' ' + str(Prec)
f2.write(data0 + data1 + data2)
f0.close()
f2.close()
if GraphOutput == 1:
cT_band = cT_layer(caso1, lfinEFA, cttest1.ResWindow(LIMITE).freq())
import pylab
pylab.figure(1)
pylab.subplot(211)
pylab.title('Transmission coeff. Results in the fitting frequency band -' + str(LIMITE)+'dB')
pylab.xlabel('Frequency (MHz)')
pylab.ylabel('TC Magnitude Spectrum (dB)')
pylab.plot(cT_band.freq()/1e6, cT_band.mag(),cttest1.ResWindow(LIMITE).freq()/1e6,cttest1.ResWindow(LIMITE).mag(),'o')
pylab.subplot(212)
pylab.xlabel('Frequency (MHz)')
pylab.ylabel('CT Phase Spectrum (rad)')
pylab.plot(cT_band.freq()/1e6, cT_band.phase(),cttest1.ResWindow(LIMITE).freq()/1e6,cttest1.ResWindow(LIMITE).phase(),'o')
pylab.figure(2)
pylab.subplot(211)
pylab.title('Transmission coeff. Results in the whole frequency range')
pylab.xlabel('Frequency (MHz)')
pylab.ylabel('TC Magnitude Spectrum (dB)')
pylab.plot(cT_fin.freq()/1e6, cT_fin.mag(),cttest1.freq()/1e6,cttest1.mag(),'o')
pylab.subplot(212)
pylab.xlabel('Frequency (MHz)')
pylab.ylabel('CT Phase Spectrum (rad)')
pylab.plot(cT_fin.freq()/1e6, cT_fin.phase(),cttest1.freq()/1e6,cttest1.phase(),'o')
pylab.show()
Ultrasounds for Biological Applications and Materials Science 