Short-Time Noise Suppression Method with Multiple Signals for Received Time Estimation Masanari Nakamura1 , Akira Tanaka1 , Hiromichi Hashizume2 and Masanori Sugimoto1 1 Hokkaido University, Kita 8, Nishi 5, Kita-ku, Sapporo, Hokkaido, 060-0808 Japan 2 National Institution for Academic Degrees and Quality Enhancement of Higher Education, 1-29-1 Gakuen-nishimachi Kodaira-shi, Tokyo 187-8587 Japan Abstract The detection of the leading edge of a signal using a threshold value can be used as a method for estimating the received time, which is not easily affected by reverberation or the sensor characteristics of the speaker and microphone. When the noise level is relatively high compared with the signal, the threshold must be set high to prevent false alarms, causing a systematic error in the received time estimation. The noise can be suppressed by transmitting the signal multiple times and averaging the results; however, it is necessary to wait for the reverberation to disappear after each transmission. Therefore, this measurement is time consuming. In this study, we propose a noise suppression method that can be applied to cases in which signals are transmitted multiple times without waiting for the disappearance of reverberation. For the proposed method, an analysis of noise suppression performance was conducted, and the limits of performance were provided under the given conditions. In addition, numerical experiments on the performance of the received time estimation in a reverberant environment were conducted to compare the proposed method with the conventional method and confirm the effectiveness of the proposed method. Keywords received time estimation, threshold, false alarm, reverberation, sensor characteristics 1. Introduction With the recent spread of mobile devices, indoor positioning techniques for these devices have been extensively studied[1]. In this study, we discuss a positioning method that uses the propagation time of acoustic signals. When using the propagation time, it is important to estimate the received time of the acoustic signal. In indoor environments, reverberation and the physical characteristics of speakers and microphones can cause large systematic errors in the estimation of received time. The detection of the leading edge of the signal with a threshold value can be used as a method that is less sensitive to reverberation and sensor characteristics. In detection using thresholds, it is necessary to set a threshold based on the noise level to avoid false noise detection. When the signal-to-noise ratio is low, the threshold must be set high, the parts of the leading edge of the signal are delayed, and a systematic error is caused in the estimation of the received IPIN 2022 WiP Proceedings, September 5 - 7, 2022, Beijing, China $ masanari@ist.hokudai.ac.jp (M. Nakamura) Β© 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) time. Therefore, to reduce systematic errors with the desired probability of a false alarm, the signal-to-noise ratio must be improved. If the noise follows a white normal distribution, it can be suppressed by transmitting the signal multiple times and averaging them. However, this method requires waiting for the disappearance of the reverberation after each transmission. Therefore, it takes a long time to send all the signals. This could be a problem, especially when positioning is required at many locations. In this study, we propose a noise suppression method that can be applied even when the transmission time is set shorter than the reverberation time to reduce the time required for signal transmission. By applying an orthogonal projection matrix that preserves the transmitted signal component of the received signal, only the noise component is suppressed. The contributions of this study are as follows: β€’ We devised a noise suppression method that can be applied when the transmission period is set to less than the reverberation time. β€’ The relationship between the transmission period and the number of transmissions and the noise suppression performance of the proposed method is discussed. The limitations of the noise suppression performance of the proposed method under given conditions are also provided. β€’ We conducted numerical experiments on the performance of received time estimation un- der multiple reverberation conditions, and demonstrated the effectiveness of the proposed method over the conventional method. In addition, we conducted numerical experiments regarding the performance of received time estimation when the required transmission time is changed, and discussed their relationship. The remainder of this study is structured as follows. In Section 2, we describe conventional methods for estimating the received time. In Section 3, we describe the proposed method and discuss the relationship between the transmission period, number of transmissions, and noise suppression performance. In Section 4, we discuss the results of the experiments comparing the proposed method with conventional methods. Finally, in Section 5, we present our conclusions. 2. Related work 2.1. Observation model First, we described the observational model discussed in this study. The start times of the transmission and the transmission signal are set to 0 and 𝑠(𝑑), respectively. The signal length was 𝑇 . A signal transmitted through a speaker is nonlinearly distorted, depending on the physical properties of the speaker. This is expressed as follows: 𝑓 (𝑑) = Ξ˜π‘  (𝑠(𝑑)) (1) 𝑓 (𝑑) and Ξ˜π‘  (Β·) represent the signal distorted by the speaker and the characteristics of the speaker that cause distortion, respectively. In indoor environments, reverberation is caused by multiple delayed waves reflected on walls, ceilings, and other surfaces. When this signal is observed using a microphone, nonlinear distortion occurs at the microphone. Let β„Ž(𝑑) and Ξ˜π‘Ÿ (Β·) be the indoor impulse response and microphone characteristics that cause distortions, respectively. The reverberation signal can be expressed as follows: 𝑔(𝑑) = Ξ˜π‘Ÿ ((β„Ž * 𝑓 )(𝑑)) (2) * is a convolution. Let Δ𝑑 be the time required for propagation between the speaker and microphone. The received signal can be expressed as: π‘Ÿ(𝑑) = 𝑔(𝑑 βˆ’ Δ𝑑 βˆ’ 𝛿) + 𝑛(𝑑) (3) 𝛿 represents the clock offset between the transmitter and receiver. 𝑛(𝑑) is white normal noise. If 𝛿 is known, the propagation time Δ𝑑 can be estimated from the received time. The distance can be estimated by multiplying Δ𝑑 by the propagation velocity. Even if 𝛿 is not known, 𝛿 can be eliminated if the transmission times of multiple speakers are known. Therefore, the distance between the speaker and the microphone can be estimated. 2.2. Received time estimation by maximum peak detection The matched filter is a well-known method for estimating the received time in noise[2]. In the following, we explain the matched filter in the case where the transmission signal 𝑠(𝑑) is a sinusoidal wave, as an example. {οΈƒ sin(2πœ‹π‘“ 𝑑) 0 ≀ 𝑑 ≀ 𝑇 𝑠(𝑑) = (4) 0 π‘œπ‘‘β„Žπ‘’π‘Ÿπ‘€π‘–π‘ π‘’ 𝑓 is frequency, which is a constant multiple of 1/𝑇 . The matched filter output is expressed as 1 ∞ ∫︁ 𝑐(𝑑) = π‘Ÿ(𝜏 )𝑒(𝜏 βˆ’ 𝑑)π‘‘πœ, (5) 𝑇 βˆ’βˆž {οΈƒ exp(𝑗2πœ‹π‘“ 𝑑) 0 ≀ 𝑑 ≀ 𝑇 𝑒(𝑑) = (6) 0 π‘œπ‘‘β„Žπ‘’π‘Ÿπ‘€π‘–π‘ π‘’ where 𝑒(𝑑) denotes the reference signal. For simplicity, 𝛿 is assumed as zero in this study. If there is no influence of reverberation or sensor characteristics, the absolute value of 𝑐(𝑑) reaches its maximum value, where 𝑑 = Δ𝑑. Therefore, in this case, the received time must be the time when the absolute value of 𝑐(𝑑) is largest. In this paper, this method is referred to as the maximum peak method. The maximum peak method is optimal when there is only white normal noise and no reverberation or sensor characteristic influence[2], and it is used in various acoustic positioning methods[3, 4, 5, 6, 7, 8, 9]. In a reverberant environment, as expressed in Eq. (5), the reference signal also responds to the delayed waves. Therefore, there is an error in the received time estimation with the maximum peak method when the delayed wave is superimposed on the direct wave or when the delayed wave is larger than the direct wave. The maximum peak method assumes that the direct waves of the received and reference signals are coincident. However, in practice, because distortions occur depending on the sensor (Eq. (1)), it has been reported that systematic errors can occur [10]. One way to avoid both the reverberation and distortion of the received signal is to measure 𝑔(𝑑) in advance and use this signal as a reference. However, because the delayed waves that constitute reverberation and the sensor characteristics vary depending on the location of the speaker and microphone, a large amount of premeasurement is required. 2.3. Received time estimation by leading-edge detection 2.3.1. False alarm probability and threshold To avoid the effects of reverberation and sensor characteristics, detection of the leading edge of the signal can be used to estimate the received time. Specifically, a threshold was set for the absolute value of the matched filter output, and the first point at which it exceeded the threshold was detected. From the matched filter equation (Eq. (5)), the time of the detected point is 𝑇 earlier than the received time. Therefore, the received time was the sum of the detected time and 𝑇 . Hereafter, this method is referred to as the leading-edge method. In this section, we discuss the probability of false noise detection. This is called the false-alarm probability. When a matched filter is applied to white normal noise, the probability density distribution of the absolute value of its output is a Rayleigh distribution. We describe the Rayleigh distribution when a discretized matched filter is used because the process with discrete data is discussed in Section 3 and later. Let 𝑁 denote the signal length 𝑇 converted to the number of samples. Let 𝑇𝑠 denote the sampling period, and let 𝑇 be a consist multiple of 𝑇𝑠 . In this case, 𝑁 = 𝑇 /𝑇𝑠 . Let 𝜏 be the threshold; then, the false alarm probability 𝑃𝐹 𝐴 is calculated as follows: ∫︁ ∞ π‘₯2 (οΈ‚ )οΈ‚ π‘₯ 𝑃𝐹 𝐴 = exp βˆ’ 2 𝑑π‘₯ (7) 𝜏 𝐸2 2𝐸 The Rayleigh parameter 𝐸 is computed using the power spectrum 𝜎 2 of the white normal noise as 𝐸 2 = 𝜎 2 /2. Using the cumulative distribution function of the Rayleigh distribution, the above equation can be expressed as 𝑃𝐹 𝐴 = exp(βˆ’πœ 2 /(2𝐸 2 )). Solving this for πœβˆšοΈ€ provides a threshold that allows signal detection at a specified false-alarm probability as 𝜏 = βˆ’2𝐸 2 log(𝑃𝐹 𝐴 ). As can be observed from these equations, setting a large threshold 𝜏 reduces the false-alarm probability. However, a larger threshold value caused a greater delay in the estimated received time. Therefore, to reduce the systematic error in the received time with a sufficiently small false-alarm probability, the signal-to-noise ratio of the received signal must be improved. 2.3.2. Improvement of signal-to-noise ratio by averaging Let 𝑛1 , 𝑛2 , . . . , 𝑛𝑍 be samples that follow βˆ‘οΈ€ a white normal distribution with mean 0 and vari- ance 𝜎 2 , and the variance of the average 𝑍 𝑖=1 𝑛𝑖 /𝑍 is 𝜎 /𝑍. Therefore, if the noise of the 2 received signal follows a white normal distribution, it can be suppressed by averaging multiple transmitted signals. Hereafter, this method will be referred to as the averaging method. In reverberant environments, if the averaging method is used with the transmission period set to less than the reverberation time, the transmitted signal components will be corrupted during averaging by the delay waves of reverberation. In this case, the received time cannot be properly estimated. Therefore, to improve the signal-to-noise ratio, the signal transmission period should be longer than the reverberation time. However, this requires a significant amount of time. In the next section, we propose a noise suppression method that can improve the signal-to-noise ratio even when the transmission period is set below the reverberation time to reduce the required transmission time. 3. Proposed method 3.1. Transmission and reception model In the proposed method, 𝑠(𝑑) is transmitted 𝐿 times in period 𝑇𝑑 . Let 𝑁𝑑 denote the signal transmission period 𝑇𝑑 converted to the number of samples. For simplicity, let 𝑇𝑑 be constant multiples of 𝑇𝑠 . In this case, 𝑁𝑑 = 𝑇𝑑 /𝑇𝑠 . Let π‘Ÿ be the discrete data of the signal that is received and sampled and let π‘π‘Ÿ be its length. Let π‘π‘Ÿ be sufficiently large. Subsequently, π‘Ÿ can be expressed as follows: ⎑ ⎀ π‘‚Ξ”π‘›βˆ’1 βˆ‘οΈπΏ ⎒ 𝑂𝑁𝑑 (π‘™βˆ’1) βŽ₯ (8) ⎒ βŽ₯ π‘Ÿ= ⎒ ⎒ 𝐼𝑀 βŽ₯ 𝑔 + π‘›π‘Ÿ . βŽ₯ 𝑙=1 ⎣ 𝑂𝑁𝑑 (πΏβˆ’π‘™) ⎦ π‘‚π‘π‘Ÿ βˆ’π‘€ βˆ’Ξ”π‘›+1βˆ’π‘π‘‘ (πΏβˆ’1) From the explanation in Section 2.3.1, Δ𝑛 is the index of the sample closest to Δ𝑑 βˆ’ 𝑇 and rounded to the nearest (Δ𝑑 βˆ’ 𝑇 )/𝑇𝑠 . Let 𝑔 = [𝑔1 𝑔2 . . . 𝑔𝑀 ]𝑇 be the sampling of Eq. (2) and let 𝑀 be its length. 𝑇 in the upper-right corner of the symbol indicates transposition. 𝐼𝑀 is the identity matrix of degree 𝑀 and 𝑂𝑖 denotes the zero matrix whose size is 𝑖 Γ— 𝑀 . π‘›π‘Ÿ is a vector of the noise components. Each element is assumed to follow a white normal distribution with a mean of zero and a variance 𝜎 2 . In the following, it is assumed that 𝜎 2 is known from the prior measurements. The noise suppression method detailed in Section 3.3 assumes that 𝑔1 is contained in the first 𝑁𝑑 elements of the vector to which the orthogonal projection matrix is applied. To satisfy this assumption, the signal was cut out of π‘Ÿ. Let 𝑦 denote the signal. Hereafter, the range from the first to the 𝑁𝑑 -th sample of 𝑦 is called the decision interval. The signal cutout procedure is described in detail. First, a discrete matched filter is applied √ to the received signal to determine the starting point of the signal cutout as 𝑐𝑖 = π‘Ÿπ‘–,𝑖+𝑁 𝑇 βˆ’1 𝑒/ 𝑁. Here, π‘Ÿπ‘–,𝑖+𝑁 βˆ’1 is the cutout of π‘Ÿ from the 𝑖-th to the (𝑖 + 𝑁 βˆ’ 1)-th. Let Δ𝑛(1) be the index that exceeds threshold 𝜏 (1) for the first time in |𝑐𝑖 |. Here, 𝜏 (1) is set (1) such that the false-alarm probability 𝑃𝐹 𝐴 is sufficiently small. Therefore, Δ𝑛(1) β‰₯ Δ𝑛. The starting point of the signal cutout is Δ𝑛(1) βˆ’ 𝑁𝑑 + 1, which allows 𝑔1 to be included in the decision interval. In this case, the cut-out signal 𝑦 can be expressed as follows: 𝑦 = 𝐴π‘₯ + 𝑛 (9) ⎑ ⎀ 𝑂 β€² 𝐿 βˆ‘οΈ ⎒ 𝑑 𝑁 (π‘™βˆ’1) 𝐴= ⎣ 𝐼𝑀 β€² ⎦ (10) βŽ₯ 𝑙=1 𝑂𝑁′ 𝑑 (πΏβˆ’π‘™) ]︀𝑇 π‘₯ = 0 𝑧1 𝑔 𝑇 0 𝑧2 (11) [οΈ€ π‘₯ denotes a vector with zeros before and after 𝑔. 𝑀 β€² denotes the length of π‘₯. 𝑂𝑖′ is the zero matrix whose size is 𝑖 Γ— 𝑀 β€² . 𝑧1 , 𝑧2 are 𝑁𝑑 βˆ’ 1 + Δ𝑛(1) βˆ’ Δ𝑛 and 𝑀 β€² βˆ’ 𝑀 βˆ’ 𝑁𝑑 + 1 βˆ’ Δ𝑛(1) + Δ𝑛, respectively. The cutout length is 𝑀 β€² + 𝑁𝑑 (𝐿 βˆ’ 1), which includes all transmitted signals and their reverberations. This equals the number of rows in 𝐴, and let 𝑅. 𝑛 is a part of π‘›π‘Ÿ that is cut out of π‘Ÿ with the signal component. Here, we discuss how to set the parameter 𝑀 β€² . The length 𝑀 of 𝑔 contained in π‘₯ is determined by Ξ˜π‘  (Β·), Ξ˜π‘Ÿ (Β·), and β„Ž(𝑑), and varies depending on the sensor characteristics and environment. Therefore, 𝑀 is assumed to be unknown. In addition, because the starting point of the cutout is Δ𝑛(1) βˆ’ 𝑁𝑑 + 1, the beginning of π‘₯ can contain fewer than 𝑁𝑑 zeros. Considering the above, 𝑀 β€² should be set to a value with some margin. 3.2. Noise suppression in the decision interval In the proposed method, the orthogonal projection matrix 𝑃𝐴 𝑃𝐴 = 𝐴(𝐴𝑇 𝐴)βˆ’1 𝐴𝑇 (12) was applied to 𝑦. To show that (𝐴𝑇 𝐴)βˆ’1 exists, it is sufficient to show that 𝐴𝑇 𝐴 has full rank, which is proved in the Appendix. The result 𝑦 β€² can be expressed as follows: 𝑦 β€² = 𝑃𝐴 𝑦 = 𝐴(𝐴𝑇 𝐴)βˆ’1 𝐴𝑇 𝐴π‘₯ + 𝑃𝐴 𝑛 = 𝐴π‘₯ + 𝑃𝐴 𝑛 (13) This equation shows that signal component 𝐴π‘₯ can be retained without corruption. For the noise component, let 𝑃𝐴 𝑛 be 𝑛′ , then, the 𝑖-th element of 𝑛′ can be expressed as 𝑛′𝑖 = 𝑝𝑇𝑖 𝑛 where 𝑝𝑖 refers to the 𝑖-th columns of 𝑃𝐴 . Because 𝑛′ is a linear sum of white normal noise, 𝑛′𝑖 follows a normal distribution, and its variance (πœŽπ‘–β€² )2 is (πœŽπ‘–β€² )2 = 𝑝𝑇𝑖 𝑝𝑖 𝜎 2 . This implies that 𝑝𝑇𝑖 𝑝𝑖 represents the noise suppression performance; if this value is less than one, a smaller threshold value can be set with the same false alarm probability. In this case, as discussed in section 2.3.1, the systematic error in the received time estimation can be reduced. The detailed process is (2) (2) as follows: The threshold value πœπ‘– is calculated from a given false-alarm probabilityβˆšπ‘ƒπΉ 𝐴 and (πœŽπ‘–β€² )2 . The absolute value of the discrete matched filter output 𝑐′𝑖 = (𝑦𝑖,𝑖+𝑁 β€² βˆ’1 ) 𝑒/ 𝑁 is 𝑇 (2) calculated, and detection is performed using the threshold πœπ‘– . 𝑦𝑖,𝑖+𝑁 β€² βˆ’1 is the cutout of 𝑦 β€² from the 𝑖-th to the (𝑖 + 𝑁 βˆ’ 1)-th. (𝑖) Let Δ𝑛(2) be the index of 𝑦 β€² first detected by the threshold 𝜏2 , then, the estimate of Δ𝑛 is expressed as Δ𝑛(1) + Δ𝑛(2) βˆ’ 𝑁𝑑 . 3.3. Relationship between transmission period, number of transmissions, and noise suppression performance In this section, we first show that to study the noise suppression performance of 𝑃𝐴 in the decision interval, it is sufficient to investigate the diagonal components of (𝐴𝑇 𝐴)βˆ’1 . Next, we show that the noise suppression performance can be expressed in a simple form using 𝐿, 𝑁𝑑 , and 𝑀 β€² under the given conditions, and we discuss the relationship between each parameter and the noise suppression performance. In addition, the performance limits of the proposed method under these conditions are discussed. 3.3.1. Analysis of 𝑃𝐴 𝑝𝑇𝑖 𝑝𝑖 representing the noise suppression performance, is the 𝑖-th diagonal component of 𝑃𝐴𝑇 𝑃𝐴 . Here, 𝑃𝐴 is a projective matrix and is symmetric. Hence, 𝑃𝐴𝑇 𝑃𝐴 = 𝑃𝐴 𝑃𝐴 = 𝑃𝐴 . Therefore, 𝑝𝑇𝑖 𝑝𝑖 is identical to the 𝑖-th diagonal component of 𝑃𝐴 . Here, 𝐴 can be expressed as 𝐴 = 𝑇 (𝐼1,𝑁 𝑑 𝐾0𝑇 )𝑇 where 𝐼1,𝑁𝑑 is the first to 𝑁𝑑 -th rows of the identity matrix of degree 𝑀 β€² , and 𝐾0 is the (𝑁𝑑 + 1)-th and the later rows of 𝐴. Therefore, 𝑃𝐴 can be written as: (οΈ‚ )οΈ‚ 𝐼 𝑃𝐴 = 1,𝑁𝑑 (𝐴𝑇 𝐴)βˆ’1 𝐼1,𝑁 (οΈ€ 𝑇 𝑇 )οΈ€ 𝐾 0 𝐾0 𝑑 (οΈ‚ 𝑇 βˆ’1 )οΈ‚ (𝐴 𝐴)1,𝑁𝑑 (οΈ€ 𝑇 𝐼1,𝑁𝑑 𝐾0𝑇 . (14) )οΈ€ = 𝑇 βˆ’1 𝐾0 (𝐴 𝐴) Thus, each element in the top-left 𝑁𝑑 -by-𝑁𝑑 matrix of 𝑃𝐴 and (𝐴𝑇 𝐴)βˆ’1 is identical. This indicates that the noise suppression performance in the decision interval can be evaluated by investigating the first to 𝑁𝑑 -th diagonal components of (𝐴𝑇 𝐴)βˆ’1 . 𝐴𝑇 𝐴 is a symmetric Toeplitz matrix, as follows: 𝑣1 𝑣2 𝑣3 Β· Β· Β· Β· Β· Β· 𝑣𝑀 β€² βŽ› ⎞ .. .. ⎟ . ⎜ ⎜ 𝑣2 𝑣1 𝑣2 . ⎟ .. ⎟ ⎜ ⎟ ⎜ .. .. .. ⎜ 𝑣3 𝑣2 . . . . ⎟ 𝑇 𝐴 𝐴=⎜ ⎜ . ⎟ (15) ⎜ .. .. .. .. ⎟ ⎜ . . . 𝑣2 𝑣3 ⎟ ⎟ ⎜ . .. ⎝ .. ⎟ . 𝑣2 𝑣1 𝑣2 ⎠ 𝑣𝑀 β€² Β· Β· Β· Β· Β· Β· 𝑣3 𝑣2 𝑣1 All the elements of this matrix are integers and are represented as follows: {οΈƒ 𝐿 βˆ’ π‘šβˆ’1 ((π‘š βˆ’ 1) mod 𝑁𝑑 ) = 0 ∧ 𝐿 > π‘šβˆ’1 π‘£π‘š = 𝑁𝑑 𝑁𝑑 . (16) 0 π‘œπ‘‘β„Žπ‘’π‘Ÿπ‘€π‘–π‘ π‘’ Note that the inverse matrix of symmetric Toeplitz matrix can be obtained by the Trench algorithm. We then examine the diagonal components of inverse matrix (𝐴𝑇 𝐴)βˆ’1 . Let each of these elements be defined as follows: βŽ› ⎞ 𝑀1,1 𝑀1,2 Β· Β· Β· 𝑀1,𝑀 β€² ⎜ 𝑀2,1 𝑀2,2 Β· Β· Β· 𝑀2,𝑀 β€² ⎟ (𝐴𝑇 𝐴)βˆ’1 = ⎜ . .. ⎟ . (17) ⎜ ⎟ . .. . . ⎝ . . . . ⎠ 𝑀𝑀 β€² ,1 𝑀𝑀 β€² ,2 Β· Β· Β· 𝑀𝑀 β€² ,𝑀 β€² This matrix satisfies (𝐴𝑇 𝐴)(𝐴𝑇 𝐴)βˆ’1 = 𝐼. We consider the first diagonal component 𝑀1,1 of (𝐴𝑇 𝐴)βˆ’1 . For the first column of (𝐴𝑇 𝐴)βˆ’1 , β€² we assume that 1 + 𝑁𝑑 (𝑖 βˆ’ 1) (𝑖 = 1, 2, . . . , ⌈ 𝑀𝑁𝑑 βŒ‰)-th rows are real numbers, and all other rows are zero. βŒˆΒ·βŒ‰ denotes the ceiling function. Here, in π‘š-th rows of 𝐴𝑇 𝐴, all columns except the ((π‘š βˆ’ 1) mod 𝑁𝑑 ) + 1 + 𝑁𝑑 (𝑖 βˆ’ 1))-th column are zero, and thus the inner product of the first column of (𝐴𝑇 𝐴)βˆ’1 and all rows except the (1 + 𝑁𝑑 (𝑖 βˆ’ 1))-th row of 𝐴𝑇 𝐴 are zero. From the above discussion, the product of 𝐴𝑇 𝐴 and the first column of (𝐴𝑇 𝐴)βˆ’1 can be expressed as: βŽ› ⎞ βŽ› ⎞ 𝑀1,1 1 ⎜ 𝑀1+𝑁𝑑 ,1 ⎟ ⎜0⎟ ⎜ ⎟ ⎜ ⎟ 𝐿 ⎜ 𝑀1+2𝑁 ,1 ⎟ 𝐡𝑄 ⎟ = ⎜0⎟ (18) ⎜ ⎟ ⎜ 𝑑 .. ⎟ ⎜.⎟ ⎠ ⎝ .. ⎠ ⎜ ⎝ . 𝑀1+(π‘„βˆ’1)𝑁𝑑 ,1 0 𝑣1 𝑣(1+𝑁𝑑 ) Β·Β·Β· Β·Β·Β· Β·Β·Β· 𝑣(1+(π‘„βˆ’1)𝑁𝑑 ) βŽ› ⎞ ⎜ 𝑣(1+𝑁𝑑 ) ⎜ 𝑣1 𝑣(1+𝑁𝑑 ) Β·Β·Β· Β·Β·Β· 𝑣(1+(π‘„βˆ’2)𝑁𝑑 ) ⎟ ⎟ ⎜ 𝑣(1+2𝑁 ) 𝑑 𝑣(1+𝑁𝑑 ) 𝑣1 Β·Β·Β· Β·Β·Β· 𝑣(1+(π‘„βˆ’3)𝑁𝑑 ) ⎟ 𝐿 (19) ⎜ ⎟ 𝐡𝑄 = ⎜ .. .. .. .. .. ⎟ ⎜ ⎜ . . . . . ⎟ ⎟ ⎜ . .. .. .. .. .. ⎟ ⎝ . . . . ⎠ 𝑣(1+(π‘„βˆ’1)𝑁𝑑 ) Β·Β·Β· Β·Β·Β· 𝑣(1+2𝑁𝑑 ) 𝑣(1+𝑁𝑑 ) 𝑣1 where 𝑄 = βŒˆπ‘€ β€² /𝑁𝑑 βŒ‰. When 𝑄 = 1, 𝑀1,1 = 1/𝐿 according to Eq. (18). 𝑄 becomes one only when the transmission period 𝑁𝑑 is greater than 𝑀 β€² , which corresponds to the averaging method. When 𝑄 > 1, 𝑀1,1 is the (1,1) component of (𝐡𝑄 𝐿 )βˆ’1 . Therefore, we investigate (𝐡 𝐿 )βˆ’1 . In 𝑄 the following, we impose the condition 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1 for 𝐿 and 𝑄. Then each element of 𝐡𝑄 𝐿 can be expressed as follows: 𝐿 βˆ’ 1 Β·Β·Β· Β·Β·Β· πΏβˆ’π‘„+1 βŽ› ⎞ 𝐿 ⎜ .. .. ⎟ ⎜ πΏβˆ’1 𝐿 . . ⎟ ⎜ ⎟ 𝐡𝑄𝐿 =⎜ ⎜ .. .. .. .. .. ⎟ (20) ⎜ . . . . . ⎟ ⎟ ⎜ .. .. ⎟ ⎝ . . 𝐿 πΏβˆ’1 ⎠ πΏβˆ’π‘„+1 Β·Β·Β· Β·Β·Β· 𝐿 βˆ’ 1 𝐿 𝐿 is as follows. Theorem 1. When 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1, the inverse matrix of 𝐡𝑄 βŽ› ⎞ 𝛽1 βˆ’0.5 0 Β·Β·Β· Β·Β·Β· Β·Β·Β· 0 𝛽2 βŽœβˆ’0.5 1 βˆ’0.5 0 Β·Β·Β· Β·Β·Β· Β·Β·Β· 0 ⎟ ⎜ ⎟ ⎜ .. ⎟ ⎜ 0 ⎜ βˆ’0.5 1 βˆ’0.5 0 Β·Β·Β· Β·Β·Β· . ⎟ ⎟ ⎜ . .. .. .. .. .. .. ⎟ ⎜ . . . . . . ⎜ . . ⎟ ⎟ 𝐿 βˆ’1 (𝐡𝑄 ) =⎜ . .. ⎟ . (21) ⎜ . .. .. .. .. .. ⎜ . . . . . . . ⎟ ⎟ ⎜ . ⎟ ⎜ . Β·Β·Β· Β·Β·Β· βˆ’0.5 βˆ’0.5 ⎟ ⎜ . 0 1 0 ⎟ ⎜ ⎟ ⎝ 0 Β·Β·Β· Β·Β·Β· Β·Β·Β· 0 βˆ’0.5 1 βˆ’0.5⎠ 𝛽2 0 Β·Β·Β· Β·Β·Β· Β·Β·Β· 0 βˆ’0.5 𝛽1 2𝐿 βˆ’ 𝑄 + 2 𝛽1 = (22) 4𝐿 βˆ’ 2𝑄 + 2 1 𝛽2 = . (23) 4𝐿 βˆ’ 2𝑄 + 2 Proof. We show that Eq. (21) is 𝐡𝑄𝐿 (𝐡 𝐿 )βˆ’1 = 𝐼. For the first column of (𝐡 𝐿 )βˆ’1 , the inner 𝑄 𝑄 product of the first row of 𝐡𝑄 𝐿 is (2𝐿 βˆ’ 𝑄 + 2)𝐿 πΏβˆ’π‘„+1 βˆ’ 0.5(𝐿 βˆ’ 1) + = 1. (24) 4𝐿 βˆ’ 2𝑄 + 2 4𝐿 βˆ’ 2𝑄 + 2 The inner product of 𝐡𝑄 𝐿 with the 𝑖 (𝑖 = 2, 3, . . . , 𝑄)-th row is: (𝐿 βˆ’ 𝑖 + 1)(2𝐿 βˆ’ 𝑄 + 2) πΏβˆ’π‘„+𝑖 βˆ’ 0.5(𝐿 βˆ’ 𝑖 + 2) + = 0. (25) 4𝐿 βˆ’ 2𝑄 + 2 4𝐿 βˆ’ 2𝑄 + 2 For the 𝑖 (𝑖 = 2, 3, . . . , 𝑄 βˆ’ 1)-th column of (𝐡𝑄 𝐿 )βˆ’1 , the inner product of the 𝑖-th row of 𝐿 is 𝐡𝑄 βˆ’ 0.5(𝐿 βˆ’ 1) + 𝐿 βˆ’ 0.5(𝐿 βˆ’ 1) = 1. (26) For the inner product of the 𝑖 (𝑖 = 2, 3, . . . , 𝑄 βˆ’ 1)-th column of (𝐡𝑄 𝐿 )βˆ’1 and 𝑗 (𝑗 ΜΈ= 𝑖)-th row of 𝐡𝑄𝐿 , if 𝑖 > 𝑗, βˆ’0.5(𝐿 βˆ’ 1 + 𝑖 βˆ’ 𝑗) + (𝐿 βˆ’ 2 + 𝑖 βˆ’ 𝑗) βˆ’ 0.5(𝐿 βˆ’ 3 + 𝑖 βˆ’ 𝑗) = 0. (27) If 𝑖 < 𝑗, then βˆ’0.5(𝐿 βˆ’ 1 + 𝑖 βˆ’ 𝑗) + (𝐿 + 𝑖 βˆ’ 𝑗) βˆ’ 0.5(𝐿 + 1 + 𝑖 βˆ’ 𝑗) = 0. (28) We considered the 𝑄-th column of (𝐡𝑄 𝐿 )βˆ’1 . Each of 𝐡 𝐿 and (𝐡 𝐿 )βˆ’1 is equal to the matrix 𝑄 𝑄 that inverts the row and column orders. Therefore, by the same reasoning as in the first column, the inner product of the 𝑄-th column of (𝐡𝑄 𝐿 )βˆ’1 with the 𝑄-th row of 𝐡 𝐿 is 1, and the other 𝑄 rows of 𝐡𝑄 𝐿 are 0, which concludes the proof. Figure 1: Examples of noise suppression performance This indicates that the noise suppression performance 𝑀1,1 can be expressed in the form 𝛽1 (Eq. (22)). The above discussion is for the first column of (𝐴𝑇 𝐴)βˆ’1 , but holds for the second and subsequent columns, and the value of each element can be calculated using (𝐡𝑄𝐿 )βˆ’1 . Therefore, we can represent the noise suppression performance in the decision interval as 𝛽1 (Eq. (22)). 3.3.2. Relationship between 𝐿, 𝑄, and noise suppression performance Figure 1 shows the results of calculating the (1,1) components of (𝐡𝑄 𝐿 )βˆ’1 for 𝐿 and 𝑄. First, as mentioned in the previous section, 𝑄 = 1 corresponds to the averaging method, and the noise suppression performance is 1/𝐿. For 𝑄, the noise suppression performance decreased as 𝑄 increased. For 𝐿, the noise suppression performance increases as 𝐿 increases but appears to converge to 0.5, when 𝑄 is not one. In the following, we show that the above properties are maintained under the condition 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1. We will discuss the case 𝐿 < 𝑄 βˆ’ 1 in future work. First, we investigate the noise suppression performance when 𝑄 increases. If 𝐿 β‰₯ π‘„βˆ’1 β‰₯ 1 is satisfied, the noise suppression performance can be expressed by Eq. (22), and thus the difference between the noise suppression performance when 𝑄′ and 𝑄′ +1 (𝐿 β‰₯ 𝑄′ +1 > 𝑄′ β‰₯ π‘„βˆ’1 β‰₯ 1) is 2𝐿 βˆ’ (𝑄′ + 1) + 2 2𝐿 βˆ’ 𝑄′ + 2 βˆ’ (29) 4𝐿 βˆ’ 2(𝑄′ + 1) + 2 4𝐿 βˆ’ 2𝑄′ + 2 This expression can be simplified as 1 . (30) 2(2𝐿 βˆ’ 𝑄′ )(2𝐿 βˆ’ 𝑄′ + 1) The equation is positive because 𝐿 β‰₯ 𝑄′ β‰₯ 1. In other words, as 𝑄 increased, the noise suppression performance decreased. The value of 𝑄 increased as the transmission period 𝑁𝑑 decreased. Therefore, the noise suppression performance decreased as the transmission period was shortened. Next, we investigate the noise suppression performance as the number of transmissions 𝐿 increases. The difference between 𝛽1 when the number of transmissions 𝐿′ (𝐿′ β‰₯ 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1) and 𝐿′ + 1 are 2(𝐿′ + 1) βˆ’ 𝑄 + 2 2𝐿′ βˆ’ 𝑄 + 2 βˆ’1 β€² βˆ’ β€² = . (31) 4(𝐿 + 1) βˆ’ 2𝑄 + 2 4𝐿 βˆ’ 2𝑄 + 2 (2𝐿 βˆ’ 𝑄 + 1)(2𝐿′ βˆ’ 𝑄 + 3) β€² The equation is negative because 𝐿′ β‰₯ 𝑄 β‰₯ 1. Therefore, the noise suppression performance improves when the number of transmissions is increased. In addition, from Eq. (22), it can be observed that 𝛽1 converges to 0.5 as 𝐿 is increased. To summarize the above, it can be said that under the condition that 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1, the noise suppression performance can be improved by increasing the number of transmissions 𝐿, but there is an lower limit of 0.5. The lower limit of 𝛽1 0.5 is the same as that of the noise suppression performance of the averaging method when the transmission period is set to more than 𝑀 and the number of transmissions is set to two, as described in Section 2.3.2. In other words, when the proposed method was used with a transmission period of 𝑀 or less, its noise suppression performance was less than or equal to that of the averaging method under the aforementioned conditions. In contrast, the proposed method has the advantage of reducing the time required for the signal transmission. Specific numerical examples and performance comparisons are presented in Section 4.3. 3.4. Relationship between the starting point of the cutout and the error caused by the projection In this section, we consider the case in which 𝑦 does not contain a part of 𝑔. There are two cases in which the cutout is made either too early or too late, and we discuss them sequentially. If the starting point of the cutout is too early, and 𝑦 does not contain the last 𝑖 components of 𝑔, 𝑦 is represented as follows: 𝑂𝑖′ (οΈ‚ )οΈ‚ 𝑦= π‘₯+𝑛 (32) 𝐴1,π‘…βˆ’π‘– 𝐴1,π‘…βˆ’π‘– is the first to (𝑅 βˆ’ 𝑖)-th rows of 𝐴. Let π‘Žπ‘– be the 𝑖-th column of matrix 𝐴, then, the matrix on the right-hand side of Eq. (32) can be expressed as (33) (οΈ€ )οΈ€ π‘Ž1+𝑖 π‘Ž2+𝑖 Β· Β· Β· π‘Žπ‘€ β€² π‘Žπ‘€ β€² +1 Β· Β· Β· π‘Žπ‘€ β€² +𝑖 Here, π‘Žπ‘— (𝑗 > 𝑀 β€² ) is a vector consisting of a 𝑖-by-1 zero vector, and the first to (𝑀 β€² βˆ’ 𝑖)-th rows of π‘Žπ‘—βˆ’π‘– . The product of this and 𝐴 is: π‘Žπ‘‡1 βŽ› ⎞ )οΈ‚ ⎜ π‘Žπ‘‡ ⎟ 𝑂𝑖′ (οΈ‚ ⎜ 2 ⎟ (οΈ€ 𝑇 (34) )οΈ€ 𝐴 = ⎜ . ⎟ π‘Ž1+𝑖 π‘Ž2+𝑖 Β· Β· Β· π‘Žπ‘€ β€² +𝑖 . 𝐴1,𝑅′ βˆ’π‘– ⎝ .. ⎠ π‘Žπ‘‡π‘€ β€² This indicates that the first to (𝑀 β€² βˆ’ 𝑖)-th columns are identical to the (𝑖 βˆ’ 1)-th to 𝑀 β€² -th columns of 𝐴𝑇 𝐴. Thus, 𝑦 β€² can be represented as: 𝑂𝑖′ (οΈ‚ )οΈ‚ β€² 𝑇 βˆ’1 𝑇 𝑦 =𝐴(𝐴 𝐴) 𝐴 π‘₯ + 𝑃𝐴 𝑛 𝐴1,π‘…βˆ’π‘– (οΈ‚ β€²β€² )οΈ‚ 𝑂 𝑖 𝐾1 =𝐴 π‘₯ + 𝑃𝐴 𝑛. (35) 𝐼𝑀 β€² βˆ’π‘– 𝐾2 𝑂𝑖′′ is the 𝑖-by-(𝑀 β€² βˆ’π‘–) zero matrix. 𝐾1 and 𝐾2 are 𝑖-by-𝑖 and (𝑀 β€² βˆ’π‘–)-by-𝑖 matrices, respectively. From Eq. (34), 𝐾1 and 𝐾2 are not zero matrices. Therefore, it can be observed that 𝑦 β€² contains an error that is the product of the last 𝑖 components of π‘₯ and 𝐾1 , 𝐾2 . Note that this error does not occur when 𝑧2 β‰₯ 𝑖 because the last 𝑧2 components of π‘₯ are zero (Eq. (11)). If the starting point of the cutout is too late, and 𝑦 does not contain the first 𝑖 components of 𝑔, then 𝑦 is represented as follows: (οΈ‚ )οΈ‚ 𝐴𝑖+1,𝑅 𝑦= π‘₯ + 𝑛. (36) 𝑂𝑖′ The matrix on the right-hand side of this equation can be written as (37) (οΈ€ )οΈ€ π‘Ž1βˆ’π‘– π‘Ž2βˆ’π‘– Β· Β· Β· π‘Ž1 Β· Β· Β· π‘Žπ‘€ β€² βˆ’π‘– . Here, π‘Žπ‘— (𝑗 < 1) is a vector comprising the 𝑖-th to 𝑀 β€² -th rows of π‘Žπ‘—+𝑖 and 𝑖-by-1 zero vectors. The product of this and 𝐴 is: βŽ› 𝑇 ⎞ π‘Ž1 (οΈ‚ )οΈ‚ ⎜ π‘Žπ‘‡ ⎟ 𝐴𝑖+1,𝑅 ⎜ 2 ⎟ (οΈ€ 𝐴𝑇 (38) )οΈ€ β€² = ⎜ . ⎟ π‘Ž1βˆ’π‘– π‘Ž2βˆ’π‘– Β· π‘Žπ‘€ β€² βˆ’π‘– . 𝑂𝑖 ⎝ .. ⎠ π‘Žπ‘‡π‘€ β€² It can be observed that (𝑖 + 1)-th to 𝑀 β€² -th columns of this matrix are identical to first to (𝑀 β€² βˆ’ 𝑖)-th columns of 𝐴𝑇 𝐴. Thus, 𝑦 β€² can be represented as: (οΈ‚ )οΈ‚ β€² 𝑇 βˆ’1 𝑇 𝐴𝑖+1,𝑅 𝑦 =𝐴(𝐴 𝐴) 𝐴 π‘₯ + 𝑃𝐴 𝑛 𝑂𝑖′ (οΈ‚ )οΈ‚ 𝐾3 𝐼𝑀 β€² βˆ’π‘– =𝐴 π‘₯ + 𝑃𝐴 𝑛. (39) 𝐾4 𝑂′′ 𝑖 𝐾3 and 𝐾4 are (𝑀 β€² βˆ’ 𝑖)-by-𝑖 and 𝑖-by-𝑖 matrices, respectively. From Eq. (38), 𝐾3 and 𝐾4 are not zero matrices. Therefore, this indicates that 𝑦 β€² contains errors owing to the product of the first 𝑖 components of π‘₯ and 𝐾3 , 𝐾4 . 4. Numerical Experiment In this section, numerical experiments were conducted to verify the effectiveness of the proposed method. In Experiment 1, we compared the performances of the conventional and proposed methods when the number of delay waves was varied. Experiment 2 was a comparative evaluation of the required transmission time and noise suppression performance. Table 1 Parameter Sampling frequency 𝑓𝑠 = 1/𝑇𝑠 48 kHz Frequency of transmission signal 𝑓 18 kHz Length of signal 𝑇 0.001 s Signal-to-noise raio 20 dB (1) False alarm probability 𝑃𝐹 𝐴 10βˆ’10 (2) False alarm probability 𝑃𝐹 𝐴 10βˆ’6 Transmission period 𝑇𝑑 0.04 s 𝑀′ 0.13×𝑓𝑠 Trial count 1000 4.1. Experimental setting This section describes the experimental setup common to Experiments 1 and 2. The list of parameters is presented in Table 1. If the parameters are not specified, the values listed in Table 1 are used. The transmitted signal was a sine wave in Eq. (2). The characteristics of the speaker and microphone were ignored here to evaluate only the influence of the delay waves. The number of samples in the received signal π‘Ÿ was set to 0.54 Γ— 𝑓𝑠 = 25, 920. For direct waves, the received time is not necessarily a constant multiple of 𝑁𝑠 when sampled. To simulate this, a uniform random number in the range of 0.15 s to 0.15 + 𝑁𝑠 s was generated for each trial and was set as the received time. The amplitude of the direct wave was set as 1. For the simulated delayed waves, the amplitudes of each delayed wave were uniformly random numbers, ranging from 0 to 0.8. The received time of each delayed wave was set to a uniform random number in a range from the received time of the direct wave to 0.08 s. (1) The threshold 𝜏 (1) was given values calculated from the false alarm probabilities 𝑃𝐹 𝐴 and the noise variance. 𝑀 β€² is set to a value of 0.13 s with a margin with respect to 𝑀 based on the discussion in Section 3.1. 4.2. Experiment 1: Comparison with conventional methods In this section, we compare the performance of each method when the number of delay waves is varied. In the proposed method, the number of signal transmissions is set to four. The following three methods were evaluated using conventional methods. The first is the maximum peak method described in Section 2. Note that for transmission, the signal length 𝑠(𝑑) is set to 𝑇 = 0.004 s, and the number of transmissions is set to one. The reason for this was to make the energy identical to that of the signal transmitted using the proposed method. The second is the leading-edge method described in Section 2. Here, 𝜏 in Section 2.3.1 is (2) calculated from the false-alarm probabilities 𝑃𝐹 𝐴 and noise variance. (2) The third is the leading-edge method using the threshold πœπ‘– . Hereafter, this method will be referred to as the low-threshold leading-edge method. The difference between this method and (2) the proposed method is that detection is performed using the threshold πœπ‘– without the noise (a) Number of delay waves: 0 (b) Number of delay waves: 25 (c) Number of delay waves: 50 Figure 2: Results of Experiment 1 Figure 3: Number of false alarm Figure 4: Results of Experiment 2 suppression process described in Section 3. This method was evaluated to show that the number of false alarms increased when the threshold was lowered, without the noise suppression process of the proposed method. The evaluation results of each method are shown in Figures 2a, 2b, and 2c. In these figures, the horizontal axis shows the error of the estimated index and the vertical axis shows the cumulative distribution function (CDF). The number of false alarms for each condition and method is shown in Figure 3. False alarms were defined as errors greater than 15. Figure 2a shows that, in the absence of delay waves, the maximum peak method yields the best performance, as described in Section 2. However, in the maximum peak method, the false alarm rate increases as the number of delay waves increases (Figures 2b, 2c, and 3). For the methods based on leading-edge detection (leading-edge, low-threshold leading-edge, and proposed methods), the performance was not significantly degraded even when the number of delay waves increased, as shown in Figures 2a, 2b, 2c, and 3. Moreover, Figures 2a, 2b, and 2c show that the proposed method reduces the error around the true value compared to the leading-edge method. Comparing the low-threshold leading-edge method and the proposed method, although the errors around the true value are not significantly different, the number of false alarms increases in the low-threshold leading-edge method compared to the proposed method, as shown in Figure3. 4.3. Experiment 2: Required transmission time and noise suppression performance This section presents a comparative study of the relationship between the required transmission time and noise suppression performance. In this experiment, the number of delay waves is set to 25. For the proposed method, the number of transmissions was set to one, two, three, and four. Their required transmission times are 0.13, 0.17, 0.21, and 0.25 s, respectively. For the averaging method, because 𝑀 is unknown, the transmission period 𝑇𝑑 is set to 𝑀 β€² Γ— 𝑇𝑠 = 0.13 s, and the number of transmissions is set to two. In this case, the required transmission time was 0.26 s. As described in Section 2.3.2, under these conditions, the noise suppression performance of the averaging method was 0.5, which is identical to the performance limit of the proposed method shown in the 3.3.2 section. The results of the experiment are shown in Figure. 4. As aforementioned, for the proposed method, the required transmission time increases as the number of transmissions increases. However, the performance of the estimation of the received time improves as the number of transmissions increases. The averaging method had the best performance in estimating the received time, but the required transmission time was the largest. These results confirm that the proposed method can provide noise suppression with reduced transmission time if the required noise suppression performance is greater than 0.5. 5. Conclusion In this study, we proposed a noise suppression method that can be applied to multiple signal transmissions without waiting for the disappearance of reverberation. Numerical experiments were conducted to compare the effectiveness of the proposed method with that of conventional methods in estimating the received time under multiple reverberation conditions. The relationship between noise suppression performance and 𝐿 and 𝑄 is discussed under the conditions 𝑄 = 1 or 𝐿 β‰₯ 𝑄 βˆ’ 1 β‰₯ 1. The case in which 𝐿 < 𝑄 βˆ’ 1 will be discussed in our future work. Acknowledgments This work was supported by JSPS KAKENHI Grant Numbers 20K23313 and 21K17797. References [1] J. Xiao, Z. Zhou, Y. Yi, L. M. Ni, A survey on wireless indoor localization from the device perspectivew, ACM Comput. Surv. 49 (2016). [2] M. A. Richards, J. A. Scheer, W. A. Holm (Eds.), Principles of Modern Radar: Basic principles, 1st ed., Scitech publishing, 2010. [3] M. O. Khyam, S. S. Ge, X. Li, M. Pickering, Orthogonal chirp-based ultrasonic positioning, Sensors 17 (2017). [4] M. O. Khyam, S. S. Ge, X. Li, M. R. Pickering, Pseudo-orthogonal chirp-based multiple ultrasonic transducer positioning, IEEE Sensors Journal 17 (2017) 3832–3843. [5] M. Hazas, A. Hopper, Broadband ultrasonic location systems for improved indoor posi- tioning, IEEE Trans. on Mobile Computing 5 (2006) 536–547. [6] F. HΓΆflinger, R. Zhang, J. Hoppe, A. Bannoura, L. M. Reindl, J. Wendeberg, M. BΓΌhrer, C. Schindelhauer, Acoustic self-calibrating system for indoor smartphone tracking (AS- SIST), in: IPIN 2012, 2012, pp. 1–9. [7] P. Lazik, A. Rowe, Indoor pseudo-ranging of mobile devices using ultrasonic chirps, in: ACM Sensys 2012, 2012, pp. 99–112. [8] M. Nakamura, H. Kameda, Indoor localization using multiple stereo speakers for smart- phones, in: IPIN 2019, volume 2498 of CEUR Workshop Proceedings, 2019, pp. 235–242. [9] M. Nakamura, H. Hashizume, M. Sugimoto, Simultaneous localization and communication method using short-time and narrow-band dual-carrier acoustic signals, IEEE Sensors Journal 22 (2022) 5163–5172. [10] T. Suzaki, M. Nakamura, H. Murakami, H. Watanabe, H. Hashizume, M. Sugimoto, Ex- panding the positioning area for acoustic localization using cots mobile devices, in: MobiQuitous2021, 2021, pp. 422–437. Appendix We show that 𝐴𝑇 𝐴 is full rank. Because 𝐴𝑇 𝐴 is a Gram matrix, it is sufficient to show that all the column vectors of 𝐴 are linearly independent. First, for the 𝑖 (1 ≀ 𝑖 ≀ 𝑁𝑑 )-th column vector, the 𝑖-th row is one, but all 𝑖-th rows of column vectors except 𝑖-th column vector are zero, and thus they are linearly independent. For the 𝑖 (𝑖 > 𝑁𝑑 )-th column vector, let Ξ› be the index set of the column vector whose 𝑖-th row is one, except for the 𝑖-th column as follows: {οΈ‚ βŒˆοΈ‚ βŒ‰οΈ‚ }οΈ‚ 𝑖 ((𝑖 βˆ’ 1) mod 𝑁𝑑 ) + 1 + 𝑁𝑑 (π‘š βˆ’ 1) | π‘š = 1, 2, . . . , βˆ’1 (40) 𝑁𝑑 In the following, let πœ†1 , πœ†2 , πœ†3 ∈ Ξ›: For πœ†1 -th rows, the value is zero in 𝑖-th columns but one only in πœ†2 (πœ†1 β‰₯ πœ†2 )-th columns. Therefore, if the 𝑖-th column vector is linearly dependent, then the linear sum of βŒˆπ‘–/𝑁𝑑 βŒ‰ βˆ’ 1 column vectors must satisfy the βŒˆπ‘–/𝑁𝑑 βŒ‰ equality constraint, where πœ†3 -th rows are zero and 𝑖-th row is one. Because this is impossible, it can be said that all 𝑖-th column vectors are linearly independent. Thus, (𝐴𝑇 𝐴) is shown to have full rank.