Robust Room Impulse Response Measurement Using Perfect Periodic Sequences for Wiener Nonlinear Filters †

: The paper discusses a measurement approach for the room impulse response (RIR), which is insensitive to the nonlinearities that affect the measurement instruments. The approach employs as measurement signals the perfect periodic sequences for Wiener nonlinear (WN) ﬁlters. Perfect periodic sequences (PPSs) are periodic sequences that guarantee the perfect orthogonality of a ﬁlter basis functions over a period. The PPSs for WN ﬁlters are appealing for RIR measurement, since their sample distribution is almost Gaussian and provides a low excitation to the highest amplitudes. RIR measurement using PPSs for WN ﬁlters is studied and its advantages and limitations are discussed. The derivation of PPSs for WN ﬁlters suitable for RIR measurement is detailed. Limitations in the identiﬁcation given by the underestimation of RIR memory, order of nonlinearity, and effect of measurement noise are analysed and estimated. Finally, experimental results, which involve both simulations using signals affected by real nonlinear devices and real RIR measurements in the presence of nonlinearities, compare the proposed approach with the ones that are based on PPSs for Legendre nonlinear ﬁlter, maximal length sequences, and exponential sweeps.


Introduction
Measuring the room impulse response (RIR) is a basic operation for acoustics and audio signal processing.
A plethora of approaches have been proposed in the literature for RIR measurement. The early approaches considered the use of impulsive signals [8,9], which directly estimate the impulse response of a linear system. Periodic pulses [10,11] were often used to contrast the effect of noise. Time-stretched pulses [12,13], composed by an expanded pulse excitation, have been proposed to overcome limitations in the amplitude of the pulses, which, if electronically generated, could damage the loudspeaker or could activate the protection circuit of the amplifier. Periodic sequences have been largely used in RIR estimation. An approach that is still popular today is maximal length sequence (MLS) technique of [14], where the excitation signal is a MLS, i.e., a binary pseudo-noise periodic sequence having period 2 n − 1 with n ∈ N [15]. The MLSs have an almost perfect autocorrelation function and, thus, they allow for estimating of the impulse response of a linear system with the cross-correlation method, i.e., computing the cross-correlation between output and input sequences. In reality, the estimation, also in noise-less conditions, is affected by a small error, which becomes negligible for large n. The perfect periodic sequences (PPSs) for linear systems [16][17][18] have a perfect autocorrelation function, a train of unit pulses, and thus, in noise-less conditions, allow the perfect estimation of the impulse response of a linear system with the cross-correlation method. Sweep signals have been largely employed for RIR: linear sweeps, exponential sweeps, and perfect sweeps. Linear sweeps are used in the time delay spectrometry technique [19,20]. The exponential sweeps, originally proposed in [21,22], are today the most popular technique for RIR measurement due to the robustness towards nonlinearities, as discussed later. Perfect sweeps are PPSs having the form of a periodic sweep and they have a perfect autocorrelation function [23]. Stepped sines, where the excitation signal is composed by pure tones in steps of increasing frequency, have also been used for RIR measurement [24].
A problem affecting many of these classical approaches is the sensitivity to the nonlinear components that are present in the measurement systems. A high volume of the measurement signal is often used to contrast the effect of noise. At this high volume, below the pain threshold, the acoustic path can still be considered as a linear system. Nevertheless, nonlinear effects can appear in the power amplifier or in the loudspeaker used for the measurement. The nonlinear effects are often the cause of artifacts in the measured signal. For example, the artifacts are particular evident in the MLS technique, where they assume the form of spikes in the measured signal [25]. The spikes are caused by the fact that the product of a MLS, b(n), with a delayed version of the same MLS, b(n − i), originates another delayed version of the same MLS, b(n − j) = b(n)b(n − i) [26]. In order to overcome these problems, RIR measurement techniques robust towards nonlinearities have been researched.
Nowadays, one of the most popular techniques is that based on exponential sweeps. The exponential sweep technique was independently developed by different researchers [21,22]. The approach has been improved in [27], where a synchronized exponential sweep technique is presented, remarking its importance for a proper analysis of higher harmonics. The popularity of the exponential sweep technique is due to its robustness towards the nonlinearities of the measurement system. It was shown in [21] that, if the measurement system can be modeled as a Hammerstein system, i.e., as a memoryless nonlinearity, followed by a linear system [28], the contribution of the nonlinear terms can be segregated at negative times in the measure RIR and it can be removed by windowing. In reality, it was later shown that the RIR measure performed according to [21,22,27] is still affected by the nonlinear kernels of the Hammerstein filter, but the measure can be corrected accounting for these nonlinear kernels [29] (In [29], but also in other papers [30][31][32], the exponential sweeps are used in order to derive an Hammerstein model to emulate nonlinear systems.). This correction is usually not performed by current measurement systems. Unfortunately, memoryless nonlinearities rarely occur and the measurement system is generally affected by nonlinearities with memory. It was shown in [33,34] that nonlinear distortions with memory also affect the exponential sweep technique and alter its RIR measurement.
In [35,36], a novel approach for RIR measurement robust towards nonlinearities with memory was proposed. In this approach, the entire measurement system, composed of the power amplifier, the loudspeaker, the acoustic path, and the microphone, is directly modeled as a Legendre nonlinear (LN) filter. LN filters are nonlinear filters with memory and they are linear combinations of polynomial basis functions. The basis functions are products of Legendre polynomials of the input samples and they are orthogonal for white uniform inputs [37]. LN filters admit PPSs, which, in the context of nonlinear filters, are periodic sequences that guarantee the perfect orthogonality of the basis functions over a sequence period. The PPSs for LN filters have samples with an approximate uniform distribution. Applying a PPS input, the coefficients of a LN filter can be estimated computing the cross-correlation between the system output and the basis functions, i.e., with the cross-correlation method. The first-order kernel of the LN filter can be estimated by simply computing the cross-correlation between the output and PSS input signals.
Another family of polynomial filters having orthogonal basis functions is the Wiener nonlinear (WN) filters, which derive from the truncation of the Wiener nonlinear series. Their basis functions are orthogonal for white Gaussian inputs. Additionally, WN filters admit PPSs [38] and, in this case, the sample distribution is approximately Gaussian. For this reason, the PPSs for WN filters appear more appealing for measuring RIRs: for the same input power, they less excite the highest amplitudes. Indeed, with a uniform distribution 42% of the samples have magnitude larger than the standard deviation, while with a Gaussian distribution the value reduces to 32% of the samples. In this paper, we extend the approach of [36] and consider the robust RIR measurement using PPSs for WN filters. This measurement was first proposed by the authors in the conference paper [39] and it is fully detailed in the current manuscript. We first review the theory of WN filters and their identification with PPS. The PPSs for WN filter were originally developed in [38], but those sequences are not suitable for RIR measurement: their period increases geometrically with the filter memory length and already for small memories it becomes prohibitively large. Thus, in this paper it is fully detailed how PPSs that are suitable for RIR identification with period proportionate to the filter memory length can be developed. Their derivation is nontrivial and is different from that of the PPSs for LN filters of [36]. A detailed analysis of RIR measurement under different non-ideal conditions, e.g., memory underestimation, order of nonlinearity underestimation, and noise effect, is also performed. The analysis of the noise effect provides novel results that were not discussed in [36]. The paper also provides novel experimental results that compare the proposed approach with competing approaches, i.e., MLSs, exponential sweeps, PPSs for LN filters. We consider both simulations that invove recorded signals affected by real nonlinear devices and convolved with a measured RIR, and real measurements performed in a room with professional equipment.
The main original contributions of this paper are the following: (i) full study of a novel RIR measurement approach based on PPSs for WN filters robust towards nonlinearities. (ii) Detailed derivation of PPSs for WN filters suitable for RIR measurement (the derivation is different from that for LN filters). (iii) Full analysis of detrimental effects that can affect the measurement, including an original study of the effect of Gaussian and non-Gaussian noises. (iv) Novel experimental results, including measurements performed in a room, when comparing the proposed approach with competing approaches (MLSs, exponential sweeps, and PPSs for LN filters).
It has to be pointed out that the proposed approach, like that of [35,36], is ineffective against the sub-harmonic distortions that may sometimes affect loudspeakers [40][41][42][43]. Single-input single-output polynomial models, like the WN, the LN, or the Volterra filters, simply cannot describe the sub-harmonic distortions [44].
The rest of the paper is organized, as follows. Section 2 reviews the WN filters and discuss their properties. Section 3 details the RIR measurement using PPSs for WN filters. Section 4 describes how PPSs that are suitable for RIR measurement can be developed and what are their limitations. Section 5 provides an analysis of RIR measurement under different non-ideal conditions. Section 6 presents the experimental results and illustrates the effectiveness of the proposed approach also in comparison with competing approaches. Section 7 provides the paper conclusions.

