SUGISAKAMASANORI

Journal of Robotics, Networking and Artificial Life
Vol. 9(4); March (2023), pp. 343–348
ISSN (Online):2352-6386; ISSN (Print):2405-9021
https://alife-robotics.org/jrnal.html

Ryuugo Mochizuki1, Yuya Nishida2, Kazuo Ishii2
1

Center for Socio-Robotic Synthesis, Kyushu Institute of Technology, 2-4 Hibikino, Wakamatsuku, Kitakyushu, 808-0196, Fukuoka, Japan
Department of Human Intelligence Systems, Kyushu Institute of Technology, 2-4 Hibikino, Wakamatsuku, Kitakyushu,808-0196, Fukuoka, Japan

2

Received 27 November 2021
Accepted 01 February 2023

Acoustic impedance
Ultrasonic
Hardness estimation
Internal reflection
Reflection coefficient

In food industry, shortage of workers is one of a serious problem. Automation of food handlingis a critical nowadays. To alleviate the damage caused by food picking operation by robotic hand,
we propose non-contact acoustic impedance estimation with ultrasonic wave, which should be
preceded before the picking for optimization of grasp stiffness. We assume the correlation
between hardness and acoustic impedance. That is, the harder the medium is, the larger the
impedance should be. As cooked food, the unevenness of hardness should be considered. We
made samples with media of different acoustic impedances, and experimented with ultrasonic tofind the relation between reflection coefficient and hardness. As the result, reflection coefficientincreased by 0.03-0.15 at the boundary of two media as the bottom medium is switched fromurethane to Aluminum.
© 2022 The Author. Published by Sugisaka Masanori at ALife Robotics Corporation Ltd.
This is an open access article distributed under the CC BY-NC 4.0 license(http://creativecommons.org/licenses/by-nc/4.0/).

of prepared dishes may run short, thus, automation offood process and packing are examples of crucial tasks.
Here, we propose to introduce non-contact hardness
estimation with ultrasonic wave to optimize grasp
stiffness, which precedes to the picking operation torealize safe food picking by a robot hand.
As cooked food, Unevenness of hardness should beconsidered. Here is an example of radish (Fig.1), which
was cut into proper size (nearly 50[mm] x 50[mm] x20[mm]) and heated with a micro wave oven for threeminutes. We used rheometer to measure the strain andstress, then, found Young’s Modulus of randomly
selected locations on the radish. Young’s Modulus varied
depending on the locations as Table 1 shows.
Furthermore, We estimated reflection coefficient by

1. Introduction
In Japan, employment opportunities for both sexes havesteadily increased up to now. According to the report [1],housewife household has been decreased by nearly two
million, while, dual-income households have increasedby the same number. Single households have also
increased (for elderly, about three million, for nonelderly, about 1.5 million) for 20 years recently [2]. As
the market of Japanese prepared dishes, in 2018, theincrease rate has grown up to 127.3[%] since 2009 [3].
From the point of view, the demand for prepared dishes
can be assumed to have relation with the increase of
single or dual-income households. If the increase of
single or dual income households continues, the supply

Corresponding author’s E-mail mochizuki.ryuugo126@mail.kyutech.jp, y-nishida@brain.kyutech.ac.jp, ishii@brain.kyutech.ac.jp
URL: www.lsse.kyutech.ac.jp343

However, to the best of our knowledge, no report ofhardness estimation considering heterogeneity of
hardness has been found. The heterogeneity usuallyorients from uneven heating of food during cooking
process, thus, we should consider the case of soft interior
with hard surface.
3. Proposed Method

Fig. 1 radish sample

3.1. Outline

Table 1 Unevenness of Hardness(Young’s Modulus)

Acoustic impedance 𝜁(∈ {𝜁1 , 𝜁2 , … }) can be expressedby Eq. (1), with 𝜁 and 𝑐 , here, ρ[kg/m3 ] indicatesdensity of sample [8].

Location Young’s Modulus
Number[kPa]
1
1.8 x 1042
2.8 x 1043
1.4x1044
2.1X10458.3x104
overwrapping two samples with different acoustic
impedance, considering the value as hardness. Acoustic
impedance is the product of density and sonic velocity.
Larger reflection coefficient should be estimated at thelocation of neighboring media as the discrepancy of
acoustic impedances is larger, which means harder
medium reflects stronger acoustic wave.

𝜁 = ρ𝑐.

(1)

To analyze received wave, we consider arbitral samples
as multiple joints of media as Fig. 2 shows (𝑧𝑘 denotes
the surface of medium k). Acoustic impedance 𝜁𝑘 ofmedia k has the relationship with reflection coefficient 𝑟𝑘
as indicated by Eq. (2)[9], which means larger reflection
can be found as the difference of hardness is larger
between sample k and k-1, under condition media k isharder than k-1.
𝜁 −𝜁𝑘−1

𝑟𝑘−1,𝑘 = 𝑘𝜁𝑘 +𝜁𝑘−1

2. Related Work.

(2)

We treat all ultrasonic waves in voltage because soundpressure cannot be measured in [Pa] with receiver.
As the calibration (Environment is shown as Fig. 3), we
estimate attenuation coefficient 𝛼0 [Np/m] and voltage
𝑃0 [V], which corresponds to sound pressure of emitted
wave, and is the reference to find 𝑟0,1 of samples, thenobtain and treat emitted wave 𝑝𝑖,0 (𝑡) in similar way as
basis function for wavelet transform [10]. Received wave
𝑝𝑟 (𝑡) can be regarded as summation of reflected waves
𝑝𝑟,𝑘 (𝑡) at different surfaces of media. 𝑝𝑟,𝑘 (𝑡) isexpressed in two ways, that is modelled reflection 𝑝̂𝑟,𝑘 (𝑡),and reflection 𝜓𝑟,𝑘 (𝑡) based on a waveform 𝜓𝑖,0 (𝑡) ,which is the estimation of 𝑝𝑖,0 (𝑡) and obtained from thedirect wave under the environment (Fig. 3). The former
expresses relationship between reflection coefficient 𝑟𝑘
at the surface of medium k and attenuation coefficient
𝛼𝑘 [Np/m] of amplitude inside medium k. based on
amplitude. The latter can express waveform because
𝜓𝑖,0 (𝑡) is stretched or compressed in two direction (timeand amplitude) to fit 𝑝𝑟,𝑘 (𝑡) in this stage. The flowchart(Fig. 4) shows the procedure to analyze 𝑝𝑟 (𝑡). First, theMachado et. al. have analyzed ultrasonic reflection toevaluate stiffness [4] of three different wood species byfinding sonic velocity in indirect method (Here,transmitter and receiver are aligned parallel and attached
on a surface of a timber) and compared the result of direct
method (Transmitter and receiver are aligned face toface). The relative error rate was less than 10 % between
the two method.
Cho et al. used five kinds of cheeses to determine
mechanical properties (Young’s modulus, hardness,
Toughness) with ultrasonic wave. They improved
stability of thickness and ultrasound velocity
measurement As the result, depth of cheeses wasestimated in 99.98[%] of accuracy, then high correlation(coefficient of determination 𝑅2 > 0.9, 0 ≤ 𝑅2 ≤ 1) was
found between acoustic velocity and the mechanicalproperties of the cheeses [5],[6].
Satyam et al. estimated varying hardness of tomatoes
over a time period from fresh to full ripen cycle of thecrops [7].

344

amplitude of 𝑝𝑟 (𝑡) is thresholded with fixed value 𝑇𝑃 for0
the judgement. If the judgement is true, 𝜓𝑖,0 (𝑡) slides inthe time axis direction, then stretched or compressed in
the ratio of 𝐴𝑘 vertically and 𝐵𝑘 horizontally so that thecorrelation between 𝜓𝑟,𝑘 (𝑡) and 𝑝𝑟 (𝑡) is maximized.

Function Generator generates one cycle of sign wave
(Fig.5(a)) to the transmitter and direct wave is received
as Fig.5(b) shows (Here, the direct wave is cropped to
exclude reverberation.). 𝑃𝑜 and 𝛼𝑜 be calculated with the
[V]

[V]

t[us]
t[us]
(a)Input of Function Generator(b)Amplified Direct Wave
Fig. 5 Waveform of Direct Wave

[V]

[V]

Fig. 2 Assumption of Sample Composition

÷ 𝑒 −𝛼0 𝐿𝑇𝑅
t[us]
t[us]
(a) Direct Wave(b) Waveform of 𝜓𝑖,0
(Amplitude: 𝑃𝐷 [𝑉])(Estimation of 𝑝i,0 (𝑡))Fig. 6 Estimation of Emitted Wave 𝑝𝑡,0Fig. 3 Calibration
Environment

amplitudes of directed waves that are observed from a
few different distances 𝐿 𝑇𝑅 [m] from the transmitter.
When 𝑃𝑜 and 𝛼𝑜 are obtained, as Fig. 6 shows,cancellation of the attenuation by 𝛼𝑜 in direct wave Fig.
6(a) results in the solution of 𝑝𝑖,0 (𝑡), that is, 𝜓𝑖,0 (𝑡) (Fig.
6(b)).

Fig. 4 Process of Analyzing 𝑝𝑟

We limit searching area to find maximum of thecorrelation to 𝑡𝑘 > 𝑡𝑘−1 (Here, 𝑡𝑘 shows estimated time
when 𝑝𝑟,𝑘 (𝑡) is received). At this stage, 𝑡𝑘 should
overwrap the maximum amplitude of 𝜓𝑖,0 (𝑡).)In the estimation of 𝑝̂𝑟,𝑘 (𝑡), 𝛼𝑘 and 𝑟𝑘 are solved so thatthe squared error between 𝑝̂𝑟,𝑘 (𝑡) and 𝜓𝑟,𝑘 (𝑡) is
minimized. 𝑝̂𝑟,𝑘 (𝑡) has another unknown parameter 𝑐𝑘[m/s] (sonic velocity in medium k), which is obtained
through the fitting of linear function (the relation between
sonic velocity 𝑐𝑘 and acoustic impedance 𝜁𝑘 , which is
mentioned in Ref. 11). Hereafter, the estimation of 𝛼𝑘and 𝑟𝑘−1,𝑘 continues as 𝑝𝑟 (𝑡) is iteratively subtracted by𝜓𝑟,𝑘 (𝑡) until amplitude of 𝑝𝑟 (𝑡) drops below 𝑇𝑃 .

