A frequency domain adaptive filtering method for identifying a low-rank acoustic system

Abstract

The application discloses a frequency domain adaptive filtering method for identifying a low-rank acoustic system, which extends the nearest Kronecker product (NKP) to the frequency domain, and establishes an NKP-based frequency domain recursive least square (NKP-FRLS) algorithm for identifying a time-varying acoustic system. The NKP is used to decompose an adaptive filter with a length of L into two groups of sub-filters with lengths of L1 and L2, and a signal model and a recursive least square cost function are established, the input signal spectrum matrix, the sub-filter error vector, the power spectrum matrix and the Kalman filter gain matrix are calculated in sequence, and on this basis, the sub-filter vector is calculated to obtain the coefficient vector of the modeling filter, so that the calculation amount of the adaptive filter is reduced, and a Gaussian noise robust adaptive algorithm is obtained. Experiments show that the application is superior to the traditional frequency domain recursive least square (FRLS) algorithm in convergence performance and calculation efficiency.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of acoustic system identification, and more particularly relates to a frequency-domain adaptive filtering method for identifying a low-rank acoustic system. BACKGROUND

[0002] Acoustic system identification plays an important role in many applications, such as acoustic echo cancellation, active noise control, sound source localization, and speech dereverberation. Since acoustic systems are usually time-varying, adaptive filtering methods are needed to identify such systems. In practical applications, the impulse response of an acoustic channel can be very long, from hundreds to thousands of sample intervals, which makes it difficult for adaptive algorithms to accurately identify such long-tailed time-varying acoustic channels.

[0003] To solve this challenging problem, adaptive filtering algorithms based on shortening the impulse response of an acoustic channel have been extensively studied in the past two decades. Among these algorithms, the recently proposed Nyman-Kaup-Peek (NKP) decomposition method has attracted wide attention from researchers. This method uses NKP decomposition to decompose the impulse response of an acoustic system and obtains a low-rank approximation. Based on this, the high-dimensional acoustic system identification problem is converted into a low-dimensional acoustic system identification problem, which greatly reduces the computational complexity of the algorithm. However, the performance of the algorithm is greatly affected by the parameter values.

[0004] Another technique to reduce computational complexity is the frequency-domain adaptive filtering method, such as the frequency-domain recursive least squares (FRLS) algorithm, which uses fast Fourier transform (FFT) to improve computational efficiency when identifying the impulse response of an acoustic channel. However, there is room for improvement in terms of computational efficiency and convergence performance. SUMMARY

[0005] The present application aims to overcome the shortcomings of the prior art and provide a frequency-domain adaptive filtering method for identifying a low-rank acoustic system to improve the computational efficiency and convergence performance of a frequency-domain adaptive filter.

[0006] To achieve the above-mentioned application purposes, the frequency-domain adaptive filtering method for identifying a low-rank acoustic system according to the present application comprises the following steps:

[0007] (1) Initialization

[0008] First, a first group of P sub-filters with a length of L1 is constructed and a second group of P sub-filters with a length of L2 is constructed m is the block time index, and then the initialization is performed:

[0009]

[0010]

[0011] where parameter η is selected in the range of greater than 0 and less than or equal to 1, P < min{L1, L2}, thus the first and second groups of sub-filter vectors of block time index m = 0 are:

[0012]

[0013]

[0014] Initialize two groups of power spectrum matrices of block time index m = 0

[0015]

[0016]

[0017] wherein and are unit matrices of sizes PL1xPL1and PL2xPL2, respectively, and parameters δ1, δ2 are selected according to specific implementation;

[0018] (2) Collect the sound signal to obtain a sample value x(n) and construct a sound signal vector x(m) corresponding to block time index m:

[0019] x(m) = [x(mL-L)x(mL-L+1)...x(mL+L-1)] T

[0020] wherein L = L1L2, and the value of the vector at a time less than or equal to 0 is 0;

[0021] (3) Starting from block time index m = 1, calculate the coefficient vector of the adaptive filter

[0022] 3.1) Calculate the spectrum matrix of the input signal

[0023]

[0024] wherein diag{·} represents expanding a vector into a diagonal matrix, F 2L×2L is a Fourier matrix of size 2Lx2L;

[0025] 3.2) Calculate the sub-filter vector

[0026] 3.2.1) Calculate the sub-filter spectrum matrix

[0027] ​

[0028] wherein denotes the Kronecker product, I L×L is an L x L identity matrix, is an L1x L1identity matrix, 0 L×L is an L x L zero matrix;