Wiener Nonlinear Filters
Like the Volterra filters derive from the double truncation with respect to order and memory of the Volterra series, the WN filters derive from the same double truncation of the Wiener series [28]. By applying the Stone-Weierstass theorem [45], it can be proved that WN filters can arbitrarily well approximate any discrete-time, time-invariant, finite-memory, causal, continuous, nonlinear system with input-output relationship where f is a continuous function from R N to R and N is the memory length of the filter [46]. WN filters are the linear combination of polynomial basis functions, which are the product of Hermite polynomials of variance σ 2 x of the delayed input samples. The Hermite polynomials can be generated with the recursive relation [47], where H j [ξ] the Hermite polynomial of degree j and the recursion starts with H 0 [ξ] = 1 and H 1 [ξ] = ξ. The first five Hermite polynomials of variance σ 2 x are reported in Table 1. For compactness, in what follows, the Hermite polynomials of order 0 and 1 will be indicated as 1 and ξ in place of H 0 [ξ] and H 1 [ξ], respectively, while the Hermite polynomials of order j = 2, 3, . . . will be indicated with H j [ξ]. The Hermite polynomials are orthogonal for ξ ∈ N (0, σ 2 x ), i.e., The WN basis functions can be formed following the same procedure of [37,48]. The Hermite polynomials are first written for ξ = x(n), x(n − 1), . . . , x(n − N + 1): , . . . and then the terms with different variables are multiplied in any possible manner, taking care to avoid repetitions. Let us define the diagonal number D as the maximum time difference between the input samples involved in the basis functions. Subsequently, Table 2 provides the WN basis functions of orders 1, 2, 3, memory length N, and diagonal number D. The basis function of order 0 is the constant 1 and for sake of simplicity it will be neglected in the following. The complete set of WN basis functions of memory N can be obtained by setting D = N − 1, but it has been shown in the literature that the basis functions with diagonal number D N − 1 are often sufficient to accurately model many real systems [49]. x(n), x(n − 1), . . . , x(n − N + 1) Order 3: A WN filter of order P, memory lenght N, and diagonal number D is the linear combination of all WN basis functions up to the order P, memory N, and diagonal number D. Neglecting the constant term, a WN filter of order 3, memory N, and diagonal number D is the linear combination of the basis functions in Table 2. The filter can also be implemented in the form of a filter bank with following relationŷ where b 2,i (n) and b 3,i,j (n) are the zero-lag basis functions of 2-nd and 3-rd order (the first elements of the rows of Table 1). Specifically, The naming convention of Volterra filters is used in the following. Thus, • h 1 (n) is the first order kernel, i.e., a length N sequence collecting the coefficients of the linear terms x(n − i). It should be noted that the nonlinear kernels are, in practice, band matrices with matrix bandwidth D [50].
We must remark that the first order kernel of WN filters does not coincide with the nonlinear filter impulse response, i.e., with the first order kernel of a Volterra filter, which iŝ whereŷ[Aδ(n)] is the filter response to a pulse of amplitude A. In fact, all odd order basis functions H 2k+1 [x(n − i)] with k ≥ 1 and i = 0, . . . , N − 1 include a linear term that contributes to the impulse response.
For the orthogonality property of Hermite polynomials in (3), the WN basis functions are orthogonal for a white Gaussian input with variance σ 2 x , i.e., when for any pair of basis functions b i (n) and b j (n) with i = j.
It was also shown in [38,51] that the WN filters admit PPSs that guarantee the same orthogonality of the WN basis functions on a finite period. The PPSs that were developed in [38,51] are periodic sequences whose samples have approximately a Gaussian distribution. With a PPS input, a WN filter can be identified with the cross-correlation method, computing the cross-correlation between the basis functions and the system output. Let < · > P denote time average over the period P. For a PPS input, for any pair of basis functions b i (n) and b j (n) with i = j, Accordingly, the coefficient of b i (n) can be obtained by computing the cross-correlation < b i (n)y(n) > P . For the RIR measurement, only the first order kernel of the WN filter has to be estimated and it is given by The cross-correlation can be performed in the time domain or, more efficiently, in the frequency domain, resorting to fast convolution.