3.3. Expression of the reflected wave
As we have mentioned in section 3.1, 𝑝𝑟,𝑘 (𝑡) is
indicated in two ways, that is, 𝑝̂𝑟,𝑘 (𝑡) and 𝜓𝑟,𝑘 (𝑡), which
is made from the compression or stretch of 𝜓𝑖,0 (𝑡).
As 𝑝̂𝑟,𝑘 (𝑡), we consider the attenuation of ultrasonicwave originates from the reflection (The amplitude of thereflected wave narrows to the ratio 𝑟𝑘−1,𝑘 at location z=𝑧𝑘 ) and transmission (The transmitted wave travels in
velocity 𝑐𝑘 [m/s] and the amplitude narrows to the ratio of𝑒 −𝛼𝑘 (𝑧𝑘+1−𝑧𝑘) at 𝑧 = 𝑧𝑘+1 since its entry at 𝑧 = 𝑧𝑘 ).
After the emission of 𝑝𝑖,0 (𝑡) , the reflection andtransmission are repeated 2k-1, 2k times, each other.
Thus, 𝑝̂𝑟,𝑘 (𝑡) is expressed in Eq. (3), (4).

3.2. Calibration of attenuation coefficient andsound pressure of the emission−2𝛼𝑘−1 (𝑧𝑘 −𝑧𝑘−1 )𝑝̂.
𝑟,𝑘 (𝑡) = 𝑟𝑘−1,𝑘 𝐶Π,k 𝜓𝑖,0 (𝑡)𝑒𝑘−1−2𝛼𝑘′−1 (𝑧𝑘′ −𝑧𝑘′ −1 )2𝐶Π,k = ∏𝑘 ′ =1(1 − 𝑟𝑘 ′ −1,𝑘 ′ )𝑒.

