• No results found

A practical two-transducer ultrasonic flow meter

N/A
N/A
Protected

Academic year: 2021

Share "A practical two-transducer ultrasonic flow meter"

Copied!
72
0
0

Loading.... (view fulltext now)

Full text

(1)2008:102. MASTER'S THESIS. A Practical Two-Transducer Ultrasonic Flow Meter. Amin.G. Saremi. Luleå University of Technology Master Thesis, Continuation Courses Electrical engineering Department of Computer Science and Electrical Engineering Division of EISLAB 2008:102 - ISSN: 1653-0187 - ISRN: LTU-PB-EX--08/102--SE.

(2) 2.

(3) Preface Nowadays ultrasonic transducers are extensively used. Many industrial systems carry out their measurement process by taking the advantage of ultrasound waves. Sound of frequencies above audible range is called ultrasound, ultrasonic waves are classified under category of mechanical waves. Likewise all other mechanical waves, ultrasound needs material environment called medium to propagate through by moving the particles. Since the end of seventies, ultrasonic flow meters have been widely designed and applied in different industrial plans mainly related to fluid mechanics, chemistry, oil industry, etc. However, the initial motive of this thesis was to design a reliable mass flow meter applicable in controling the precise amount of energy transferred through gas pipe lines. The mass flow meter presented in this thesis is based on a two-transducer ultrasonic measurement set up where both transducers emit ultrasonic beam simultaneously. In this thesis we will describe how by checking the amplitudes and timings of the received echoes many valuable physical properties of the flowing liquid can be revealed. In other words, we are able to calculate the speed of sound in the liquid, the flow velocity as well as the liquid’s acoustic impedance and its density by only investigating the amplitudes and timings of the received ultrasonic echoes. Moreover, we show that frequency behaviour of the system and the liquid can be illuminated thoroughly. As a master thesis, this project partially stands on previous works in this field and mainly on the works done by Dr. Jan van Deventer and Dr. Jonny Johansson and various scientific publications. In coming discussions not only do we systematically take advantage of previous thoughts, designs and results in order to explore further achievements, but also we present many new solutions. This thesis is aiming towards raising new electronic solutions and signal processing techniques to reduce the uncertainty of the estimations and to speed up the computation.. Acknowlegements Firstly, I’d like to appreciate my lovely parents for all the sincere sacrifices they made all through their lives in order to provide me with good educational chances. Secondly, I must thank the nice staff at department of computer science and electrical engineering, Lule˚ a University of technology, Sweden. Speciall thanks to Dr. Jan van Deventer, Dr. Johan Carlson, Dr. Jonny Johansson and all those who generously helped the author accomplish this project.. i.

(4) ii.

(5) Contents Preface 0.1. Acknowlegements . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. i i 1. 1 Physical and electrical properties of ultrasonic transducers 3 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 A glance at the acoustic world . . . . . . . . . . . . . . . . . . . 3 1.2.1 Acoustic impedance and reflection . . . . . . . . . . . . . 3 1.2.2 Attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3 Diffraction . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.4 Frequency dependence . . . . . . . . . . . . . . . . . . . . 4 1.3 The piezoelectric transducer . . . . . . . . . . . . . . . . . . . . . 5 1.3.1 Mechanical-Electrical properties . . . . . . . . . . . . . . 5 1.3.2 Ultrasonic transducers versus the microphones . . . . . . 6 1.3.3 Sensor design . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Modeling of the ultrasound system . . . . . . . . . . . . . . . . . 7 1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2. 2.1 2.2 2.3. 2.4 2.5. 2.6. The experimental set up 11 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 The physical scheme of the measurement system . . . . . . . . . 11 The electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1 Exciting electronic circuit . . . . . . . . . . . . . . . . . . 12 2.3.2 Listening electronic circuit . . . . . . . . . . . . . . . . . . 14 2.3.3 Combined transmission and reception . . . . . . . . . . . 14 2.3.4 ASIC: application specific integrated circuit . . . . . . . . 14 A/D conversion and the list of the measured signals in this project 16 Calculation of the flow velocity and the sound speed in the flow . 17 2.5.1 Calculation of the flow velocity . . . . . . . . . . . . . . . 18 2.5.2 Calculation of the sound speed . . . . . . . . . . . . . . . 19 Calculation of the acoustic impedance, density and mass flow rate 20 2.6.1 Reflection coefficients and Calculation of Acoustic impedance 21 2.6.2 Calculation of the density . . . . . . . . . . . . . . . . . . 23 2.6.3 Calculation of the mass flow rate . . . . . . . . . . . . . . 23 iii.

(6) iv 3. 3.1 3.2. 3.3. 3.4. 3.5. 3.6. Signal processing solutions to estimate the flight times 25 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Block diagram of the signal processing unit . . . . . . . . . . . . 25 3.2.1 Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.2 The methods . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.3 Compute and display . . . . . . . . . . . . . . . . . . . . 29 Method 1: Cross correlation solution . . . . . . . . . . . . . . . . 30 3.3.1 Time-domain cross correlation . . . . . . . . . . . . . . . 30 3.3.2 Frequency domain Cross correlation . . . . . . . . . . . . 31 3.3.3 Frequency domain versus Time domain . . . . . . . . . . 31 3.3.4 Sub-sample precision . . . . . . . . . . . . . . . . . . . . . 32 Method 2: Optimization solution . . . . . . . . . . . . . . . . . . 33 3.4.1 The ideal model of the ultrasonic signals . . . . . . . . . . 33 3.4.2 Generation of the ideal model in MATLAB . . . . . . . . 34 3.4.3 Least mean square algorithm . . . . . . . . . . . . . . . . 35 3.4.4 Performance, Advantages and drawbacks of the optimization solution . . . . . . . . . . . . . . . . . . . . . . . . . 36 Method 3: Phase Unwrapping solution . . . . . . . . . . . . . . . 36 3.5.1 Frequency response of non-ideal ultrasonic system . . . . 36 3.5.2 Phase unwrapping . . . . . . . . . . . . . . . . . . . . . . 37 3.5.3 Advantages, drawbacks and corrections . . . . . . . . . . 38 Results and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 39. 4. 4.1 4.2. 4.3. 4.4. frequency investigation of the system 41 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Calculation of the frequency response . . . . . . . . . . . . . . . 44 4.2.1 Computation of RP L , the reflection factor . . . . . . . . . 44 4.2.2 Simulation in MATLAB . . . . . . . . . . . . . . . . . . . 44 Valuable information obtained from the frequency response . . . 46 4.3.1 Amplitude of the frequency response . . . . . . . . . . . . 46 4.3.2 Phase of the frequency response and the flight time . . . . 46 4.3.3 Non-linearity of the phase curve . . . . . . . . . . . . . . 46 4.3.4 Frequency dependence of the sound speed in the liquid . . 47 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48. 5. 5.1 5.2 5.3. 5.4. Density estimation and the mass flow meter 49 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Estimation of the acoustic impedance of the liquid . . . . . . . . 49 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.3.1 Why do the probes give out different results? . . . . . . . 51 5.3.2 Averaging the results from each of the probes . . . . . . . 52 5.3.3 Computaion of the flow rate, m ˙ . . . . . . . . . . . . . . . 53 Electronic impementation . . . . . . . . . . . . . . . . . . . . . . 53.

(7) v 5.5. Advatages, drawbacks and conclusions . . . . . . . . . . . . . . .. 54. 6.1 6.2 6.3. Results, Conclusion and future works 57 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58. 6.

(8) vi.