[0029] Thus, the combined sub-filter spectrum matrix

[0030]

[0031] 3.2.2), calculating the frequency-domain sub-filter error vector ∈ 1(m):

[0032]

[0033] wherein y (m) = F L×L y(m), y(m) = [y(mL) y(mL+1)... y(mL+L-1)] T is the observation signal vector corresponding to block time index m, y(n) is the observation signal sample; F L×L is an L x L Fourier matrix;

[0034] 3.2.3), calculating the power spectrum matrix

[0035]

[0036] wherein λ1is a forgetting factor, satisfying 0 < λ1< 1, the superscript H denotes the conjugate transpose, Re denotes the real part;

[0037] 3.2.4), calculating the Kalman filter gain matrix κ2(m):

[0038]

[0039] 3.2.5), calculating the sub-filter vector

[0040]

[0041] wherein μ1is an adaptive step size, satisfying 0 < μ1< 2,

[0042] 3.3), calculating the sub-filter vector

[0043] 3.3.1), calculating the sub-filter spectrum matrix

[0044]

[0045] wherein, is a unit matrix of size L2 x L2;

[0046] Thus, the combined sub-filter spectrum matrix

[0047]

[0048] 3.3.2), calculating the frequency-domain sub-filter error vector ∈ 2(m):

[0049]

[0050] 3.3.3), calculating the power spectrum matrix

[0051]

[0052] wherein λ2 is a forgetting factor, satisfying 0 < λ2 < 1;

[0053] 3.3.4), calculating the Kalman filter gain matrix κ1(m):

[0054]

[0055] 3.3.5), calculating the sub-filter vector

[0056]

[0057] wherein μ2 is an adaptive step, satisfying 0 < μ2 < 2;

[0058] 3.4), calculating the filter coefficient vector

[0059]

[0060] (4), obtaining the filter coefficient vector from the block time index m.

[0061] The invention aims to achieve:

[0062] The application provides a frequency domain adaptive filter method for identifying a low-rank acoustic system, which extends a nearest Kronecker product (NKP) to a frequency domain, and establishes a NKP-based frequency domain recursive least square (NKP-FRLS) algorithm for identifying a time-varying acoustic system and Thus, a signal model and a recursive least square cost function are established, and an input signal spectrum matrix, a sub-filter error vector, a power spectrum matrix and a Kalman filter gain matrix are calculated in sequence, and on this basis, a sub-filter vector is calculated and further a coefficient vector of a modeling filter is obtained Thus, the calculation amount of the adaptive filter is reduced, and a Gaussian noise robust adaptive algorithm is obtained. Experiments show that the application is superior to a traditional frequency domain recursive least square (FRLS) algorithm in terms of convergence performance and calculation efficiency. BRIEF DESCRIPTION OF DRAWINGS

[0063] Figure 1 is a specific embodiment flow chart of the frequency domain adaptive filter method for identifying a low-rank acoustic system of the application;

[0064] Figure 2 is a convergence comparison chart of FRLS and NKP-FRLS algorithms (the application) in identifying a time-varying acoustic system when a white Gaussian sequence is used as an input signal under the condition of a signal-to-noise ratio SNR=15dB;

[0065] Figure 3 is a convergence comparison chart of FRLS and NKP-FRLS algorithms in identifying a time-varying acoustic system when speech is used as an input signal under the condition of a signal-to-noise ratio SNR=15dB. DETAILED DESCRIPTION

[0066] The specific embodiments of the application are described below in conjunction with the accompanying drawings, so that those skilled in the art can better understand the application. It should be particularly noted that in the following description, when the detailed description of known functions and designs may obscure the main content of the application, these descriptions will be omitted here.

[0067] 1, Time domain signal model

[0068] Without loss of generality, the application designs a typical adaptive filter in the context of acoustic system identification. The error signal between the unknown linear system at time n and the output of the model filter can be expressed as:

[0069]

[0070] where ∈(n) is the real-valued error signal, y(n) is the observed signal, which is composed of the output d(n) = x T (n)h of the unknown system and additive noise v(n), i.e., y(n) = d(n) + v(n), h is the unknown system impulse response of length L,

[0071] x(n) = [x(n) x(n - 1)... x(n - L + 1)] T (2)

[0072] is the input signal vector composed of the last L samples, the superscript T is the transpose operator,

[0073]

[0074] is an estimate of the model filter h of length L.

[0075] In practical acoustic applications, most linear systems are essentially low-rank due to the redundancy caused by the reflections of sound waves off the walls of a room. The impulse response of the acoustic channel of such a system can be approximated by a low-rank model related to a Kronecker product of a set of shortening impulse responses. Therefore, the adaptive filter is decomposed using a Kronecker product as follows:

[0076]

[0077] where and are sub-filters of length L1and L2, respectively, denotes the Kronecker product, assuming L = L1L2, P < 2 min{L1, L2}. By exploiting the following relation:

[0078]

[0079] where and are the identity matrices of dimension L1x L1and L2x L2, respectively, the error signal in (1) can be expressed in two equivalent forms as follows:

[0080]

[0081]

[0082] where:

[0083]

[0084]

[0085]

[0086]

[0087]

[0088]

[0089] 2. Frequency-domain signal model

[0090] In order to derive the adaptive filtering algorithm in the frequency domain, the time-domain signal model needs to be represented in blocks. To this end, the error signal model in (6) and (7) is rewritten in vector form as

[0091]

[0092]

[0093] where m is the block time index,

[0094] ∈1(m) = [∈1(mL) ∈1(mL+1)... ∈1(mL+L-1)] T (16)

[0095]

[0096] X(m) = [x(mL) x(mL+1)... x(mL+L-1)] (18)

[0097] v(m) = [v(mL) v(mL+1)... v(mL+L-1)] T (19)

[0098]

[0099] Q2(m) = [Q 2，1 (m) Q 2，2 (m)... Q 2，P (m)] (21)

[0100]

[0101] Q1(m) = [Q 1，1 (m) Q 1，2 (m)... Q 1，P (m)] (23)

[0102] Note that in (14) and (15), the block length represented by and is equal to the length L of the adaptive filter to ensure good filtering performance of the inventive method.

[0103] To reduce the computational complexity of the adaptive filtering algorithm, fast block operations are needed to make the computation more efficient. Because the matrix Q 1，p (m) and Q 2，p (m) are of dimensions L x L2and L x L1, respectively, it is not possible to diagonalize them. However, the matrix X T (m) is a Toeplitz matrix and can be diagonalized as:

[0104]

[0105] where 0 L×L is a zero matrix of size L x L, R 2L×2L is a Fourier matrix of size 2L x 2L whose elements are defined as:

[0106]

[0107] denotes the imaginary unit;

[0108]

[0109] diag{•} denotes the spreading of its vector into a diagonal matrix;

[0110] x(m) = [x(mL-L) x(mL-L+1)... x(mL+L-1)] T (27)

[0111] The block error signals in (14) and (15) can then be expressed in the frequency domain as:

[0112]

[0113]

[0114] where:

[0115] ∈ 1(m) = F L×L ∈1(m) (30)

[0116] y 1(m) = F L×L y(m) (31)

[0117]

[0118]

[0119]

[0120]

[0121] ∈ 2(m) = F L×L ∈2(m) (36)

[0122]

[0123]

[0124]

[0125] From the above equations, it can be seen that the dimensions of and are 2L x PL2and 2L x PL1, respectively, while the dimension of is 2L x 2L. Generally, PL1< 2L and PL2< 2L, and the dimensions of are smaller than that of , which is beneficial to reduce the computational complexity of the adaptive filtering algorithm.

[0126] 3. Optimization criterion and adaptive algorithm

[0127] Since the present application focuses on the NKP decomposition of the frequency-domain adaptive filter, we consider a typical optimization criterion, i.e., the recursive least square (RLS) criterion, which can generate an adaptive algorithm robust to Gaussian noise. Thus, a set of cost functions based on the above frequency-domain signal model can be defined as follows:

[0128]

[0129] where λ1and λ2are the forgetting factors satisfying 0 < λ1, λ2< 1;

[0130]

[0131]

[0132] The superscript H denotes the conjugate transpose.

[0133] According to (40), the gradient with respect to can be derived as:

[0134]

[0135] where the superscript * denotes the conjugate. Let the following normal equation can be obtained:

[0136]

[0137] where:

[0138]

[0139] is the power spectrum matrix, Re{·} denotes taking the real part;

[0140]

[0141]

[0142] Similarly, another normal equation can be obtained from (40) as follows:

[0143]

[0144] where:

[0145]

[0146] is the power spectrum matrix;

[0147]

[0148] Based on (44), (45), (47) and (48)-(50), the NKP-FRLS algorithm can be derived, the flow of which is shown in Figure 1 .

[0149] Specifically, as shown in Figure 1 , the frequency domain adaptive filtering method for identifying a low-rank acoustic system according to the present application comprises the following steps:

[0150] Step S1: initialization