Sound Pressure 𝑃𝑜 cannot be estimated in [Pa], then we
assume that the pressure is proportional to output voltage[V] of the receiver. In calibration process, we find 𝑃𝑜[V] ,𝛼𝑜 [Np/m] to use in the analysis hereafter.

(3)
(4)

However, 𝑝̂𝑟,𝑘 (𝑡) only expresses the relationship of two
voltage levels between emission and reflection at t[sec]345

As Eq. (6) shows, voltage level 𝑃𝐷 of attenuates sharply
as 𝐿 𝑇𝑅 increases near the origin of the emission, thus, weselected 𝐿 𝑇𝑅 as short distance as possible (=10. . ,50[𝑚𝑚])The attenuation by 𝛼0 is cancelled as Eq.
(7) shows.
𝑃0 = 𝑃𝐷 (𝐿 𝑇𝑅,1 )𝑒 𝛼0 𝐿𝑇𝑅,1 .
(7)

under attenuation. Though the Function Generator
outputs a cycle of sine wave, the direct wave is viewed asa group of waves with continuously varying amplitude,
and frequency 𝑓𝑐 [Hz] (same as central frequency of thetransducer). For the reason, we form 𝜓𝑟,𝑘 (𝑡) bystretching 𝜓𝑖,0 (𝑡) 𝐴𝑘 times vertically and 1/ 𝐵𝑘 times
horizontally so that the correlation between 𝑝̂𝑟,𝑘 (𝑡) and𝑝𝑟 (𝑡) is maximized. Here, we consider the attenuation of
amplitude and slight change of frequency after thereflection as Eq. (5) shows.
𝑝̂
𝑟,𝑘 (𝑡) = 𝐴𝑘 𝜓𝑖,0 (

𝑡−𝑡0
𝐵𝑘

)

