Strategy of tactile estimation modeling
The structure of this study is shown in Figure 1 to demonstrate the effectiveness of the nonlinear modeling of tactile estimation. First, we conducted a sensory evaluation test. In parallel, a tactile sensing system is developed, including a newly designed tactile sensor. The vibration data while the sensor traces on a sample are acquired, followed by feature extraction based on the perceptibility of mechanoreceptors. Then, we develop several model candidates with or without nonlinearities, whose inputs and outputs are the extracted features and sensory evaluation scores, respectively. The mean squared error between the model predictions and the actual values is compared to determine the most effective model.
Target samples
As shown in Figure 2 (a), eight plastic plates were used in this study as target samples. The plates had different surface patterns. The differences are shown in Table S1 as the arithmetic average roughness of the samples, Ra, and the arithmetic average swell, Wa, measured using DektakXT (Bruker Corporation, Billerica, MA, USA). The dynamic friction coefficient, \(\mu \text{'}\), was measured using KESSE with a 10 mm2 pianowire sensor (Kato Tech Co. Ltd., Kyoto, Japan), as shown in Figure 2 (b).
Sensory evaluation test
To quantify the tactile sensations, a sensory evaluation test was conducted with human subjects, performed under 24.6 ± 1.0 ℃ and 41.3 ± 4.5% relative humidity with the participation of 35 healthy adults (21 males and 14 females) aged 22.2 ± 1.1 (between 21 and 25) years old. We employed a semantic differential method with a sevenstep unipolar scale for the Japanese adjectives listed in Table S2. Before the evaluation, the subjects were free to touch all the samples to understand the variety of the samples. During the test, each sample was placed in a box to exclude visual information. The subjects were allowed to evaluate the samples in random order. In addition, the evaluation words were given to each participant in random order to prevent a possible order effect from occurring. The test protocol was approved by the Bioethics Board of the Faculty of Science and Technology, Keio University. The subjects received a thorough explanation of the test methods in advance and then signed an informed consent form before participating in the study.
To identify trends in the subjects’ responses, the subjects were classified using cluster analysis using the Python library, Scipy. The scores for all evaluated words were considered in the cluster analysis. The Euclidean distance was used as the distance function, and the Ward method was employed.
For the analysis of tactile sensation, the evaluation words were analyzed by principal component analysis (PCA) for each cluster of subjects based on the evaluation scores using SPSS (Version 22, International Business Machines Corp., Armonk, NY, USA). The conditions for extracting the principal components (PCs) include the criteria that the eigenvalue of each PC should be greater than unity.
Vibration measurement system and procedure
We developed a tactile sensing system capable of detecting vibrations while a tactile sensor runs over a sample. Figure 3 shows the images of the tactile sensor and an overview of the tactile sensing system. A cylindrical shaft converts vertical displacement into the strain of a leaf spring with two strain gauges (KFGS03120C123N30C2, Kyowa Electronic Instruments Co. Ltd., Tokyo, Japan) glued on both sides to detect vibrations when a silicone rubber pad traces a sample surface. The hardness of the silicone rubber pad was designed to be equivalent to that of a human finger. Figure S1 shows a comparison of hardness between the forefinger pad and the developed sensor measured using a durometer TYPE OO (GS754G, Teclock Co. Ltd., Nagano, Japan). The results of the Student ttest showed no significant differences between the two.
A coating material (X9317551, ShinEtsu Chemical Co. Ltd., Tokyo, Japan) was adhered to the surface of the silicone rubber pad. The outputs from the strain gauges were acquired using a dynamic strain amplifier (DPM913B, Kyowa Electronic Instruments Co. Ltd., Tokyo, Japan). The relationship between the output of the strain gauge (V) and the vertical deformation of the tactile sensor (d) is shown in Figure S2. From the figure, we can obtain the transformation equation with linear regression as
$$\begin{array}{c}d= 1409.4V+710.35\#\left(1\right)\end{array}$$
As shown in Figure 3 (c), the tactile sensor was fixed to a traction arm of the static and dynamic friction measuring instruments (TL201Ts, TrinityLab. Inc., Tokyo, Japan). Upon testing, the normal force N between the tactile sensor and the sample can be adjusted by placing weights on the traction arm. As the sample table of the TL201Ts moves horizontally, the tactile sensor runs over a sample. In addition, a force sensor connected to the traction arm detects the tangential force F.
The vibration information measurement conditions were as follows: the tracing speed and distance of the tactile sensor were 10 mm/s and 30 mm, respectively. The normal force, N, applied between the tactile sensor and sample was 0.49 N. Measurements were repeated 11 times per sample.
Data processing methods
There are four mechanoreceptors in the glabrous skin, as shown in Figure S3: Meissner corpuscles (FA I), Pacinian corpuscles (FA II), Merkel disks (SA I), and Ruffini endings (SA II) [2224]. They respond to mechanical stimuli applied to the skin and then fire nerve impulses to the neuron. The physiological threshold of the amplitude of stimulation for firing against the frequency has been reported for each receptor, as summarized in Figure S4. Based on this, we can approximate the threshold line, L, on the logarithmic chart for each mechanoreceptor as
where \({L}_{FAI}\), \({L}_{FAII}\), \({L}_{SAI}\), and \({L}_{SAII}\) are the thresholds for FA I, FA II, SA I, and SAII, respectively, and f is the frequency of the vibration stimulus. Each mechanoreceptor fires when the intensity of the mechanical stimulus surpasses the corresponding threshold line.
The vibration data acquired by the vibration measurement system (Figure 3) were compared with the abovementioned threshold lines in the frequency domain to extract meaningful information for tactile estimation as follows: First, the acquired output from the strain gauges was converted to vertical displacement using Equation (1), resulting in the vibration data in the time domain. Then, we transformed the vibration data from the time domain to an amplitude spectrum in the frequency domain using fast Fourier transformation (FFT), implemented in MATLAB (MATLAB 2020a, Math Works Inc., Natick, MA, USA) at a sampling frequency of 10 kHz mediated with a Hamming window, to obtain the vibration data in the frequency domain. Figure 4 shows a conceptual diagram of the extraction of feature values. The colored area between the lowest threshold and the measured vibration data corresponds to the firing status of the mechanoreceptors. This study considers the combination of firing receptors and divides the entire area into eight subareas, Di, as shown in Figure 4. The subscript i represents the mechanoreceptors that are supposed to fire in the corresponding frequency range. Note that Di could be zero when the vibration data are always lower than the thresholds. Detailed formulas for calculating each feature are provided in the Supplementary Material.
Tactile estimation models
We performed a regression analysis to predict the principal component scores for each cluster using the features extracted from the vibration data, Di, and the dynamic friction coefficient, \(\mu \text{'}\), using the Python library and state models. Considering that the human tactile perception nature has nonlinearity [25], we developed four types of linear/nonlinear regression models:
$$\begin{array}{c}Linear :y={\beta }_{0}+{\sum }_{i=1}^{p}{\beta }_{i}{x}_{i}\left(6\right)\end{array}$$
$$\begin{array}{c}Logarithmic :y={\beta }_{0}+{\sum }_{i=1}^{p}{\beta }_{i}\text{log}\left({x}_{i}\right)\left(7\right)\end{array}$$
$$\begin{array}{c}Interaction :y={\beta }_{0}+{\beta }_{1}{x}_{i}+{\beta }_{2}{x}_{j}+{\beta }_{3}{x}_{i}{x}_{j}\left(8\right)\end{array}$$
$$\begin{array}{c}Polynomial :y={\beta }_{0}+{\sum }_{i=1}^{a}{\beta }_{1}{x}_{1}^{i}\left(9\right)\end{array}$$
where xi and xj are the explanatory variables, that is, Di and µ’. y is the objective variable, that is, the PC score. \({\beta }_{i}\) represents the coefficients to be determined. In the interaction model, any two explanatory variables were chosen to build the model, that is, 6C2 = 15 types of models were built for one objective variable.
The explanatory variables were introduced using the bruteforce method. A regression model was constructed using data obtained from seven out of eight samples. The data for the remaining samples were used to validate the developed regression model. This process was repeated eight times, that is, all samples were used for validation. The model with the lowest error was selected as the best model for each PC in each cluster.
For comparison, we also conducted regression analyses using the feature extraction methods reported in a previous study [17]. In other words, a total of eight regression equations were constructed for one objective value by combining two types of feature extraction methods and four types of models, and the name of each regression model was defined, as shown in Table 1. Previous research [17] only considered three features extracted from the vibration data and used linear regression analysis. Thus, A1 is the method reported in a previous study [17], and the other models were newly constructed in this study.
Table 1
Classification of the regression models based upon the feature extraction method and model type

Feature extraction method

Previously reported
method [17]

Proposed method

Model type

Linear

A1

B1

Logarithmic

A2

B2

Interaction

A3

B3

Polynomial

A4

B4
