Ultrasounds for Biological Applications and Materials Science

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

One layer plate

 

two_layers

Layered plate

 

 

One Layer Model Fitting

IPython notebook (html) example

Shared project on github


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.

github-button
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()

Advertisements
0034 915618806-058
t.gomez@csic.es
Skype: usbiomat_csic
wordpress stats plugin