Room Impulse Response Measurement
This section discusses the RIR measurement while using PPSs. The approach is the same of [35,36], with the difference that now WN filters are used in place of LN filters to model the measurement chain, and PPSs for WN filters are applied to measure the RIR. Figure 1 describes the measurement system, which is composed of a power amplifier, a loudspeaker, the room acoustic path, and a microphone. At the sound pressure levels (SPLs) that are normally used in the measurements, the room acoustic path is modelled as a linear system (harmonic distortions appears in the acoustic path only for SPLs greater than 140 dB SPL [52]). Given the low levels of the acquired signals, the microphone nonlinearities can also be neglected. The main sources of nonlinearities in the measurement system are to be found in the power amplifier and the loudspeaker, which are often over-driven during measurements to contrast the effect of noise. The measurement objective is to estimate the RIR h R (n), which is assumed to have length M.
The power amplifier and the loudspeaker are first considered as linear systems with impulse response h 1 (n) having memory length N, with h 1 (n) that can also include the microphone impulse response. The entire measurement chain is a linear system with impulse response having memory length In these ideal conditions, any measurement system can obtain a good estimate of h T (n), which is generally considered the measured room impulse response. According to (7), the measure of the true RIR h R (n) is, in reality, affected by the convolution with h 1 (n), the impulse response of the power amplifier and loudspeaker. h 1 (n) can be measured in an anechoic chamber and, according to some researchers, it can be used to improve the estimate of h R (n) by equalizing the equipment. For this purpose, in [53], the use of the Kirkeby algorithm was suggested. The algorithm is implemented, as follows: where FFT[·] and IFFT[·] are direct and inverse fast Fourier transform (FFT) operators, respectively, the division is performed at the single frequencies, and (ω) is a frequency-dependent regularization parameter. It has to be pointed out that the approach in (8) can obtain the RIR h R (n) only in the case of a perfectly omnidirectional loudspeaker. Because this is never the case, it is usually considered to be satisfactory the knowledge that h 1 (n) affects the measurement mildly and in a known manner. The power amplifier and the loudspeaker are next considered the source of nonlinear effects. The input-output relationship of the power amplifier and loudspeaker system is assumed to be a WN filter of order K, memory length N, and diagonal number D. For the simplicity of presentation the case of K = 3 is considered in the following and the input-output relationship of the power amplifier and loudspeaker system is given by (4). h 1 (n) is the first order kernel of the WN filter and it has length N. The input output relationship of the entire measurement system becomeŝ which is the input-output relationship of a WN filter with order K, diagonal number D, and memory Using as input a PPS for WN filters of order K, memory length P, and diagonal number D, h T (n) can be measured with (6). The orthogonality of the WN basis functions for a PPS input signal guarantees that the measurement of h T (n) is not affected by the nonlinear kernels h 2,i (n) and h 3,i,j (n) for all i, j. As argued for the linear case, h T (n) can be regarded as a measure of the true RIR h R (n). The measure is again affected by the convolution with h 1 (n), the first order kernel of the power amplifier and loudspeaker system. h 1 (n) can be measured in an anechoic chamber by using the same PPS and the same signal power used in measurement. Experiments showed that, in most cases, its frequency response is similar to that of the loudspeaker. In case the loudspeaker is omnidirectional, h 1 (n) could be used to estimate the true RIR h T (n) by equalizing the equipment with the Kirkeby algorithm in (8). For the properties of PPSs, the resulting RIR will not be affected by the power amplifier and loudspeaker nonlinearities in any way.