4.2. Environment
The environment for experiment and the samples forreflection estimation is shown in Fig. 7, Fig. 8(a), (b),(c), (d), (e), each other. Air was removed from the liquid
to facilitate ultrasonic transmission. Dents and slightsteep were seen in the sample surfaces. The truth of𝜁 and 𝑐 is summarized in Table 2. Urethane Gel is madethrough(5)

Maximum sound pressure is assumed to be measured at
𝑡 = 𝑡0 in case transmitter and receiver contact each other.
4. Experiment
As we mentioned in Chapter 1, hardness s is regarded as
acoustic impedance, then 𝜁𝑘 is considered to be uniform
because the hardness is even in medium k, here, s is basedon Asker Type C [12] and labelled on the package ofpurchases [13]. We regard a sample as a medium in theexperiment. Now finding complicate hardness changewith ultrasonic is difficult, thus, we simplify hardness
heterogeneity with two different media of different
hardness, jointed vertically.
The object of experiment is to configure the increase of
𝑟1,2 as 𝜁2 is far from 𝜁1 regardless of different sample
depth 𝐻𝑆 (𝐻𝑆1 , 𝐻𝑆2 for sample 1, 2 each other). Though𝑟𝑘−1,𝑘 is negative when 𝜁𝑘 is lower than 𝜁𝑘−1 ,we only
consider positive reflection coefficient. We moved thesample randomly so that the center of the sample is near
the center between the receiver and transmitter. The
displacement and ultrasonic emission/analysis wererepeated alternatively for fifty times.

(a)Al 𝐻𝑆 = 10[𝑚𝑚](b)s=7 HS=36[mm] (c)s=0 HS=35[mm]Fig. 7 Experimental
Environment

Table 2 True Value of Sonic Speed c and Acoustic
Impedance  (at Temperature 30[deg])
Sample Name
HS [mm]c [m/s]ζ [Pa s/m]

4.1. Calibration
To find 𝑟0,1 from the acoustic model, 𝑃𝑜 and 𝛼𝑜 are
required. We recorded amplitude of direct wave 𝑃𝐷 using
the environment shown in Fig. 3 and changing 𝐿 𝑇𝑅 . We
expect that 𝑃𝐷 attenuates during the transmission
between the 𝐿 𝑇𝑅 [m] of interval as Eq. (6) shows. Under,
𝐿 𝑇𝑅 = 𝐿 𝑇𝑅,1 in the equation, 𝛼0 and 𝑃0 are iteratively
estimated for 20 times.
𝑃𝐷 (𝐿 𝑇𝑅 ) = 𝑃𝐷 (𝐿 𝑇𝑅,1 )𝑒 −𝛼0 (𝐿𝑇𝑅−𝐿𝑇𝑅,1) .

(d)s=7 HS=30[mm] (e)s=0 HS=29[mm]Fig. 8 Samples for the Environment

Al

10

6420

17.3x106

Gel
(s=0)

29
35

2237
2233

2.5x106
2.4x106

Gel
(s=7)

30
36