(9) List of Figures 1.1 1.2 1.3 1.4 1.5. 1.6 2.1. 2.2 2.3 2.4. 2.5 2.6 2.7. 2.8 2.9 3.1 3.2. Diffraction effect . . . . . . . . . . . . . . . . . . . . . . . . . . . Side view of a single element transducer. The arrows show the displacement movement of the transducer and the cap. . . . . . . Capacitive microphone, the moving plate is displaced due to strength of the existing acoustic pressure field . . . . . . . . . . . . . . . A piezoelectric transducer made up in a plastic container and together with the backing layer . . . . . . . . . . . . . . . . . . . Echo signal produced by a transducer with poorly matched backing (a) and with well-matched backing (b). Note that the ringing swings have been canceled in (b) however the amplitude of the passed pressure field is lower than in (a) . . . . . . . . . . . . . the Leach equivalent circuit of a piezoelectric device . . . . . . . Cross section of the system. The left probe is complete whilst the right probe is a cross section view. The liquid flows through the inlet-outlet pipe along the path BC (or reversely CB) . . . . . . . The outer look of the system. . . . . . . . . . . . . . . . . . . . . A popular push-pull system for generating the square driving pulses [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Block shcematic of the complete chip outlining the main blocks and external connections. Control signal inputs can be bonded directly on the chip. [2] . . . . . . . . . . . . . . . . . . . . . . . Chip functionality illustrated by the behaviour of the voltage on the transducer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . The useful signals in our samples numbered according to the text . Cross section of the system. The left probe is complete whilst the right probe is a cross section view. The liquid flows through the inlet-outlet pipe along the path BC (or reversely CB) It’s obvious that if the other direction is the case then all our results will be simply multiplied by a minus sign. . . . . . . . . . . . . . . . . . The reflection of the interface between environment1 and 2. . . . One of the probes, the steps inside it and the probe-liquid interface The diagram of the signal processing unit . . . . . . . . . . . . . Pre-processed samples obtained by sampling the measured signals at each of the channels (transducers) . . . . . . . . . . . . . . . .. vii. 4 5 6 8. 8 9. 12 13 14. 15 16 17. 19 21 23 26 26.

(10) viii 3.3. 3.4. 3.5. 3.6 3.7. 3.8. 3.9 3.10. 3.11 3.12. 3.13 3.14. 4.1. 4.2. 4.3 4.4. The Amplitude of the Fourier transform (a representetive of Power Spectrum Density (PSD)) of each of the received signals by transducer(a) and by transducer(b). This shows how strong are the low frequency terms. our favourite ultrasonic data exits around 2Mhz The amplitude of the Fourier transform of each of the measurements being zoomed around 2Mhz versus our designed 6th order FIR filter. Our band pass filter is with central frequency fc = 2.062M hz and -3db cut off frequencies of fl = 1.9375M hz and fh = 2.1875M hz.) . . . . . . . . . . . . . . . . . . . . . . . . The output of the bandpass filter seen in figure 3.4. The unwanted terms are cancelled and the ultrasonic pulses look clearly outstanding. Compare this figure with figure 2.6 to see how we name each of the pulses.) . . . . . . . . . . . . . . . . . . . . . . The extracted (gated) ultrasonic pulses , e1a , e1b , e2a , e2b , e3a , e3b , e4a , e4b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The result of time domain cross correlation of e4a with e4b . The maximum represents the number of shifts by which the sequences have the match with each other. Note that the length of the resulting correlation of x and y is (length of (x))+(length of (y)) 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The result of frequency domain cross correlation of e4a with e4b . The maximum represents the number of shifts by which the sequences have the match with each other. Note that it is identical to previous figure, the only difference is it is performed 126 times faster! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The impulse response of the ideal model defined in equation (3.4) the equivalent of figure 3.8 in frequency domain, H(ω) is called the frequency response and is the Fourier transform of the impulse response, h[n] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x[n] versus y[n] where y is shifted 344.86 samples and attenuated 0.75 times. In other words, y[n] = 0.75 x[n − 344.86] . . . . . . . The initial guesses are τ = 310 and A = 1. The optimization algorithm converges to tau=344.86 and A=0.75. The real y[n] versus the estimated one by optimization solution. They perfectly match! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The blok diagram of our non-ideal system. h[n] is the impulse response whilst H(jω) is called the frequency response of the system The time of flight for ideally-modeled x and y in samples, which is precisely true . . . . . . . . . . . . . . . . . . . . . . . . . . . . The block diagram of our ultrasonic system (shown in fig2.1) in frequency domain when transducer(a) is being excited and transducer(b) listening . . . . . . . . . . . . . . . . . . . . . . . . . . . The block diagram of our ultrasonic system (shown in fig2.1) in frequency domain when transducer(b) is being excited and transducer(a) listening . . . . . . . . . . . . . . . . . . . . . . . . . . . The beam while being reflected from the probe-liquid interface . . Frequency response of the liquid for the first configuration shown in figure4.2, the frequency response conveys very valuable information which will is discussed down here. . . . . . . . . . . . . .. 27. 28. 29 30. 32. 33 34. 34 35. 37 37 39. 42. 42 44. 45.

(11) ix 4.5. 4.6 5.1 5.2 5.3 5.4 5.5. Due to the figure 4.4, the liner approximation of the frequency response of the liquid will look like this. The attenuation factor is 0.13 and the delay is 1.968e-6 [s]. However the true frequency response is the one in figure 4.4 in blue. . . . . . . . . . . . . . Frequency dependance of the sound speed . . . . . . . . . . . . . . The density estimation by data from probe(a) versus probe(b) . . The Output of the Schmitt trigger circuit . . . . . . . . . . . . . The estimated density of water, the red line is the average of the results from the probes which looks reasonable for water at 23◦C. The electronic schematic of the designed densitometer . . . . . . The output of the amplifier stage . . . . . . . . . . . . . . . . . .. 46 47 51 52 52 53 54.

(12) x.

(13) 1. 0.1. Introduction. One can think of many applications for ultrasonic systems in modern industry, ranging from toothpaste manufacturing to space engineering. However the initial motivation for starting current project was the application of an ultrasonic flowmeter in gas industry. In harsh winters of Lapland the temperature easily reaches -30◦ C. Some days in summer on the other hand, the temperature exceeds +30◦ C. This enormous temperature difference has a drastic influence on density of oil products and consequently on oil pricing. In other words, filling our car tanks at gas statations in summers costs more than doing it during winters, for the same amount of energy! If we supply the gas stations with flow meters we will be able to pay for the actual mass of gas instead of paying for the temperature-dependant volume of it. Note that it is the mass of the gas which acounts for the absolute amount of energy it can produce not necessarily its volume. In a much larger scale, the ultrasonic mass flowmeters can confidently control the flow rate of the oil passing through the pipe lines from oil-producing countries to the consumers. Without mass flowmeters most of oil producers might have serious arguements with their costumers about the real amount of the delivered energy. Using a reliable flowmeter can directly improve the oil pricing formulas. In this project, the measurement is carried out by means of a set up based on a two-transducer measurement system illustrated in figure 2.1. The two transducers simultaneously emit ultrasonic beams aiming towards the flow. According to the timing and amplitude of the echoes received, we manage to compute the flow rate of the liquid together with many valuable information on physical properties of the liquid. To achieve the above goals, in chapter 1, we raise some brief theoretical aspects of the transducer as an electronic component. In chapter 2, we look into the specific measurement system that we exploit in order to obtain our samples and the hardware used. It also shows how we can determine the density and the flow rate based on mechanical formulas. Chapter 3 discusses the signal processing techniques that are applied to extract the desirable information from the signals which are necessary for computing the flow speed and the sound speed according to the discovered formulas in previous chapter. Chapter 4 investigates the frequency response of the system. Chapter 5 shows how we practically estimate the density and flow rate of the liquid. Eventually chapter 6 summerizes the results, performances, properties and eventually presents the final conclusions. Moreover, it proposes some future works that can be done in this area..

(14) 2.

(15) Chapter 1. Physical and electrical properties of ultrasonic transducers 1.1. Overview. In this chapter, some physical properties of an ultrasonic system have briefly been introduced. Besides, the electro-mechanical behaviour of such a piezoelectric and its equivalent electrical circuit is justified. Since ultrasound wave is categorized under the class acoustic waves, one needs to have a sufficient background of physical aspects of the wave propagation and acoustic environment. Secondly, a practical investigation of the piezoelectric disc and its mechanicalelectrical interpretation must be considered. Unlikely to light beam, both sound and ultrasound are mechanical waves that need a medium to propagate through. In this chapter some essential aspects of such propagations are clarified.. 1.2. A glance at the acoustic world. An acoustic wave is an energy disturbance that travels through a medium. Mass is displaced and the medium tries to return to its normal or undisturbed state. As it does that, the disturbance moves on in the medium.. 1.2.1. Acoustic impedance and reflection. Once a propagating ultrasonic wave crosses an interface seperating two media, a proportion of the wave will be reflected, another part will pass through the interface towards the next medium and a percentage of the wave energy will be absorbed. Since many years ago, scientist have been measuring the acoustic impedance, determining acoustic reflection factord, the absorption factor, a, for different materials at different frequencies. The results have been published in tables and are available. The acoustic impedance of a material Za can be calculated as below:. 3.