PPSs for RIR Measurement
This section explains how PPSs that are suitable for RIR identification with period proportional with the memory length can be developed.
We are interested in developing a PPS x p (n) of period P for a WN filter of order K, memory N T , diagonal number D, and Gaussian input variance σ 2 x . In order to allow reproduction with digital to analog converters, the PPS should be bounded by 1, i.e., |x p (n)| < 1 ∀n.
PPSs for WN filters were developed in [38,51]. The approach that is used for their derivation is very different from that of LN filters [35,36]. In LN filters, PPSs are derived by imposing the orthogonality of all basis functions over a period of the PPS. In contrast, with WN filters, the approach imposes that all joint moments estimated over a period involved in the filter estimation (i.e., up to order 2K) assume the ideal values of a Gaussian noise N (0, σ 2 x ). Specifically, for a full WN filter of order K, memory N, the following system of nonlinear equations is imposed for all r 0 , . . . , r N T −1 ∈ N with r 0 > 0 (for the periodicity of the sequence), r 0 + . . . + r N T −1 ≤ 2K. In (10) µ r is the r-th moment of the Gaussian process N (0, σ 2 x ), and q!! the double factorial, q!! = q · (q − 2) · (q − 4) · . . . · 1.
The system in (10) is a nonlinear equation system in the variables x p (n), for n = 0, . . . , L − 1. For sufficiently large P, the nonlinear equation system is underdetermined and it may have infinite solutions in the variables x p (n). Applying the Newton-Raphson method, a solution to the system in (10) has always been found. The Newton-Raphson method has been implemented, as in [54] (ch. 9.7). The iterations start from a random distribution of x p (n) in N (0, σ 2 x ). To obtain a sequence in [−1, +1], the samples x p (n) are reflected in that range every time they exceed it. The Newton-Raphson method converged in all simulations, provided σ 2 x was sufficiently small. This property can be understood when considering that the sample distribution is similar to a Gaussian and the probability of having samples outside [−1, +1] should be sufficiently small to guarantee the convergence. Solutions to (10) were found for L = 3Q ÷ 4Q, with Q the number of nonlinear equations, and σ 2 x ≤ 1/10 [51]. The PPSs that were developed with the system in (10) are not suitable for RIR measurement. The period P increases geometrically with the memory N T of the filter, and already for small N T (e.g., of 30-35 samples) it can be prohibitively large. The RIR measurement requires memory lengths of thousands of samples and the PPSs of [38,51] cannot be used. In reality, in (10) P depends geometrically on N T , because the PPS has to allow the estimation of all kernels of WN filter. In RIR measurement, only the first order kernel needs to be computed, as shown in the previous Section. Thus, we can relax our constraints and impose that only those joint moments involved in the measurement of the first order kernel assume the ideal values of a Gaussian noise N (0, σ 2 x ). The resulting nonlinear equation for all r 0 , . . . , r N T −1 ∈ N with r 0 > 0 (for the periodicity of the sequence), r 0 + . . . + r N T −1 ≤ K + 1, and with x r 0 The number of nonlinear equations in the reduced system (12), increases exponentially with the order K, geometrically with the diagonal number D, but only linearly with the memory length N T . The Newton-Raphson method has still computational and memory requirements that can be prohibitively large, since they increase with the cube of number of equations, i.e., Q 3 . As discussed in [37,55,56], it is thus very useful to impose specific structures to the PPS in order to further reduce the number of equations and variables. The periodic sequence can be a priori structured to exploit the properties of: • Symmetry: for any N-tuple of samples a 1 , a 2 , . . . , a N , there is also the reversed one a N , a N−1 , . . . , a 1 . For every couple of symmetric joint moments (e.g., < x(n)x 2 (n − 2) > P and < x 2 (n)x(n − 2) > P ) only one has to be considered in (12). • Oddness: for any N-tuple a 1 , a 2 , . . . , a N , there is also the opposite one −a 1 , −a 2 , . . . , −a N . All odd joint moments are a priori zero. • Oddness-1: for any N-tuple a 1 , a 2 , . . . , a N , there is also the one obtained by alternatively negating the samples a 1 , −a 2 , a 3 , . . . , −a N . All odd-1 joint moments are a priori zero. • . . . , • Oddness-2 R : for any N-tuple there is also the one obtained by alternatively negating blocks of 2 R samples. All odd-2 R joint moments are a priori zero (A basis function is defined odd-2 R if alternatively negating blocks of 2 R samples inverts the sign of its output.).
Each of these properties allows for almost halving the number of equations and variables in (12). Structured sequences, which can also be used for any PPS, can be found in [56]. The reduction in the number of equations comes at the cost of a longer period of the resulting PPS, since a structured PPSs is formed by a sequence of independent variables and several of its replica, flipped, or modified in sign. Nevertheless, the reduction of computational cost and memory requirements of the Newton-Raphson algorithm is remarkable and instrumental to efficiently solve the system in (12).
Many PPSs that are suitable for the RIR identification have been developed solving (12) for structured periodic sequences. The developed sequences can be freely downloaded from [57]. Table 3 summarizes the characteristics of some PPSs for WN filters suitable for RIR measurement. In the last column, the symbols S, O, O-R, mean that the sequence exploits symmetry, oddness, oddness-i with i ranging from 1 to R, respectively. All of the sequences of Table 3 have samples in the interval [−1, +1] and power 1/12.