2185
1643

2.5x106
1.8x106

𝑇𝑃

Table 3 Setting for the Experiment
𝑃0 /100𝑓𝑐
200[kHz]𝐻𝐸𝜃1

50[mm]
15[deg]

(6)

346

𝑃0

±40V
±50V

3.07[V]
3.27[V]

𝛼0

±40V
±50V

2.43[Np/m]
1.52[Np/m]

τ±40V (Calib.)21[deg](Ref.)21[deg]±50V (Calib.)21[deg](Ref.)24[deg]𝑎1𝑎2

471.1[m /(Pa s2/m)]
349.8[m/s]transmitter into ±50[V] to increase sound pressure. The
result of 𝑟1,2 are shown in Fig. 12(pattern 1,2), andTable 4 Combination of Samples
PatternSample 1
Sample
HS1
Name
[mm]

Sample 2
Sample
HS2
Name
[mm]

Gel
(s=0)
Gel
(s=0)

29

Gel
(s=7)
Al

30

3

Gel
(s=0)

35

Gel
(s=7)

36

4

Gel
(s=0)

35

Al

10

1
2

29

𝑟0,1
10

𝑟0,1

𝑟1,2
𝑟1,2

Fig. 10 Result of 𝑟𝑘−1,𝑘 Estimation
for Pattern 1, 2 (Input:±40[V])

Isolated

𝑟0,1

Fig. 11 Result of 𝑟𝑘−1,𝑘 Estimation
for Pattern 3, 4 (Input:±40[V])

Outlier
𝑟0,1

𝑟1,2

𝑟1,2
Fig. 9 Estimation of 𝑐 with input of ζ

Fig. 12 Result of 𝑟𝑘−1,𝑘 Estimation Fig. 13 Result of 𝑟𝑘−1,𝑘 Estimation
for Pattern 1, 2 (Input:±50[V])for Pattern 3, 4 (Input:±50[V])the mixture and coagulation of two different types ofliquid. We solved true c by measuring elapsed time fromwave emission to reception and dividing 𝐻𝑆 with thetime. As the true 𝜁, we substituted c in Eq. (1) to obtain.
However, as 𝜁 of aluminum, we referred to Ref. 11
because any metal has extremely large 𝜁. We repeatedthe calibration 20 times and averaged the voltage of thereception to find 𝛼0 and 𝑃0 . Table 3 shows the setting of
this experiment (central frequency 𝑓𝑐 , distance fromtransmitter and sample surface 𝐻𝐸 , stilt of transducers
𝜃1 , temperature τ ), and estimation of 𝛼0 , 𝑃0 . The
amplitude of the input to transmitter was ±40,±50[V].
To find sonic velocity in media, we found simplified
relationship of linear function between 𝜁 and c, which isexpressed with coefficients 𝑎1 , 𝑎2 (Fig.9), 𝑎1 differsfrom ρ in that 𝑎1 is regarded uniform, while ρsynchronously changes as c changes. Table 4 shows thepatterns of the sample combination, that is, material
name, s, 𝐻𝑆 .

Table 5 Summary of Reflection Coefficient Estimation Result

Average of
𝑟0,1(50 times)
Standard
Deviationof 𝑟0,1Average of
𝑟1,2(50 times)
Standard
Deviationof 𝑟1,24.3. Result and Discussion
The comparison of 𝑟1,2 for pattern 1 and 2 is shown as
Fig. 10, for pattern 3 and 4, as Fig. 11, which wererecorded under condition that ±40[V] is input totransmitter. After that, we changed the amplitude of theInput±40[V]

Input±50[V]

Pattern1

0.60

0.71

Pattern2
Pattern3

0.79
0.36

0.76
0.56

Pattern4
Pattern1

0.45
0.05

0.64
0.04

Pattern2
Pattern3

0.05
0.05

0.05
0.06

Pattern4
Pattern1

0.05
0.22

0.08
0.29

Pattern2
Pattern3

0.37
0.15

0.32
0.24

Pattern4
Pattern1

0.18
0.05

0.36
0.04

Pattern2
Pattern3

0.06
0.08

0.10
0.06

Pattern4

0.03

0.05

Fig. 13 (pattern 3,4). True times of reflected wave
reception are also shown in Fig. 10-13, which were347