(16) 4. Figure 1.1: Diffraction effect .. Za = ρ×c. (1.1). where ρ is the density and c is the speed of sound in the medium.. 1.2.2. Attenuation. As it happens with any other mechanical wave, when an acoustic wave propagates, energy is converted into heat due to the friction in the medium which is being compressed and expanded. This is also known as viscoelatic loss, which is the major loss contributor. Other, less significant, effects are losses due to relaxation and thermal conductivity. The pressure amplitude of a plane wave − propagating in the → r direction decays exponentially. Therefore one could state → − | Ar |=| A0 | e−α.| r |. (1.2). − Where A0 is the initial pressure amplitude in the origin → r = 0 , α is the attenuation constant which differs for different media and is a certain value and − | Ar | is the pressure amplitude in typical point of r=→ r.. 1.2.3. Diffraction. Another source of attenuation in an ultrasonic system is diffraction or beam spreading. This causes the energy in the ultrasonic beam to spread over an increasing frontal area as it propagates through the medium as shown in figure1.1.. 1.2.4. Frequency dependence. It will be shortly shown in the thesis that many of the physical attributes that we are dealing with in this project are totally frequency dependent. Sound speed for instance, varies a lot with frequency. Because it has been proved that the high frequency harmonics of a sound travel faster than the low-frequency ones. Besides, acoustic impedance is also relevantly frequency varying. Later on, we will investigate this phenomenon deeply in this thesis..

(17) 5. Figure 1.2: Side view of a single element transducer. The arrows show the displacement movement of the transducer and the cap.. 1.3. The piezoelectric transducer. 1.3.1. Mechanical-Electrical properties. A device which converts one type of energy to the other is a sensor. By definition, a transducer is a device that is not only capable of converting one type of energy or physical attribute to another but also can do it vice versa. In fact, a transducer can be regarded as a double-way sensor. A piezoelectric material has an asymmetric atomic lattice likewise figure 1.2. When it is subject to an electric field, the material changes its dimension. The process is reversible, such that if the material is strained or compressed, an electrical field is generated. Hence, there is coupling between mechanical and electrical properties. According to G. S. Kino in [7] this coupling is described by the piezoelectric constitutive relation. T = C E S − eE S. D = ε E + eS. (1.3) (1.4). Equation (1.4) couples the stress T in the material to mechanical strain S. Due to Johan P. Bentley in [5], stress is defined by force/area while the applied stress produces strained in the body which is officially defined by (change in Stress length) over (original unstressed length). Elastic modulus is Strain which gives certain value for certain materials including our piezoelectric transducer. Back to the equation (1.3), E stands for the applied Electrical field. Equation (1.4) couples the electrical displacement D in the material to the applied electrical field and the strain. The constant C E is the elastic stiffness at constant electrical field, εS is the free permittivity, and e is the piezoelectric coupling constant..

(18) 6. Figure 1.3: Capacitive microphone, the moving plate is displaced due to strength of the existing acoustic pressure field. 1.3.2. Ultrasonic transducers versus the microphones. Microphones are a type of sensors that do convert acoustic pressure fields to electrical energy. Due to this definition, microphones are similar to ultrasonic transducers since both deal with acoustic pressure fields. However, they are very different devices in the sense that an ultrasonic transducer, as it goes by the name transducer, is a double-way sensing device which not only does convert the acoustic pressure field to electrical energy but also is able to conversely convert electrical signals to acoustic pressure fields. Moreover, microphones are meant to work with musical (sound) acoustic fields in the audible range of 20HZ to 20KHZ whilst ultrasonic devices work with frequencies way above audible frequencies. Microphones are mainly classified into two different categories of dynamic microphones and capacitive microphones. The capacitive microphones are able to respond to very high audio frequencies, and they are usually much more sensitive than their dynamic counterparts. Therefore we discuss the capacitor microphones here, only. A capacitive microphone depicted in figure 1.3 is a variable capacitor that consists of a fixed plate in parallel with a diaphragm which accounts for its variable plate and a small air gap between these two plates. The diaphragm which is just some microns thick and often fabricated from mylar moves backwards and forwards relative to variation of the acoustic pressure in the environment. This change in the distance between the plates of the capacitor ∆d, causes the capacitance to vary due to equation (1.5). ∆c =. ε.A ∆d. (1.5). Where A is the area of the plates and ε0 is permittivity of the air. If a fixed charge V is placed across the capsule, the change in the capacitance leads to generation of a small electrical current, i, according to equation (1.6), equation(1.7)..

(19) 7. ∆Q = V ×∆c. (1.6). ∆Q (1.7) ∆t Where ∆Q is the charge variation, V is the fixed voltage applied on the capacitor. ∆t, on the other hand, is the time duration under which these changes have occurred. Therefore, one can state that the change in acoustic pressure fields leads to generation of electrical current and consequently electrical energy. Besides, here another subtle difference between such microphones and ultrasonic transducer becomes clear. As stated above getting the microphone to work we need to apply a fixed voltage, V across the plates. This so called fixed voltage is also interpreted as bias voltage in electronic world. However, for an ultrasonic transducer, there is no need for such a bias voltage. That’s because as explained in 1.3.1, in an ultrasonic sensing element the conversion happens by change of the shape of the asymmetric lattices. i=. 1.3.3. Sensor design. A common way to build a transducer is to use a disc made of an artificial ceramic such as lead zirconium titanate (PZT). The disc is equipped with thin metal electrodes on each surface to form a transducer. Once a voltage source is connected to the electrodes an E field is generated, and the disc either increases or decreases its thickness depending upon the polarity of the applied field. The transducer will also work the other way around, i.e. when subject to a mechanical deformation, a voltage will develop across the electrodes of the transducer. This is used for reception or sensing the ultrasonic pressure fields. For both generation and reception of ultrasonic pulses, the vibration must be coupled between the transducer and the medium in which the ultrasound is supposed to propagate. If the medium has got acceptably similar acoustic impedance as of the transducer, a fairly good coupling can be accomplished by creating a good mechanical contact between the surfaces using a couplant such as glycerin. Otherwise, if the acoustic impedance of the transmitter and the medium, a matching layer is needed which must have an acoustic impedance in between the transmitter and the medium. A complete sketch of the ultrasonic transducer with matching layer and backing is depicted in figure 1.4 where the backing layer has been used to tailor the ringing of the transducer after excitation. However, exploiting a backing layer will result in fewer reflections passing towards the medium, it has this advantage of killing the ringing after the excitation as illustrated in figure 1.5.. 1.4. Modeling of the ultrasound system. It is crystal clear that before beginning any electronic design for our project, we should be able to interpret the electrical-mechanical properties of the piezoelectric into a suitable model. The aim of coming up with such a model is to be able to predict the performance of a system or transducer before it is built. This will enable us to find out the right electronic circuit and finally import it into our powerful simulation software, SPICE. By means of SPICE, we will be able to know almost everything we are looking for and create a firm backbone for our.

(20) 8. Figure 1.4: A piezoelectric transducer made up in a plastic container and together with the backing layer. Figure 1.5: Echo signal produced by a transducer with poorly matched backing (a) and with well-matched backing (b). Note that the ringing swings have been canceled in (b) however the amplitude of the passed pressure field is lower than in (a). system design. Ultrasound system can be simulated using various approaches. Using mechanical simulation engines such as FEM or FID, mathematical software tools such as MATLAB and even some ultrasonic-dedicated software are all possible options on the table. However, for an electronic engineer working on electronics for an ultrasound system none of the above looks attractive, as they do not simulate electronics. Therefore, the focus on this project and similar ones has been mainly on SPICE. The analogy between mechanical structure and electrical circuits is widely used for the analysis of mechanical systems. Using the analogy, forces are replaced by voltage sources and velocities are replaced by electrical currents. Then, an equivalent circuit of the system is constructed. This method even becomes a more powerful tool for the analysis of electromechanical systems where some parts of the system are already in electrical domain. For example, equivalent circuit analysis was successfully employed for piezoelectric transducers for their design and optimization..