Analysis of RIR Measurement Using PPSs
The sequences of Table 3 are specifically tailored for WN filters of a certain order K, memory length N T , and diagonal number D. This section provides an analysis of what happens when the measurement system has memory length, order, or diagonal number larger than those that are considered in the development of the PPS. Moreover, an analysis of the estimation error in different measurement noise conditions is also provided.

Memory Underestimation
We first consider the case where the measurement system can be modeled with a WN filter of order K and diagonal number D equal to those of the PPS, but has memory length N T larger than the memory length N T considered for the PPS. The measured first order kernel h 1 (n) is then affected by an aliasing error. The error mainly alters the first coefficients of h 1 (n). This can be easily understood considering N T = N T + 1. The coefficient h 1 (0) is affected by aliasing since < x(n)x(n − N T ) > P = 0. In contrast, the other coefficients are correctly measured since the conditions in (12) means < x(n − i)x(n − N T ) > P = 0, for all i = 1, . . . , N T .

Order or Diagonal Number Underestimation
We next consider the case where the measurement system has memory length equal to that of the PPS, but order K or diagonal number D greater than the order K or the diagonal number D considered for the PPS. The measurement is again affected by an aliasing error, which now involves all measured coefficients. Nevertheless, this aliasing error is often negligible when K and D are sufficiently large and if the PPS exploits some structure. In fact, structured PPSs have many basis functions with order and diagonal number greater than K and D, respectively, which have a priori zero cross-correlation with the linear basis functions x(n − i).
Even though this error is deterministic, working with a large period P and a large number of neglected basis functions, for the law of large numbers, the aliasing error can be considered to be stochastic and Gaussian distributed. Thus, the effect of the error can be considered similar to a measurement noise.