calculated with 𝐻𝑆 and c. Table 5 shows the estimation
result (average and standard deviation of 𝑟0,1 , 𝑟1,2 ). The
result of 𝑟0,1 , 𝑟1,2 is surrounded by a blue, red oval ineach figure (Fig. 10-13). As Table 5 shows, pattern 2shows higher 𝑟1,2 than pattern 1 and, pattern 4 showshigher 𝑟1,2 than pattern 3, which means, increase of
𝜁2 resulted in rising 𝑟1,2 . However, note that the averageof 𝑟1,2 is calculated without outlier in Fig. 12. Also, forFig. 12, the distribution of 𝑟1,2 was isolated (shown withan orange circle) apart. The 30[us] of delay in 𝑡1 wasshown in pattern4 of Fig. 13, which was mainly causedby heterogeneity of sample depth. To reduce the timedifference between measured value and truth, bump orinclination of sample surface should be eliminated undercoagulation.

8. J. Carlson, J. Deventer, A. Scolan and C. Carlander,“Frequency and Temperature Dependence of Acoustic
Properties of Polymers Used in Pulse-Echo Systems“,Proceedings in IEEE Symposium on Ultrasonics, pp. 885888, 2003.
9. S. Takeuchi, Tutorial on the underwater or medical
ultrasound transducer, (In Japanese) J. Acoust. Soc. ofJapan, Vol. 72, No. 5, pp. 264-272, 2016.
10. R. K. Young, N. K. Bose, “Wavelet Theory and Its
Applica-tions”, Kluwer Academic Publishers,1993.
11. Honda Electronics Co., Ltd. “Ultrasound Handbook,” (InJapanese), 2008.
12. “Kobunshi Keiki Co., Ltd. ASKER Durometer Type C”.
13. “Exseal Co., Ltd. Hitohada no Gel Gen’eki” (in Japanese).

Authors Introduction
Dr. Ryuugo Mochizuki
He is a researcher at Graduate
School of Life Science andSystem Engineering, Kyushu
Institute of Technology, Japan.
His research topic during thePhDcoursewas
image
processing. Currently, he isinvolved in signal processing of
ultrasonic.

5. Conclusion
We evaluated the effect of changing bottom media ofsamples into harder one (from urethane gel to aluminum).
As the result, we succeeded in finding larger 𝑟1,2
according to the increase of 𝜁2 , which may show thepossibility of estimating the variety of hardness inside
samples, however, the time of wave reception was notstable in the analysis.

Dr. Yuya Nishida
He is an Associate Professor at
Graduate School of Life Scienceand System Engineering, Kyushu
Institute of Technology, Japan.
His research area is about filed
robotics, its application, and dataprocessing.

References
1. “The Japan Institute for Labour Policy and Training, “(in Japanese)
2. “Yahoo News” https://news.yahoo.co.jp/byline/
fuwaraizo/20130909-00027589 (in Japanese)
3. N. Kinoshita, “The Challenge for Sodium Reduction by
Family Mart: Promotion of Unnoticed Sodium Reduction”,(In Japanese), Japanese Society of Health Education andPromotion, Vol. 29, No. 3, 2021.
4. J Machado, P. Palma, S. Simoes, “Ultrasonic Indirect
Method for Evaluating Clear Wood Strength andStiffness”, Non-Destructive Testing in Civil Engineering,pp. 1-6, 2009,
5. B. Cho and J. Irudayararaj, “Design and application of anon-contact ultrasound velocity measurement system with
air instability compensation”, Transactions of theAmerican Society of Agricultural Engineers, Vol.46, No.
9, pp. 901-909, 2003.
6. B. Cho and J. Irudayararaj, “A Noncontact UltrasoundApproach for Mechanical Property Determination ofCheeses”, Food Engineering and Properties, Vol.68, No. 7,
pp. 2243-2247, 2003.
7. S. Srivastava, S. Vaddadi and S. Sadistap, “Non-Contact
Ultrasonic Based Stiffness Evaluation System forTomatoes during Shelf-Life Storage”, Journal of Nutrition
and Food Sciences, Vol. 4, Issue 3, pp 1-6, 2014.

Prof. Kazuo Ishii
He is a Professor at Graduate
School of Life Science andSystem Engineering, Kyushu
Institute of Technology, Japan.
His research area is about field
robots and intelligent robot
system348