(21) 9. Figure 1.6: the Leach equivalent circuit of a piezoelectric device. The development of models for ultrasound system dates back to 1940 in Masons work. Since then, different teams have used different ways to look at the case ranging from pure mathematical ways to highly-electronic ones. For instance, Lorenzo Capineri in his work [13] looks at the transducer as a black box, estimates the zero-pole transform of it. The true reference responses of the impedance and transfer function are derived by taking Fourier transform of the measured signal. Then final step optimizes the model parameters by minimizing a mean-square-error cost function between the model functions and the reference ones. The above optimization algorithm is a very famous one but is purely mathematical. However, we might be more interested in Leach model which results in an Electrical-Mechanical model shown in figure 1.6 according to Jonny Johansson in [2]. For briefness, we do not explain the proof in here yet we will use the result. Standard SPICE does not allow frequency dependent sources. Thus, the frequency dependence has been modeled using a different approach in the implementation showed in the above figure. When the transducer operates in receiving (listening) mode the impedance of an oscilloscope probe or an amplifier is normally high. In transmitting mode however, the transducer is connected to a voltage source that supplies the excitation voltage. This voltage source normally has low output impedance.. 1.5. Conclusion. In this chapter we took a brief look on concepts such as acoustic impedance and reflection, attenuation and diffraction. Besides, the reader was introduced to ultrasonic transducer and the way of designing a suitable transducer. We also illuminated differences between a microphone and an ultrasonic transducer..

(22) 10.

(23) Chapter 2. The experimental set up 2.1. Overview. As discussed in previous chapter sound, unlike light, needs a medium to propagate through. We take advantage of this phenomenon when we use it to calculate both the velocity and the speed of sound in the flow using our ultrasonic meter. The questions this chapter addresses are firstly using what sort of set ups ultrasound can be generated and perceived. Secondly, how the measured ultrasonic echoes could be exploited to compute desirable properties of the flow. Therefore the first section of the chapter will introduce the measurement set up followed by explanation of the electronics of integrated circuits used in order to excite the piezoelectric element and also in order to listen to the ultrasonic wave fields and echoes. The third section of the chapter will be dedicated to showing how appropriate mechanical equations can be written so that the measured ultrasonic data become handy in order to calculate the velocity of the passing flow, sound speed in it. It also illuminate a way to compute the density and flow rate of the liquid.. 2.2. The physical scheme of the measurement system. The main system shown in figure 2.1 below illustrates a cross section of the device which was built up as the back bone of our ultrasonic system in this project. The system which is made up of an aluminum block consists of two probes in which the piezoelectric discs are located. Its inlet and outlet have a diameter of 39 mm and the reduced part has a diameter of 20 mm. The liquid (flow) passes through this inlet-outlet pipe. The reflectors are inserted at the bottom of the figure centers of which are 140 mm apart from each other. (dm = 140 mm) The two probes on top of which there are the transducers is made up of polyetheretherketone(PEEK) polymer rods with a diameter of 20.8 mm. The transducers are the two piezoelectric disks fabricated by a Danish provider, Feroperm with the official central resonance frequency of 2.19 Mhz. We used an appropriate backing layer to support the transducer inside the probe. Some 10 mm prior to the probe-liquid interface, as seen in figure 2.1, there is a step 11.

(24) 12. Figure 2.1: Cross section of the system. The left probe is complete whilst the right probe is a cross section view. The liquid flows through the inlet-outlet pipe along the path BC (or reversely CB) .. (which can also be regarded as a probe-air interface) that can reflect a small proportion of the beam back to the transducers. The probe is based on a design by Lynnworth used in a four ultrasonic sensors mass flow meter. Both transducers are fired (excited) simultaneously aiming downwards. Each transducer receives back two echoes, one from the step and the other one from the probe-liquid interface. In adition to those two echoes, each transducer receives a through signal which has been transmitted by the other transducer and has passed through the flow. Before describing how using the timings of these signals we derive the velocity and sound speed in the flow, we should firstly investigate the integrated circuits by which we carry out the signal acquisition of the system.. 2.3 2.3.1. The electronics Exciting electronic circuit. Exciting an ultrasonic transducer is known as a process that causes the piezoelectric disk to oscillate and generate the ultrasonic pressure field which known as ultrasonic pulse. As stated earlier in equations (1.3) and (1.4) a drastic voltage change on the piezoelectric element will cause its asymmetric material to change its dimension which sequentially builds up an ultrasonic mechanical energy field..

(25) 13. Figure 2.2: The outer look of the system.. Therefore, to be able to excite the transducer by applying these voltage changes we need a driving circuit which is supposed to provide the driving signal to the disk. The shape of the resulting ultrasonic pulse is firstly influenced by the shape of this driving signal and secondly the matching layer. A short transmitted ultrasonic pulse with high bandwidth is often desired in pulse-echo system as of ours. Short driving pulse signals are commonly generated using the fundamental oscillating frequency of the transducer. One often used method is to apply a transient electric pulse to the transducer.This driving pulse is generated by the discharge of a capacitor across the connectors of the piezoelectric disc. The voltage and the size of the capacitor set the amount of energy transferred into the transducer, deciding the pressure amplitude of the resulting ultrasonic pulse. Because the driving signal in above case would be the result of a voltage discharge over a capacitor, it will have different unpredictable harmonics. To be able to better control the frequency response of the system one may use the pulses with controlled shapes such as square waves or parts of sinusoidal waveform. A square wave excitation can be achieved with a set of power-stage transistors which connect the transducer to the positive and negetive supply voltage as shown in figure below. The logic control unit changes the amplitude of the gate-source bias which makes one of the transistors switch on and the other one swith off. VDD is the voltage source which determines the amplitude of the exciting pulses..

(26) 14. Figure 2.3: A popular push-pull system for generating the square driving pulses [2] .. 2.3.2. Listening electronic circuit. As previously mentioned in (1.3.2) for converting a received ultrasonic beam to electronic signal we don’t need any biasing circuit. The conversion is carried out by the movement of the asymetric piezoelectric elements. However, the resulting electronic signal is sometimes too weak and should be amplified carefully. However in this project we do not need to amplify the received signals.. 2.3.3. Combined transmission and reception. The transmission and reception can be done at one transducer over a single wire. The only thing to be noted is that the reception period and the excitation phase should not interfere. Otherwise, a part of the data will be corrupted. For instance, while the transducer is excited any recieved echo will be interfering with the excitation signal in our samples. Therefore, the transmissions should be scheduled so that it takes olace at appropriate times and does not collide with any of vital echoes.. 2.3.4. ASIC: application specific integrated circuit. The driving circuit of the transducers is based on a design by Dr. Jonny Johansson illustrated in figure 2.4. The design has been performed in a 0.8µm, 50 V CMOS technology provided by austriamicrosystem (AMS).[12]. The produced AISIC is a thumb-size one with a very low power consumption. The design is based on four main blocks (Control block, Boost converter, Discharge unit and.

(27) 15. Figure 2.4: Block shcematic of the complete chip outlining the main blocks and external connections. Control signal inputs can be bonded directly on the chip. [2]. Amplifier ) as depicted in figure 2.4. The control block is responsible to control the functunality of the system according to the paarameters set by the user via the digital input pads. The boost converter generates the high voltage on the piezoelectric transducer prior to excitation. Discharge unit is the most important unit which carries out the excitation of the piezoelectric disks by a rapid discharge. The amplifier is not used in the current project. The operation of the chip is examplified by the behaviour of the voltage on the transducer as shown in figure 2.5. The main phases are: 1. Charge: The transducer is charged to high voltage with the boost converter. The charging is done with several pump cycles. 2. Hold : The piezoelectric is held at the high voltage level 3. Discharge: The transducer is rapidly discharged to generate the ultrasonic pulse. Basically, the excitation of the piezoelectric is a consequence of this discharge. 4. Wait: is set zero in our project. Amplify: The amplifier powers up the received echoes. In this project we do not use the built-in amplifier. Therefore, we set the amplification.