Measurement Noise Effect
Let us now assume the measurement system has memory, order, and diagonal number lower or equal to those considered in developing the PPS, and let us consider the effect of an additive measurement noise ν(n). The measured signal iŝ where h 1 (n) is the true first order kernel that we want to measure. We separately analyze the case of a Gaussian and of a non-Gaussian noise.

Gaussian Noise
We first consider the case ν(n) is a colored Gaussian noise, where ν(n) is a zero-mean, variance σ 2 ν , white Gaussian noise, and h ν (n) is the impulse response of the forming filter with memory N ν . Without loss of generality it is assumed Note that, when h ν (n) = δ(n), ν(n) is a zero-mean, variance σ 2 ν , white Gaussian noise. The first order kernel h 1 (n) obtained applying (6) (with y(n) =m(n)), for the properties of PPSs is Since ν(n) is zero mean, E[h 1 (n)] = h 1 (n) and the estimation is unbiased. We can estimate the Mean Square Deviation (MSD) of the coefficients that are defined by For the PPSs, from (6) it is and expanding the numerator it can be proved that When the measurement system is affected by a white Gaussian noise, i.e., when h ν (n) = δ(n), the MSD reduces to which highlights that the MSD improves proportionally to the period of the sequence, i.e., to the number of samples used for the measurement.

Non-Gaussian Noise
It is also interesting to note how the measurement behaves when the noise is zero mean, but non-Gaussian. In this case, Because ν(n) is zero mean, the estimation is still unbiased, since E[h 1 (n)] = h 1 (n). According to (6) and (13), for the central limit theorem, the coefficient error will have a Gaussian distribution with which depends on the autocorrelation function of the noise and on the PPS.

Experimental Results
The section provides experimental results that illustrate the effectiveness and robustness of the proposed approach in comparison with competing approaches. Specifically, the results obtained with the proposed approach are compared with those obtained with the PPSs for LN filters [35,36], the MLS [14], and the exponential sweeps technique [21,22].