[0151] First, a first group of P sub-filter vectors with a length of L1 and a second group of P sub-filter vectors with a length of L2 are constructed, where m is a block time index, and then initialization is performed:

[0152]

[0153]

[0154] where the parameter η is selected in the range of greater than 0 and less than or equal to 1, and P < min{L1, L2}, so that the first and second group of sub-filter vectors with the block time index m = 0 are:

[0155]

[0156]

[0157] The two groups of power spectrum matrices with the block time index m = 0 are initialized​

[0158]

[0159]

[0160] wherein, and are identity matrices of size PL1xPL1 and PL2xPL2, respectively, and parameters δ1, δ2 are selected according to specific implementation;

[0161] Step S2: constructing an acoustic signal vector x(m)

[0162] The acoustic signal is collected to obtain a sample value x(n), and an acoustic signal vector x(m) corresponding to a block time index m is constructed:

[0163] x(m) = [x(mL-L) x(mL-L+1) … x(mL+L-1)] T

[0164] wherein, L = L1L2, and a value at a time less than or equal to 0 in the vector is 0.

[0165] Step S3: calculating a coefficient vector of an adaptive filter

[0166] Starting from a block time index m = 1, a coefficient vector of an adaptive filter is calculated

[0167] Step S3.1: calculating a spectral matrix of an input signal

[0168]

[0169] wherein, diag{·} represents expanding a vector into a diagonal matrix, F 2L×2L is a Fourier matrix of size 2Lx2L.

[0170] Step S3.2: calculating a sub-filter vector

[0171] Step S3.2.1: calculating a sub-filter spectral matrix

[0172]

[0173] wherein, represents a Kronecker product, I L×L is an identity matrix of size LxL, is an identity matrix of size L1 L×L is a zero matrix of size LxL.

[0174] Thus, the combined sub-filter spectrum matrix is obtained

[0175]

[0176] Step S3.2.2: Calculate the frequency-domain sub-filter error vector ∈ 1(m) :

[0177]

[0178] wherein, y (m) = F L×L y(m), y(m) = [y(mL) y(mL+1)... y(mL+L-1)] T is the observation signal vector corresponding to block time index m, y(n) is the observation signal sample; F L×L is the Fourier matrix of size L x L.

[0179] Step S3.2.3: Calculate the power spectrum matrix

[0180]

[0181] wherein, λ1 is a forgetting factor satisfying 0 < λ1 < 1, the superscript H represents conjugate transpose, and Re represents taking real part.

[0182] Step S3.2.4: Calculate the Kalman filter gain matrix κ2(m):

[0183]

[0184] Step S3.2.5: Calculate the sub-filter vector

[0185]

[0186] wherein, μ1 is an adaptive step size satisfying 0 < μ1 < 2,

[0187] Step S3.3: Calculate the sub-filter vector

[0188] Step S3.3.1: Calculate the sub-filter spectrum matrix

[0189]

[0190] wherein, is the unit matrix of size L2 x L2.

[0191] Thus, the combined sub-filter spectrum matrix is obtained

[0192]

[0193] Step S3.3.2: Calculate the frequency domain sub-filter error vector ∈ 2(m) = R2(m) - R1(m) * κ1(m) * R2(m) (2)

[0194]

[0195] Step S3.3.3: Calculate the power spectrum matrix

[0196]

[0197] where λ2 is a forgetting factor, satisfying 0 < λ2 < 1.

[0198] Step S3.3.4: Calculate the Kalman filter gain matrix κ1(m):

[0199]

[0200] Step S3.3.5: Calculate the sub-filter vector

[0201]

[0202] where μ2 is an adaptive step size, satisfying 0 < μ2 < 2.

[0203] Step S3.4: Calculate the filter coefficient vector

[0204]

[0205] Step S4: Filter coefficient vector Transmit filter

[0206] The filter coefficient vector obtained by the block time index m is Transmit filter, filter the sound signal.

[0207] 3. Experiment

[0208] To investigate the performance of the proposed adaptive filtering algorithm, an impulse response measured in a real room was used as the true acoustic impulse response of a linear system in the experiment. The sparsity of the impulse response was about 0.75. The sampling rate of the impulse response and the excitation speech signal was 16 kHz. The acoustic system under Gaussian white noise was identified by using the frequency-domain recursive least squares (FRLS) algorithm and the proposed NKP-FRLS algorithm, and the parameter settings were L = 1024, L1 = 512, L2 = 2, μ1 = μ2 = 0.5, λ1 = λ2 = 0.97, P ∈ [2, 4]. In the experiment, the normalized mean square deviation (NMSD) was used to evaluate the performance of the adaptive acoustic system identification algorithm. The related algorithm was implemented on a Lenovo ThinkPad P15 laptop computer running Microsoft Windows 10 operating system.