(28) 16. Figure 2.5: Chip functionality illustrated by the behaviour of the voltage on the transducer.. factor one. 5. Off : System is turned off untill next pulse generation. Using the inductive boost converter, the inductor in figure 3.3 is charged. The charge is transferred on the peizoelectric disk When the charge reaches the highest level, the pumping phase is stopped. As seen in figure 3.4 the piezoelectric is discharged suddenly. According to (1.3), this drastic voltage change will cause the piezoelectric to generate ultrasonic pulse. Then the transducer will be in listening phase during which the echoes are received. The whole process starts after 10 milli second. In other words, the flowmeter updates its data every 10 ms. Therefore, we can think of our flowmeter as a real time system.. 2.4. A/D conversion and the list of the measured signals in this project. The Analogue to digital converter (ADC) of our set up is a regular RTI7054 oscilloscope A/D device with a sample rate of 156 Mhz. The samples are saved in the cash memory of the device and are later imported to MATLAB with a .mat extension. To refer to the pulses correctly we used the following notation exy is the xth pulse received by the y th transducer. For instance e3a is the third signal (in the figure) recieved by transducer(a). According to figure 2.6 this would be the signal number 3 in the plot of fig 2.4 which is the reflection of probe-liquid reflection. The following useful pieces of data are extracted from it and saved. They are the only data sources we have in order to obtain knowledge about anything related to this ultrasonic system. 1. e1a and e1b , the excitation pulses..

(29) 17. Figure 2.6: The useful signals in our samples numbered according to the text .. 2. e2a and e2b , the echoes which are reflected by the edge of the step (probeair interface) in each of the probes. 3. e3a ,e3b The echoes which are reflection of each of the probe-liquid interfaces. 4. e4a and e4b , the through signals that have been transmited by each of the transducers and received by the other. (e4b is a proportion of the transmitted ultrasonic beam by transducer (a) makes it to travel along the path ABCD to transducer (b) and vice versa a proportion of the emitted ultrasonic beam by transducer (b) which travels along the path DCBA and makes it to reach transducer(a)), is called e4a ).. 2.5. Calculation of the flow velocity and the sound speed in the flow. In this section we will present how the velocity of the flowing liquid and the speed of sound in the liquid can be computed by only having the mentioned echoes and through signals..

(30) 18 The only essential pieces of data that are necessary to have for these calculations are summarized down here as a remider: 1. Echoes of the probe-liquid interfaces (e3a and e3b ) 2. Through signals (e4a and e4b) These two bunches of signals match number 3 and 4 in figure (2.3). In the following subsection we see how we can derive the speed of sound and the velocity of the flow by only knowing the timing of the above pieces of data.. 2.5.1. Calculation of the flow velocity. Figure 2.7 depicts the overall scheme of our ultrasonic system. The purpose of → − our calculation is to find the velocity of the flow shown V in figure (2.5). In all calculatiuons in this project we assume that the direction of the flow is along −−→ the BC. It’s obvious that if the other direction is the case then all our results will be simply multiplied by a minus sign. Apparently, when the ultrasound beam being emitted by transducer (a) travels from point B to point C (downstream direction) its speed which equals the speed of sound added to the speed of the flow. Conversely, when the other ultrasound beam emitted by transducer(b) is traveling from point C to point B (upstream direction)it is opposed by the stream of the flow so its speed is the speed of the flow subtracted from the spead of sound. The main idea is that subtracting the downstream velocity from the upstream velocity gives us the velocity of the flow due to the following equation.  −−−−−−−→ → − → − Vdonstream = C + V (2.1) −−−−−−→ → − → − Vupstream = C − V By subtracting the lines in the above set of equations from each other we → − can easily derive V the velocity of the flow. dm /2 dm /2 → − −−−−−−−−→ −−−−−−→ → − 2 × V = Vdownstream − Vupstream =⇒ V = − (2.2) tBC tCB Where dm is the distance between the reflectors (beween ponits B and C)and tBC is the time it takes for the beam to travel the distance from point B to C. Vice versa tCB is the time it takes for the ultrasonic beam to travel from C to B. It’s very vivid that the tBC can be seen as the result of subtracting tAB and tCD from the time corresponding to the whole distance from A to D or tABCD . The same analysis would be correct for the reverse direction for tCB the time it takes for the beam to go from point C to B. In other words we can state  tBC = tABCD − tAB − tCD (2.3) tCB = tDCBA − tBA − tDC According to the figure 2.7 the following geometrical equation will be true where dm is the distance of the reflectors between points B and C. dl is the whole distance that the beam travels in the liquid, ABCD.  AB + CD = dl − dm (2.4) m tAB + tCD = dl −d C.

(31) 19. Figure 2.7: Cross section of the system. The left probe is complete whilst the right probe is a cross section view. The liquid flows through the inlet-outlet pipe along the path BC (or reversely CB) It’s obvious that if the other direction is the case then all our results will be simply multiplied by a minus sign.. To refer to the pulses correctly we use a similar notation to the one we introduced in section 2.4 for the arrival times: txy is the time of the xth pulse is received by the y th transducer. Having those time of flights one could state:  t3b tABCD = t4b − t3a 2 − 2 (2.5) t3b t3a tDCBA = t4a − 2 − 2 Finally, inserting the values from equations (2.4) and (2.5) into the main equation (2.3) will lead to obtaining the result below.  m 3b tBC = t4b − t3a +t − dl −d 2 C (2.6) t3b +t3a dl −dm tCB = t4a − 2 − C If we replace tBC and tCB from the above equation (2.6) in equation (2.2) the velocity of the liquid will be precisely estimated. → − V =. 2.5.2. t4b −. dm 2 t3a +t3b 2. −. dl −dm C. −. t4a −. dm 2 t3b +t3a 2. −. dl −dm C. (2.7). Calculation of the sound speed. Equation 2.7 states that in order to calculate the flow velocity we need to have → − the sound speed C in first place. Back to figure 2.7 one way to compute the.

(32) 20 sound speed is as seen in equation below. 2 ∗ dl ABCD + DCBA → − → − =⇒ C = C = tABCD + tDCBA tABCD + tDCBA. (2.8). Very identical equations to the ones in previous section can be re-written in order to eventually end up in the following result. → − C =. 2.6. t4b −. 2 ∗ dl + t4a −. t3a +t3b 2. t3b +t3a 2. (2.9). Calculation of the acoustic impedance, density and mass flow rate. → − As seen in last section, in order to calculate the velocity of the flow, V and → − sound speed in that liquid, C we only look at the timing of the pulses and we do not care about their amplitudes at all. However, in this section when it comes to calculating the acoustic impedance of the liquid, Zl the amplitudes of other pulses also come to play. In this section we describe how we can derive the acoustic impedance of the flow (liquid) independantly for each of the probes having the following peices of data: Using the data from probe(a): 1. The excitation pulses of transducer(a), e1a . 2. The echo which has been reflected by the edge of the step (probe-air interface) in probe(a), e2a . 3. The echo which has been reflected by the probe-liquid interface of probe(a), e3a . Using the data from Probe(b): 1. The excitation pulse of transducer(b), e1b . 2. The echo which has been reflected by the edge of the step (probe-air interface) in probe(b), e2b . 3. The echo which has been reflected by the probe-liquid interface of probe(b), e3b . (Please refer to figure 2.6 to see which pulses are meant above) Another subtle yet important difference is that calculation of the acoustic impedance can be performed by having the above data from one probe indepen→ − → − dantly. However to calculate V and C we need to consider the through signals too and somehow combine the received data of both transducers. According to equation (1.1), having the sound speed C and the acoustic impedance leads us.