First Experiment
An emulated environment is first considered in order to perform an exact estimation of the measured impulse responses. The measured system is composed of a nonlinearity, emulating the distortion of an amplifier or a loudspeaker, cascaded with a previously measured RIR. A real device, a Behringer MIC100 vacuum tube preamplifier, was used to introduce nonlinear effects in the measurement chain. Working at a 48 kHz sampling frequency, different PPSs for WN and LN filters, MLSs, and exponential sweeps were applied to the pre-amplifier and the corresponding output has been recorded. In order to exploit the maximum signal to noise ratio of the DAC that drives the amplifier, the different input sequences had the same peak amplitude. Acting on a potentiometer of the preamplifier, it was possible to control the amount of nonlinear distortion introduced and ten different settings were considered. Figure 2 shows the second, third, and total harmonic distortion in percent on a 1 kHz tone having the same peak amplitude of the input sequences. The harmonic distortion is defined as the percentage ratio between the power of the selected harmonic and that of the fundamental frequency. Clearly, many harmonic distortions of Figure 2 are much larger than those of a regular RIR measurement system, but they have been selected so large to stress the robustness of the proposed approach and the differences between the compared methods. The recorded output of the preamplifier was then convolved with a known RIR, measured in an auditorium (as shown in Figure 3), and a white Gaussian noise was added to the output. The added Gaussian noise has the same power for all methods and the signal to noise ratio is of 40 dB for the MLSs, 37 dB for the exponential sweep, 35.2 dB for the PPS for LN filters, and 29.2 dB for the PPS for WN filters. The auditorium RIR has a memory length that is lower than 8000 samples. Thus, the PPSs for WN filters 1-5 of Table 3 were considered in the simulations, with sequence orders ranging from 17 to 21. The sequence order is defined as the base-2 logarithm of the period or length of the sequence. Similar PPSs have been considered also for the LN filter. The RIR was identified also with MLSs and exponential sweeps with orders that range from 15 to 21. The sequence orders 15 and 16 are not available for the PPSs, but they were also included because they are sufficient to identify the RIR and they allow us to better assess the effect of nonlinearities on these identification techniques. All of the measurements were performed on one period of the periodic sequences or on a single exponential sweep. From the measurement, an impulse response of 8192 samples was extracted and all other samples, including the harmonic distortion components of the exponential sweeps, were windowed out. The simulation results are compared in terms of log-spectral distance (LSD), which allows for us to evaluate the differences between two audio signals [3,58]. The LSD is defined as the root mean square value in a band B of the difference between the spectra of the two signals expressed in dB.
When considering B = [k 1 , with F S the sampling frequency and N DFT the number of the samples of the discrete Fourier Transform (DFT), the LSD is defined as where |H R (k)| is the actual system magnitude response in the DFT domain and |Ĥ R (k)| is the estimated system magnitude response, normalized in order to have the same mean of |H R (k)| in the band B.
In the experiments, the LSD was computed in the band B = [100, 18,000] Hz, falling strictly inside the measurement system passband. Figure 4 shows the LSD between the identified and actual impulse response of the measurement chain for the first experiment. Panel (a) of Figure 4 shows the LSD that was obtained using the PPSs for WN filters, panel (b) using the PPSs for LN filters, panel (c) using the MLSs, and finally panel (d) using the exponential sweeps. The actual impulse response of the measurement chain has been estimated in noise absence at the lowest distortion setting (setting 1) while using an exponential sweep of order 21. Figure 4 shows that, for large distortions, the performance of the MLSs is largely affected. Better results are obtained with the exponential sweeps, especially for order 17 or greater, with a smooth deterioration of the LSD increasing the distortion. The PPSs for LN filters also provide very good results, with a deterioration of the LSD performance that increases with the distortion, but at a lower rate than with the exponential sweeps. In particular, the larger the diagonal number D of the PPS, i.e., the larger the order log 2 (P) and the protection against the nonlinear terms, the better the performance at highest distortion levels. In contrast, in the considered conditions, the PPSs for WN filters provide an almost constant LSD for all settings showing great robustness towards nonlinearities. This improvement is due to the fact that the WN filters with their almost Gaussian sample distribution excite the nonlinearities less. At low distortion levels, the LSD obtained with PPS for WN filters is slightly worse than for the LN filters, because the output signal to noise ratio is 6 dB lower. If the same signal to noise ratio is imposed for all sequences, the LSD of PPSs for WN filters is at least as good as that of the other methods. Next, the effect of the equalization of the equipment, in particular of the distortion nonlinearity, is studied. The response of the MIC100 preamplifier is identified and it is used to equalize the measured RIR with the Kirkeby algorithm in (8). Figure 5 shows the LSD between the auditorium RIR of Figure 3, and the RIR estimated with the equalization of the equipment. Different panels report, as in the previous Figure, the LSD obtained, respectively, with PPSs for WN filter, PPSs for LN filter, MLSs, and finally with exponential sweeps. For large distortion levels, the estimations with MLSs are still affected by the nonlinearities. The LSDs that were obtained with the exponential sweeps improves at the highest settings. The improvement can be explained with the fact that the nonlinearity affects, in a similar manner, the measurement of the impulse response of the equipment and of the room, and the equalization of the equipment compensates in part the effect of the nonlinear distortion. Additionally, the LSD that was obtained with the PPS for LN filters improves and, for large orders of the periodic sequence, the measurement appears unaffected by the nonlinearity. The PPSs for WN filters provide an almost constant LSD for all settings, confirming the good robustness against nonlinear effects. In the considered conditions, in the case of PPSs for WN filters, the equalization of the equipment compensates the frequency response of the preamplifier, but does not take to any further improvement in terms of reduction of nonlinear effects. From panel (a) and (b) of Figures 4 and 5, we can also appreciate the improvement in the log-spectral distance obtained by increasing the diagonal number and the period of the PPSs. As a rule of thumb, the researcher or technician performing the RIR measurement should first choose the memory length N T according to the considered acoustic scenario, and then he should choose the largest order K and diagonal number D that still provide a period of the PPS acceptable for performing the measurement.

Second Experiment
The second experiment considers the measurement of a real impulse response in the presence of nonlinearities. The measurement was performed in a living room with dimensions of 3.9 m × 4 m × 3.2 m employing a loudspeaker RCF Ayra Five, a microphone Behringer ECM8000, and a Focusrite Scarlett 2i2 (2nd gen) audio interface. The nonlinear distortion was generated with the guitar pedal Behringer DM100 (set in D+ mode). The pedal was used in order to introduce a measurable strong nonlinear distortion in the measurement chain to stress the differences between the different measurement methods. The room impulse response was measured with a 48 kHz sampling frequency using the same PPSs for WN and LN filters, MLSs, and exponential sweeps of the previous experiment. The measurement was performed with three different settings. In the first setting, the test sequences are directly applied to the loudspeaker without passing through the guitar pedal. In this condition, the measurement system introduces a very low nonlinear distortion and the measurement is used as reference for assessing the performance when the strong nonlinearity is introduced. The second setting measures the same room impulse response without altering the position of the loudspeaker or the microphone, but feeding the test signals through guitar pedal before applying it to the loudspeaker. The third setting measures the response of the guitar pedal with the same test signals and it allows for the equalization of the distortion pedal. To give an indication of the nonlinear distortions during the measurements, Table 4 provides for each setting the second, third, and total harmonic distortions on a sinusoidal signal with peak amplitude that is equal to that of the test signals.  Figure 6 provides the LSD between the reference RIR and the RIRs that were measured in the presence of the pedal distortion, as a function of the sequence period. The reference RIR is the impulse response measured in the first setting with an exponential sweep of length 2 21 . The LSD has been estimated in the bands B = [100, 18,000] Hz strictly inside the passband. Panel (a) of Figure 6 reports the LSD values without the Kirkeby equalization of the pedal, whose effect is shown in panel (b). 15 Figure 6. Log-spectral distance in the band [100, 18, 000] Hz versus the sequence order log 2 (P) for the different measurement sequences (a) without equalization of the equipment, (b) with equalization of the equipment.
In this experiment, without the equalization of the equipment, the measurements with PPSs always provide the best results. Similar results are obtained here both with PPSs for LN and WN filters. The MLSs and exponential sweeps only provide slightly worse results with a LSD increase of around 0.2 dB. In this experiment, the MLSs provide much better results than in the previous one: the MLSs appear much sensitive to the third order nonlinearities, but less sensitive to second order ones. With the equalization of the equipment, all of the methods provide very similar results that are much better than those without the equalization. When considering that the equalization of the equipment is rarely performed, the proposed RIR measurement system based on PPSs for WN filters appears to be a valid candidate for contrasting the effect of nonlinearities in the measurement chain.

Conclusions
The paper introduces a novel technique for RIR measurement that is based on PPSs for WN filters. The proposed technique is robust towards the nonlinearities that may affect the power amplifier or the loudspeaker of the measurement system. The approach is based on the estimate of the first order kernel of a WN filter modeling the acoustic path. The first order kernel is estimated with the cross-correlation method while using a PPS input signal. The identification of this first-order kernel has a computational cost similar to the RIR measurement with MLSs, exponential sweeps, or PPSs for LN filters (The identification costs around PN T operations, i.e., multiplications and additions, if the cross-correlation is computed in time domain. It costs (2 log 2 (P) + 1)P operations if the cross-correlation is computed in the DFT domain, assuming a FFT cost P log 2 (P) operations.). The proposed approach provides an improvement in comparison with that based on PPSs for LN filters, since the almost Gaussian distribution of the PPSs samples for WN filters less excites the nonlinearities of the measurement system. PPSs for WN filters that are suitable for RIR measurement have also been developed within the paper. Their properties have been analyzed in theory studying the effect of an underestimation of the order, diagonal number, of memory length of the WN filter modeling the acoustic path. The effect of the measurement noise has also been studied. The experimental results, when considering both an emulated scenario, involving signals that are affected by a real nonlinear device, and a real room impulse response measurement, highlight the effectiveness and the robustness towards nonlinearities of the proposed RIR measurement technique. Future research will be dedicated to the study and development of further periodic sequences suitable for RIR measurement, as the recently introduced orthogonal periodic sequences [59].