[0209] The first experiment was used to verify the convergence of the proposed algorithm when the white Gaussian sequence was used as the input signal. The acoustic impulse response of the linear system changed from h to -h at the 38th second of the input sequence, and the experimental results are shown in Figure 2 Figure 2 It can be seen from

[0210] Algorithm FRLS NKP-FRLS (the invention) Average time consumption / s 1899 708(P＝2) Average time consumption / s 1906 1845(P＝4)

[0211] Table 1

[0212] Table 1 is the time consumption of two algorithms in identifying the same acoustic system in the Lenovo ThinkPad P15 laptop computer. From Table 1, it can be seen that the proposed NKP-FRLS algorithm has higher computational efficiency than the FRLS algorithm when P = 2. Even when P = 4, the computational complexity of the NKP-FRLS algorithm is still low. From Figure 2 and the results of Table 1, it can be seen that the proposed algorithm is superior to the traditional FRLS algorithm in terms of convergence performance and computational efficiency.

[0213] The second experiment was used to verify the convergence of the proposed NKP-FRLS algorithm when speech was used as the input signal, and the experimental results are shown in Figure 3 ​It can be seen that the convergence performance of both FRLS algorithm and NKP-FRLS algorithm is reduced due to the non-stationary characteristics of speech. Compared with FRLS algorithm, the NKP-FRLS algorithm is more robust in Gaussian white noise environment, especially when P=4.

[0214] Although the foregoing describes the specific embodiments of the present application in detail, it should be clear that the present application is not limited to the scope of the specific embodiments, and all the applications that utilize the concept of the present application are within the protection scope of the present application as long as various changes are within the spirit and scope of the present application defined and determined by the appended claims.

Claims

1. A frequency-domain adaptive filtering method for identifying a low-rank acoustic system, characterized in that, The method comprises the following steps: (1) initialization First, construct the first group. A length of sub-filter And the second group Length For sub-filters , Create a block time index, then initialize it: ； ； Wherein, the parameter It can be selected in the range of greater than 0 and less than or equal to 1, In this way, the block time index The first and second groups of sub-filter vectors are: ； ； Initialization of block time index Two sets of power spectrum matrices : ； ； wherein and are identity matrices of size and respectively, the parameters δ1, δ2 being selected according to the specific implementation. (2) Collecting the sound signal to obtain a sample value and constructing a block time index a corresponding sound signal vector : ； wherein , the value of the vector at a time instant less than or equal to the time instant 0 is 0; (3) From the block time index = 1, calculate the coefficient vector of the adaptive filter ; 3.1) calculating a spectral matrix of the input signal ; ； wherein, denotes the unfolding of the vector into a diagonal matrix, is a Fourier matrix of size n x n. 3.2), computing the sub-filter vector ; 3.2.1), calculating a sub-filtered spectrum matrix ; ； wherein denotes the Kronecker product, is a unit matrix of size is a unit matrix of size is a zero matrix of size is a zero matrix of size is a zero matrix of size is a zero matrix of size Thus, the combined sub-filter spectrum matrix is obtained ; ； 3.2.2), computing the frequency domain sub-filtered error vector : ； wherein is a block time index is a corresponding observation signal vector, is an observation signal sample; is a Fourier matrix of size ​ 3.2.3), computing the power spectral matrix ; ； where λ1is a forgetting factor satisfying 0 < λ1< 1, the superscript denotes conjugate transpose, and Re denotes taking real part; 3.2.4), computing Kalman filter gain matrix ; ； 3.2.5), computing the sub-filter vector ; ； wherein is an adaptive step size satisfying , ; 3.3), computing the sub-filter vector ; 3.3.1), computing the sub-filtered spectrum matrix ; ； wherein is a unit matrix of size is a unit matrix of size Thus, the combined sub-filter spectrum matrix is obtained ; ； 3.3.2), computing the frequency domain sub-filtered error vector : ； 3.3.3), computing the power spectral matrix ; ； wherein λ2is a forgetting factor satisfying ; 3.3.4), computing Kalman filter gain matrix ; ； 3.3.5), computing the sub-filter vector ; ； wherein is an adaptive step size satisfying ; 3.4), computing the filter coefficient vector ; ； (4) block time index the resulting filter coefficient vector transmit filter, to filter the acoustic signal.

Images