(33) 21. Figure 2.8: The reflection of the interface between environment1 and 2.. to knowing the density of the liquid. Knowing the velocity of the flow, density of it and the dimensions of our measurement device paves the way towards estimating the mass flow, m ˙ which was the final target of this project. The calculations regarding the mass flow metering are introduced in the last sub section of this chapter.. 2.6.1. Reflection coefficients and Calculation of Acoustic impedance. 2.6.1.1. Pressure reflection coefficient VS. Intensity reflection coefficient. According to Lawrence Kinsler[3], reflection coeffiecient is defined as the ratio of the pressure of the reflected wave to the pressure of incident wave and is shown by R. However, there is also the intensity reflection coefficient which is introduced as the ratio of intensity instead of the ratio of pressure, this latter definition is represented by Rπ . For the interface shown in figure 2.8 the reflection factor would be: ( amp(X2 ) 2 R = pressure( X X1 ) =⇒ R = amp(X1 ) (2.10) std(X 2) 2 2 Rπ = intensity( X X1 ) =⇒ Rπ = ( std(X1 ) ) To justify the above equation we should start by knowing thatAcoustic Intensity of a sound is officially the average rate of flow of energy through a unit area normal to the direction of propagation. If we have discrete electronic samples of the ultrasonic pressure field as we really do (let’s call it x[n]), then according to Ziomek [4], the intensity of it will be: Intensity(x[n]) =. 1 N. X. x2 [i]. (2.11). i=1,toN. On the other hand, std which stands for standard deviation of the signals tells us how much the amplitude of the elements in that sequence are dispersed.

(34) 22 from each other and is defined as below. s std(x[n]) = δx =. 1 X 2 (x[i] − x) N i=1. (2.12). In above equation, x refers to the mean value of the signal. As in our case, once the mean value of the signal is zero, then we can claim that the intensity of the ultrasonic samples square of the standard deviation. P equals to the 2 Intensity(x[n]) = δx2 = N1 i=1 (x[i] − x) . (the square of the std, δx2 is also known as variance in many statistical texts). Therefore the equation (2.10) looks correct. Back to equation (1.2), one could see that the pressure reflection for figure 2.8 can be written as follows. R12 =. Za2 − Za1 Za2 + Za1. (2.13). In above equation R12 refers to the pressure reflection coefficient for the case the beam is transmitted inside the medium 1 propagating normally towards the interface with medium 2. (As illustrated in figure 2.8). According to equation (2.10) the latter equation of (2.13) can be re-written as ( Za2 −Za1 R12 = amp(x2) amp(x1) = Za2 +Za1 (2.14) R21 = −R12 2.6.1.2. Pressure transmission coefficient and intensity transmission coefficient. The ratio of the sound pressure passing through the interface to the sound pressure of the incident beam will be called the pressure transmission coefficient, T12 . Similarly, the intensity transmission coefficient will be the ratio of the intensity, TI12 . If we neglect the scattering phenomonon at the interface this intensity transmission coeffiecient will be also seen as the power transmission intensity, Tπ12 . If the absorption of the interface can be ignored then the sum of the reflected and transmited energy should be equal to the enegry of the incident beam. Besides the law of Continuity of pressure is also true for the interface. These two fundamental laws end up in following set of formulas:  Rπ + Tπ = 1 (2.15) 1+R=T Now that we have discussed the reflection and transmission factors we can get to calculate the acoustic impedance. Figure 2.9 shows the probe, the steps inside it and of course the probe-liquid interface. When the transducer(a) fires e1a it flies the distance d untill it reaches the step, therefore X1 which is the incident signal on the step is x1 = e1a .eαp .d As a result x2 = R ∗ x1 =⇒ Re1a .eαp .d . When the echoe gets back to the transducer its amplitude will be estimated as below.  e2a = x2 .eαp .d =⇒ e2a = R12 .e1a eαp .2d (2.16) =⇒ amp( ee12 aa ) = R12 .ξ.e−2αp d where ξ has been multiplied to our previous formulas in order to represent the percentage of beam the step is exposed to, in compararison to the proportion.

(35) 23. Figure 2.9: One of the probes, the steps inside it and the probe-liquid interface. the probe-liquid is. This is mainly due to the physical attributes of the probe( figure 2.9). Obviously the rest of the beam will hit the probe-liquid interface. If we replace the reflection coefficients according to Equation (2.14) for the both the step reflection and the probe-liquid reflection we will have the following set of equations. ( Z −Z amp( ee2a ) = ξ × Zaa +Zpp .e−2αpd 1a (2.17) Z −Z amp( ee3a ) = (1 − ξ) × Zll +Zpp .e−2αpdp 1a In the latter equation above, Za , the acoustic impedance of air, zp , αp , (besides e1a , e2a , e3a are all extracted) are all known. This will lead to knowing the correction ratio ξ. Replacing it in the second equation will enable us to derive zl the acoustic impedance of the liquid.. 2.6.2. Calculation of the density. Having the acoustic impedance of the liquid Zl from the previous subsection and c, the speed of sound in the liquid from subsection 2.5.2, replacing it in the very begining equation of 1.1 will result in knowing the density of the liquid.. 2.6.3. Calculation of the mass flow rate. → − −→ → − The mass flow rate is simply defined by m ˙ = ρ( A .Vm ) where A is the cross −→ section area of the pipe which the liquid flows through. Vm is the volumetric velocity which is the velocity of the flow (in 2.5.1) per meter square. The ’.’ is the inner product operand. Knowing the flow velocity from 2.5.1 and the density from 2.6.2 and measuring the A, cross sectional area of our measurement device, we can easily compute the mass flow rate..

(36) 24.

(37) Chapter 3. Signal processing solutions to estimate the flight times 3.1. Overview. In this chapter, we discuss how we handle the acquaired signal to obtain the timing of the pulses (flight times). For estimation of flight times, we evaluate three techniques. We will see that cross-correlation is a very straightforward solution. Optimization solution will be introduced and later an exact phase unwrapping solution will be brought up in details. This phase unwrapping analysis paves the way towards a very useful frequency analysis of the system. This frequency behaviour of our ultrasonic system is explained in the next chapter. The aim of the signal processing unit is to: Step 1: Filter measurement and extrac the favourite pieces of ultrasonic signal. These gated pieces of information account for the echoes and the through signal pulses and almost all information we have about the system. Step 2: Perform the actual signal processing methods to find out timings needed in calculation of the flow speed and the sound speed. Capture the execution time for each of these methods to see which one is faster to run on a CPU. Step 3: Use the flight times found in previous steps and the formulas mentioned before in section we compute the flow speed and the speed of sound in the flow. The above steps can be illustrated in the following block diagram.. 3.2. Block diagram of the signal processing unit. 3.2.1. Filtering. 3.2.1.1. The design of the appropriate band pass filter. The first block in diagram 1 is supposed to purify the samples from the unwanted terms. The crude pre-processed measurements of both channels look similar to fig 3.2 where the first big spikes in the beginning part account for 25.

(38) 26. Figure 3.1: The diagram of the signal processing unit. Figure 3.2: Pre-processed samples obtained by sampling the measured signals at each of the channels (transducers). the excitation. Remember that as stated previously during excitation a large voltage is discharged on the transducer. This will look like the spikes seen in figure 3.2 besides one could see that our measurements have an offset of 1 volt. Obviously these spikes are all of low frequency and can be easily deleted. On the other hand, one could see in figure 3.1 that high-frequency noises have contaminated the measurements in various parts of it. By now, our choice is a band-pass filter which can cancel both the high frequency noise terms and the low frequency spikes at the same time. We already know that our transducers are excited by a broad band discharge pulse therefore the strength of the spectra of our samples will have strong low frequency terms as seen in figure 3.3. Note that we are interested in extracting the ultrasonic pulses which we already know that occur around the frequency 2MHZ (the resonance frequency of the transducers is 2.019 Mhz written on them by the manufacturer). We need to discard the frequencies either lower or higher.

(39) 27. Figure 3.3: The Amplitude of the Fourier transform (a representetive of Power Spectrum Density (PSD)) of each of the received signals by transducer(a) and by transducer(b). This shows how strong are the low frequency terms. our favourite ultrasonic data exits around 2Mhz. than it. This is typically done by a bandpass filter in electronic world. To specify precisely what the appropriate band should be for our filter, we need to zoom at the power spectral density or equivalently the amplitude of the Fourier transform of the samples around the centre frequency of 2.19Mhz to find the frequency band at which the spectra of the samples are relatively stronger. This band would correspond to the exact frequency range at which our ultrasonic pulses exist. We then design our filter so that it passes the frequencies within this band and stops (attenuates) the other frequencies outside of that band. Note that the filter above is a common Butterworth filter whose coefficients are estimated by butter instruction in MATLAB. For design of bandpass filters, there is instruction fir in MATLAB which generates even more ideal frequency response. However the subtle reason behind choosing the butterworth filter is because electronic implimentation of a butterworth is easily done by capacitors and resistors using available tables. This delicate point is important once we decide to build an electronic embedded system in future. The output of the filter is depicted in figure 3.5 where one could see that all the unwanted terms in the measured samples (already shown in figure(3.2)) have been removed successfully. Ultrasonic pulses and echoes are vivdly seen in the figure, the SNR (signal to noise ratio) is pleasantly high. The next step would be gating these pieces of data..

(40) 28. Figure 3.4: The amplitude of the Fourier transform of each of the measurements being zoomed around 2Mhz versus our designed 6th order FIR filter. Our band pass filter is with central frequency fc = 2.062M hz and -3db cut off frequencies of fl = 1.9375M hz and fh = 2.1875M hz.). 3.2.1.2. Extraction of the useful data. This latter figure 3.5 is basically the same as former figure 2.6 where the echoes are labeled. Here, our task is to gate those ultrasonic pulses properly in order to be used in later. For the specific fs , sample frequency of our A/D unit we can forsee approximately at what sample these favourite ultrasonic pulses show up in our measurements. What we can easily do is to chop our samples and extract all the favourite pulses into suitable frames. Accept for the excitation pulse which comes at the first sample, we can claim it is logical to assume that the step reflections (e2a and e2b ) arrive sample offset-edge=5500, probe-liquid echoes (e3a and e3b ) approximately arive around the sample offset1 =7000 and the through signals (e4a and e4b ) is received around the sample offset2 =25000. Then we define a value of win-size=3000 samples as the window size for our frames. The result would have been depicted in series of figure bellow.. 3.2.2. The methods. In order to find the time delay between two similar signals in sequences x[n] and y[n], the most straight forward method that one could come up with immediately is the cross correlation. Detecting the maximum of the correlation result which corresponds to maximum similarity between two sequences will give us the time.

(41) 29. Figure 3.5: The output of the bandpass filter seen in figure 3.4. The unwanted terms are cancelled and the ultrasonic pulses look clearly outstanding. Compare this figure with figure 2.6 to see how we name each of the pulses.). delay we were looking for. Besides, cross correlation can be also interpreted in frequency domain. Another method is optimization. Its backbone idea is to create yb (an ideal model for y[n]) based on parameters from x[n] and, then using the Least Mean Square algorithm minimize the error between this modeled yb and the real y[n]. Once the error has been minimized all the parameters including the time delay between y[n] and x[n] are known. The third method which we call it Phase unwrapping solution - unlikely to the two methods mentioned previously- solves the problem for a more general case where the system is assumed to be a frequency dependent system consequently x[n] and y[n] can not be modelled as easily as in the previous optimization solution. The idea is to take Fourier transform of both sequences x[n] and y[n] and then investigating the amplitude and phase of the resulting spectra, we will end up in estimating the time delay as well as exploring some very valuable information which will be used in the next chapter.. 3.2.3. Compute and display. Once the timing of the ultrasonic pulses has been discovered, we will simply → − replace them into the equations in (2.5) to compute the sound speed C , the → − velocity of the flow V . The values are then displayed on MATLAB command window. The flight times estimated by those methods will also be used in the future chapters to estimate the acoustic impedances, the density of the liquid.

(42) 30. Figure 3.6: The extracted (gated) ultrasonic pulses , e1a , e1b , e2a , e2b , e3a , e3b , e4a , e4b .. and finally the mass flow rate.. 3.3 3.3.1. Method 1: Cross correlation solution Time-domain cross correlation. This method is performed in time domain and takes its root from the fact that if we have two sequences x and y, cross-correlation of their absolute values xcorr(x,y), has the maximum value for a certain shift at which the similar forms in x and y overlap each other. Briefly, according to [8] convolution of two sequences x and y is defined as follows: X x[n] ∗ y[n] = Con[n] = x[l]×y[n − l] (3.1) l. Thinking about the formula above a bit thoroughly makes it clear that convolution of two sequences x and y each one having N elements results in a sequence.

(43) 31 of size 2N − 1 whose elements are calculated this way:  Con[1] = x[1]×y[1]     Con[2] = x[1]×y[2] + x[2]×y[1]     Con[3] = x[1]×y[3] + x[2]×y[2] + x[3]×y[1]  ...   Con[n] = x[1]×y[n] + x[2]×y[n − 1] + . . . + x[n]×y[1]     ...    Con[2N ] = x[1]×y[1]. (3.2). By definition cross correlation of two sequences has the following relation with convolution of them Xcorr(x[n],y[n])=Conv(x[-n],y[n]) Therefore, cross correlation is also a shift / multiply / accumulate arithmetic operation. The function takes two sequences and the sample frequency as inputs, detects the similar forms which occur in those two sequences, and returns the number of delay samples by which the common form is apart in the input sequences. Having the sample frequency we can figure out the actual time corresponding to this number of samples. The time in the original measurement which corresponding to a number of samples in the quantized information is calculated in bellow. N (3.3) T = fs Which gives us the equivalent time in seconds corresponding to N number of samples. In our measurement set up, sample frequency, fs is 156,248Mhz. For instance if the through signal by transducer 1 and the one received by transducer 2 be sent to the function, the result of the cross correlation of their absolute values would look like figure 3.7. Note at what sample the maximum occurs.. 3.3.2. Frequency domain Cross correlation. As stated earlier, the time domain cross correlation method is based on a shift / multiply / accumulate algorithm shown previously. Multiplication and accumulation are relatively fast operations on a CPU, though shifting is a very time consuming one. Because of this fact, we look for some other ways to skip shifting operation in the cross-correlation method. A very famous math theorem comes handy in this respect [8] . ’Convolution of two vectors in time domain is equivalent to multiplication of their spectra in frequency domain’. We wrote a function which takes the two sequemces x and y and the sampling rate as inputs. Then it, takes the Fourier transform (FFT) of suitable length from both input vectors, multiplies of both spectra together in a sample-wise way. Taking the inverse Fourier transform (IFFT) from the result, as we supposed before hand, gives us a very identical output as in previous method depicted in figure 3.8. Finally, applying the same idea as before, we figure out the corresponding sample delays and the time.. 3.3.3. Frequency domain versus Time domain. Obviously, this method carries out an FFT/ multiply/ IFFT operation. Since the FFT and IFFT algorithms are acceptably fast, the current method is at least.

(44) 32. Figure 3.7: The result of time domain cross correlation of e4a with e4b . The maximum represents the number of shifts by which the sequences have the match with each other. Note that the length of the resulting correlation of x and y is (length of (x))+(length of (y)) - 1. 100 times faster in comparison with the previous time domain solution based on a Shift/ multiply/ Accumulate operation. That’s mainly because the Shift operation is relatively a very time taking arithmetic operation on any CPU. Remember that all the micro computers ranging from the old Z80 microprocessors to the most recent ADSP systems or even the desk computers obey very identical TTL logics. Therefore, the above comparison between the two solutions will be true regardless of what electronic device we choose.. 3.3.4. Sub-sample precision. To find the fraction of sample where the maximum occurs in figure 3.7 or 3.8 a second order polynomial can be used. After determining the maximum, we extract two of its neighbouring samples from each side. Then by using the function polyfit we find the second-order coefficients a, b. c in y = ax2 + bx + c . The maximum will be happening in M ax =. −b 2a. . This maximum will convey fractions of the sample as well..

(45) 33. Figure 3.8: The result of frequency domain cross correlation of e4a with e4b . The maximum represents the number of shifts by which the sequences have the match with each other. Note that it is identical to previous figure, the only difference is it is performed 126 times faster!. 3.4 3.4.1. Method 2: Optimization solution The ideal model of the ultrasonic signals. At our measurement set up shown in figure 2.7, if E[n] corresponds to the discrete signal which excites the transducers then the signal received at any of the points A, B, C ,D and also the through signals recieved by any of the transucers (e1a ,e1b ,e3a ,e3b ,e4a ,e4b which are technically all of the gated signals shown in fig 3.5 accept for the step reflections) will be a delayed and attenuated version of each other. It’s vivid that once the transducer(a) has transmitted, the signal observed at point D is attenuated version of signal observed at point B with a delay. In other words, for a typical x, the received y will look like: yb[n] = A×x[n − τ ]. (3.4). where A is the attenuation factor, τ is the number of sample delays due to the flight time (both τ and A are scalars). Apparently in such ideal model, A and τ will vary relatively to the distance of the flight. The longer the beam flies the bigger τ grows and the smaller A becomes. Ideally, A and τ are both constants and do not vary with frequency. If we want to see the system above in frequency domain instead, we should take Fourier transform of the signals and the impulse response. The Fourier.

(46) 34. Figure 3.9: The impulse response of the ideal model defined in equation (3.4). Figure 3.10: the equivalent of figure 3.8 in frequency domain, H(ω) is called the frequency response and is the Fourier transform of the impulse response, h[n]. transform of the impulse response H(ω) is called the frequency response of the system. H(ω) = f f t(h[n]) . This ideal model, yb (defined in equation 3.4) is compared to the measured y[n] in order to estimate the optimal τ and A which fit the best.. 3.4.2. Generation of the ideal model in MATLAB. We take one of our measured ultrasonic signals as input x[n] and by defining constant values for A and τ in equation (3.4), we shift and attenuate the input so that we get the signal y[n]. Since shifting/multiplying operation is a time consuming logical operation, we play with their spectra instead. In the code bellow we first divide τ by the sample frequency to have the shift value in seconds instead of samples, second we take FFT of the input and save it in X, then create the model H(ω) of figure 3.10 which is as in equation 3.4 bellow. yb[n] = A×e−jω×τ (3.5). Inverse Fourier transform of Y [n] gives us the y[n] in figure 3.11 which is the ideally-modeled received signal. The Matlab code: tau=-344.86; shifted samples A=0.75; Attenuation factor tau =tau*Ts; shifted time in second H = A*exp(-j*wcs*tau); The shifting model Y = H.*X; y = real(ifft(Y)); the ideal model for y[n] plot(t,y,’r–’,t,x,’k-’);Title(’theinputsignals’);xlabel(’sample’);ylabel(’Amplitude’);legend(’y[n]’,’x[n]’).

(47) 35. Figure 3.11: x[n] versus y[n] where y is shifted 344.86 samples and attenuated 0.75 times. In other words, y[n] = 0.75 x[n − 344.86]. If we take the first echo of probe-liquid of transducer a as the input x[n], and have tau=344.86 samples and A=0.75 then y[n] will look like the following figure 3.10.. 3.4.3. Least mean square algorithm. The backbone idea here is to create the ideal model yb introduced by equation(3.4) (or equvalently 3.5 in frequency domain) and compare it with the y[n] to find the best model parameters (τ and A) that fits the best in the model and make it as realistic as possible. By using multidimensional unconstrained nonlinear minimization (Nadler-Mead) we try to minimize the difference between the modeled signal, yb[n] and the true signal y[n].  yb[n] = A×X[n − τ ] P (3.6) J = n |(y[n] − yb[n])2 | The algorithm would be as follows: 1. Make an initial guess for the flight time (τ ) and the attenuation factor, A. 2. Create the cost function, J as defined in equation (3.6) which is a sum of square of absolute value of difference between the guessed yb[n] and true y[n]. 3. Update our previous guessed parameters using Nadler-Mead minimization (applied by fminsearch insruction in MATLAB) so that the cost function J becomes smaller. 4. Keep on running step 3 (updating guessed values) a certain number of times to minimize the cost function J as much as possible. If the initial guesses are not too unrealistic, the algorithm above will finally.

(48) 36 converge to an optimal τ and A by which yb[n] and y[n] have the most similarity to each other. The Nadler-Mead minimization in our code is set to 100000 epochs and is performed by exploiting a MATLAB function, fminsearch which returns the optimal τ which minimizes the cost function. The initial guess for τ is assumed as the answer from the cross correlation method and initial A is set to one.. 3.4.4. Performance, Advantages and drawbacks of the optimization solution. If the initial guess is chosen reasonably, not too far from the real results, the algorithm always converges to the very actual answer with a very high subsample precision. The simulated signal y=0.75 x[n − 344.86] introduced in 3.4.2 and depicted in fig 3.11 is fed into our system to check the functionality of our code. The initial guesses are τ = 310 and A = 1 The results are shown in figure 3.12, where one could easily see that finally τ converges to 344.8 and A to 0.75. This method is able to determine the flight time with a high sub-sample precision. The only drawback is that it gives a constant value for τ and A. However we know that in dealing with real signals in this project τ and A are constant. This means that the optimization solution at its best case is a more exact solution than the cross correlation which also gives a constant value for flight time .. 3.5 3.5.1. Method 3: Phase Unwrapping solution Frequency response of non-ideal ultrasonic system. Referring to the physical sketch of our measurement system shown in figure 2.7, It’s crystal clear that once the transducer(a) has transmitted, the signal observed at point D is attenuated version of signal observed at point B with a delay. We modelled the delay τ and the attenuation factor A in equation (3.4). We then came up with optimization technique. However, If we look at the spectrum of both signals X and Y , it’s proved that our ultrasonic system does not attenuate all the frequency terms equally. Besides, the delay also differs for different frequencies. For our general case we can draw the following block diagram shown bellow where h[n] is the impulse response of the system. (note that the block digram is more general version of the one shown in figure 3.9). For the system in figure 3.13 One could write the famous equation y[n] = h[n] ∗ x[n] where the ’*’ reperesnts the convolution operation. Due to a basic mathematical theorem the above equation could be re-written as in Equation (3.7). The term |H(ω)| is called Amplitude of the frequency response and the ϕ(ω) is known as the phase of the frequency response. Soon we will see that the phase conveys very important data such as the timing of the signals and the delay. Investigating the phase of the frequency responses is the fundamental idea of the phase unwrapping solution..

(49) 37. Figure 3.12: The initial guesses are τ = 310 and A = 1. The optimization algorithm converges to tau=344.86 and A=0.75. The real y[n] versus the estimated one by optimization solution. They perfectly match!. Figure 3.13: The blok diagram of our non-ideal system. h[n] is the impulse response whilst H(jω) is called the frequency response of the system. . 3.5.2. Y (ω) = X(ω)×H(ω) H(ω) = |H(ω)|ejϕ(ω). (3.7). Phase unwrapping. We claimed in previous sub section that the phase of the frequency response conveys the information about the delay time. Let us assume that we have x[n], the ultrasonic signal observed at point 1 and y[n] the signal after having passed.

References

Related documents

46 Konkreta exempel skulle kunna vara främjandeinsatser för affärsänglar/affärsängelnätverk, skapa arenor där aktörer från utbuds- och efterfrågesidan kan mötas eller

Both Brazil and Sweden have made bilateral cooperation in areas of technology and innovation a top priority. It has been formalized in a series of agreements and made explicit

För att uppskatta den totala effekten av reformerna måste dock hänsyn tas till såväl samt- liga priseffekter som sammansättningseffekter, till följd av ökad försäljningsandel

The increasing availability of data and attention to services has increased the understanding of the contribution of services to innovation and productivity in

Generella styrmedel kan ha varit mindre verksamma än man har trott De generella styrmedlen, till skillnad från de specifika styrmedlen, har kommit att användas i större

A similar problem occurs in function approximation using, for instance, polynomial basis functions.. In

However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to

Industrial Emissions Directive, supplemented by horizontal legislation (e.g., Framework Directives on Waste and Water, Emissions Trading System, etc) and guidance on operating