Method for obtaining extremely simple convolution non-negative dictionary from spectrogram based on maximum marginal likelihood estimation

Abstract

The invention discloses a method for obtaining a very simple convolution non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation, and belongs to the technical field of convolution learning. And based on maximum marginal likelihood estimation, constructing an objective function with MMLE under a convolution non-negative framework, and based on variation with maximum expectation, inferring and solving an optimization problem so as to obtain a dictionary W. According to the method, the three-dimensional auditory dictionary W with short-time frequency spectrum information can be directly obtained from the observation data V under the condition that the activation H does not need to be estimated, and the problem of under-fitting caused by insufficient model complexity and the problem of over-fitting caused by excessive model primary functions are avoided.

Description

Technical Field

[0001] The present invention belongs to the technical field of convolutional learning, and particularly relates to a method for learning a minimalist convolutional non-negative dictionary from spectrograms based on maximum marginal likelihood estimation. Background Art

[0002] Hearing is one of the most important senses of humans. In daily life, the auditory system can encode the sound waves received by the ears into meaningful representations. Since the sensory system has only very limited memory to store the detailed information of sound waves or images, the efficient coding hypothesis in neuroscience believes that neurons efficiently encode the signals input to the senses by making full use of the statistical structure in natural signals, so as to maximize the transmission of sensory information, minimize the energy consumption, and reduce the redundancy existing in the memory.

[0003] Many studies have shown that the efficient coding hypothesis holds in the initial stages of the visual and auditory systems. The efficient coding hypothesis provides another paradigm for understanding sensory representations. The information maximization criterion in efficient coding has always been the preferred choice for research because many statistical models and methods can be used to analyze this problem under this model and criterion. From the perspective of information theory, when encoding sensory information, neural representations may have reached the optimal point under information theory, so that under the condition of hardly losing information, the redundancy of the model can be greatly reduced. Therefore, probability models can be used to model the efficient coding problem and study how the nervous system deals with the contradiction between limited resources and the pressure brought by representing sound information.

[0004] Although the hypothesis about efficient coding can be traced back decades, there are still few mathematical tools with interpretable time-frequency features for audio modeling. Because there are many possible forms of high-order statistical correlations in time-series signals, obtaining new methods that are empowered by interpretable artificial intelligence and at the same time have the effectiveness of efficient coding and model interpretability is still an open research topic. In addition, the specific dependencies in the signals are not obvious, making the work of obtaining these correlations extremely challenging. Moreover, the proposed computational models should also have the characteristic of being easy to implement in practical applications.

[0005] The core issue of efficient auditory coding is to learn a minimalist and interpretable auditory dictionary. Based on this dictionary, sounds can be analyzed, represented, and efficiently encoded in subsequent tasks. The representative mathematical model for enhancing model interpretability is non-negative matrix factorization (NMF) and its variants. Mathematically, let V represent the observed data of size F×N (each column represents a frame). The goal of NMF is to estimate two non-negative matrices W and H (sizes F×K and K×N respectively), and the product of the two is used to approximate the observed matrix V.

[0006] The columns in the dictionary W are the learned basis functions or atoms, which can be regarded as the intrinsic composition patterns of spectro-temporal data during audition. H is usually called the activation matrix, and the parameters or activations stored in each row represent the contributions of the basis functions in reconstructing the sound. Usually, a sparsity constraint is added to matrix H to reflect that only a small number of neurons are active simultaneously. In this way, W can learn features with statistical significance.

[0007] In auditory dictionary learning, NMF has been successfully extended. Specifically, convolutional non-negative sparse coding introduces a convolutional operation in W and adds a sparsity constraint to matrix H to enhance interpretability during the auditory dictionary learning process. By introducing a convolutional framework, CNSC can learn the basic spectro-temporal basis functions W. The literature refers to them as spectro-temporal kernels (STKs), which are used to simulate the spectro-temporal receptive fields (STRFs) in the early auditory nerve. CNSC has been used to learn mid-level auditory codes from the statistical features of natural sounds and to study the ecological origin of grouping principles in the auditory system. However, CNSC and other existing NMF methods still have the following problems when learning auditory dictionaries:

[0008] (1) There is sub-optimality in model definition;

[0009] (2) It is difficult to select an appropriate model order;

[0010] (3) Hyperparameters cannot be automatically estimated from the data. Summary of the Invention

[0011] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a maximum marginal likelihood estimation method for learning a minimalist convolutional non-negative dictionary based on audio spectrograms.

[0012] To achieve the above purpose, the technical solution adopted by the present invention is:

[0013] A method for learning a minimalist convolutional non - negative dictionary from spectrograms based on maximum marginal likelihood estimation, comprising the following steps:

[0014] S1. Based on joint sparsity and conjugate prior, connect the activation function and the basis function, and the expressions of the activation function and the basis function are as follows respectively:

[0015]

[0016] In the formula, is the Gamma distribution, and w(t) f,k and h k,n represent the elements in W and H respectively;

[0017] S2. Assume that the observed data V is generated by a Poisson distribution with parameters W and H. Based on the expressions of the activation function and the basis function obtained in step S1, the generation model of the observed data V is:

[0018]

[0019] In the formula, when n - t < 1, h k,n-t is zero, v f,n is an element in V, represents the Poisson distribution, and Γ(·) represents the Gamma distribution, w(t) f,k represents an element in W;

[0020] S3. Based on the generation model obtained in step S2, construct an objective function and transform the objective function into the following equation;

[0021] p(V|W) = ∫ C,H p(V,C,H|W)dCdH

[0022] In the formula, p(V|W) is the probability of V given W, and p(V,C,H|W) is the joint probability distribution given W;

[0023] S4. Use Jensen's inequality to derive the lower bound of the equation obtained in step S3, and obtain the derived lower bound. The expression of the lower bound is:

[0024] Q = ∫ C,H q(C,H)logp(V,C,H|W)dCdH

[0025] In the formula, q(C,H) is an easy - to - handle instrumental distribution, and logp(V,C,H|W) represents the logarithmic distribution to be solved;

[0026] S5. Simplify the lower boundary expression obtained in step S4 to get the following simplified formula:

[0027]

[0028] In the formula, q(C f,n ) is a polynomial distribution with probability p t,k,f,n , q(h k,n ) is a Gamma distribution, and C f,n represents a matrix of size T×K;

[0029] S6. According to the principle of the variational Bayesian EM algorithm, calculate the functional derivatives of the parameter factors q(h k,n ) and q(C f,n ) in the simplified formula in step S5, and respectively derive the update formulas for the optimal variational distributions of q(h k,n ) and q(C f,n ). The update formulas are as follows:

[0030]

[0031] In the formula, ∝ indicates that the two are equal under the condition of ignoring the constant term; the symbol <·> represents the expectation under a given probability distribution condition, c t,k,f,n represents the element in the tensor C, w(t) f,k and h k,n respectively represent the elements in W and H, Γ(·) represents the Gamma distribution, α and β are the parameters in the Gamma distribution, and α>0, β>0. When n + t > N, <c t,k,f,n+t > is zero, corresponding to the left shift operation of the matrix, and h k,n-t is the element in the matrix obtained by shifting the matrix H t units to the right along the n dimension.

[0032] And the q(h k,n ) that satisfies the optimal variational distribution has the following sufficient statistics:

[0033]

[0034] In the formula, and respectively represent the shape and scale parameters of the variational distribution q(h k,n ) that follows the Gamma distribution; c t,k,f,n represents the element in the tensor C, c t,k,f,n+t represents the element obtained by shifting the tensor C t units to the left along the n dimension, and w(t) f,k represents the element in W.

[0035] S7. Based on Bayes' rule, the lower bound of the equation obtained in step S3 is derived to obtain the lower bound of the objective function to be derived. The expression of the lower bound of the objective function is as follows:

[0036]

[0037] In the formula, Q represents the lower bound of the objective function, <logp(C|W,H)> q(C,H) represents the expectation of logp(C|W,H) under the condition of the given probability distribution q(C,H), <logp(H)> q(C,H) represents the expectation of logp(H) under the condition of the given probability distribution q(C,H), w(t) f,k and h k,n represent the elements in W and H respectively, h k,n-t is the element in the matrix obtained by shifting the matrix H t units to the right along the n dimension.

[0038] S8. According to the expression of the lower bound of the objective function obtained in step S7, the update criterion for the elements in W is obtained. The update criterion is as follows:

[0039]

[0040] In the formula, the symbol <·> represents the expectation under a given probability distribution, w(t) f,k represents the element in W, c t,k,f,n represents the element in the tensor C, h k,n represents the element in H, h k,n-t is the element in the matrix obtained by shifting the matrix H t units to the right along the n dimension.

[0041] S9. The adjustment of the hyperparameters α k and β k in the sufficient statistics obtained in step S6 is transformed into the adjustment of the shape and scale hyperparameters a and b in the inverse Gamma distribution. According to the principle of the variational Bayesian EM algorithm, the update formulas for the hyperparameters a and b are inferred, and the update formulas for the hyperparameters a and b are as follows:

[0042]

[0043] In the formula, a and b are the parameters controlling the related parameter β k of the inverse Gamma distribution, K is the number of basis functions in the dictionary W, ψ is the digamma function defined as , and represent the shape and scale parameters of the inverse Gamma distribution q(h k,n ) respectively.

[0044] S10. Automatically update the hyperparameters a and b according to the update formula obtained in step S9, and learn the dictionary W according to the updated hyperparameters a and b, the update formula obtained in step S6, and the update criterion obtained in step S8.

[0045] Preferably, step S1 includes the following steps:

[0046] S11. Combine the simultaneously occurring activation functions into a basis function;

[0047] S12. Connect the activation function and the basis function obtained in step S1 using a Gamma prior with the same hyperparameters. The expressions of the activation function and the basis function are as follows:

[0048]

[0049] In the formula, is a Gamma distribution, and and the parameters α > 0, β > 0, w(t) f,k and h k,n represent the elements in W and H respectively.

[0050] Preferably, step S3 includes the following steps:

[0051] S31. Take the likelihood function of the generative model obtained in step S2 as the objective function, and divide the activation function by integration to obtain the expression of the marginal likelihood as follows:

[0052] p(V|W) = ∫ H p(V,H|W)dH

[0053] In the formula, p(V|W,H) is the likelihood function of the generative model, and p(H) is the prior of the matrix H;

[0054] S32. Introduce a tensor C with four dimensions as an intermediate variable into the generative model obtained in step S2, and partition the observed data V, so as to decouple the coupling relationship between different basis functions and time slices. Then, transform the generative model obtained in step S2 into the following combined form:

[0055]

[0056] In the formula, c t,k,f,n is a hidden variable of the four-dimensional tensor C, w(t) f,k represents the element in W, and h k,n-t is the element in the matrix obtained by shifting the matrix H t units to the right along the n dimension.

[0057] S33. Based on step S32, transform the expression of the marginal likelihood obtained in step S31 into the following equation;

[0058] p(V|W) = ∫ C,H p(V,C,H|W)dCdH

[0059] In the formula, C is an intermediate variable with four dimensions, and p(V,C,H|W) is the joint probability distribution under the condition of given W.

[0060] Preferably, step S6 includes the following steps:

[0061] S61. The parameter factors q(h k,n ) and q(C f,n ) in the simplified formula in step S5 satisfy the following expectations:

[0062]

[0063] In the formula, <·> π represents the expectation under the probability distribution π, and the symbol A -i,j represents the elements in A except for a i,j logp(V,C,H|W) represents the logarithmic distribution to be solved;

[0064] S62. Expand logp(V,C,H|W) in step S61 to obtain the expansion formula:

[0065]

[0066] In the formula, δ(·) represents the Kronecker delta function;

[0067] S63. According to the expectation expression in step S61 and the expansion formula obtained in step S62, derive the update formula for the optimal variational distributions of q(h k,n ) and q(C f,n ), and the update formula is as follows:

[0068]

[0069] In the formula, ∝ represents the equality of the two under the condition of ignoring the constant term; the symbol <·> represents the expectation under a given probability distribution, c t,k,f,n represents the element in the tensor C, c t,k,f,n+t represents the element obtained by shifting the tensor C t units to the left along the n dimension, w(t) f,k and h k,n represent the elements in W and H respectively, Γ(·) represents the Gamma distribution, and α and β are the parameters in the Gamma distribution, and α > 0, β > 0. When \(n + t\gt N\), \(\lt c t,k,f,n+t \gt\) is zero, corresponding to the left shift operation of the matrix, \(h k,n-t is an element in the matrix obtained by shifting the matrix \(H\) \(t\) units to the right along the dimension of \(n\).

[0070] And \(q(h k,n )\) that satisfies the optimal variational distribution has the following sufficient statistics:

[0071]

[0072] In the formula, and respectively represent the shape and scale parameters of the variational distribution \(q(h k,n )\) that follows the Gamma distribution; \(c t,k,f,n represents an element in the tensor \(C\), and \(c t,k,f,n+t represents the element obtained by shifting the tensor \(C\) \(t\) units to the left along the dimension of \(n\), and \(w(t) f,k represents an element in \(W\).

[0073] Preferably, step S8 includes the following steps:

[0074] S81. According to the expression of the lower bound of the objective function obtained in step S7, calculate the partial derivative of the lower bound with respect to the dictionary \(W\) to obtain the following relationship:

[0075]

[0076] In the formula, \(w(t) f,k represents an element in \(W\), \(c t,k,f,n represents an element in the tensor \(C\), and \(h k,n-t is an element in the matrix obtained by shifting the matrix \(H\) \(t\) units to the right along the dimension of \(n\).

[0077] S82. According to the relationship obtained in step S81, obtain the update criterion for the elements in \(W\), and the update criterion is as follows:

[0078]

[0079] In the formula, the symbol \(\lt\cdot\gt\) represents the expectation under a given probability distribution condition, \(w(t) f,k represents an element in \(W\), \(c t,k,f,n represents an element in the tensor \(C\), and \(h k,n-t is an element in the matrix obtained by shifting the matrix \(H\) \(t\) units to the right along the dimension of \(n\).

[0080] Preferably, step S9 includes the following steps:

[0081] S91. Comply with the prior of ARD NMF and set \(\alpha k= 1, transform the Gamma distribution into an exponential distribution, for β k Add a prior of inverse Gamma distribution, then the hyperparameters α in the sufficient statistic obtained in step S6 k and β k The adjustment of is transformed into the adjustment of the shape and scale hyperparameters a and b in the inverse Gamma distribution;

[0082] S92. Estimate the hyperparameters a and b in step S91 to obtain the likelihood function p(V|a,b), and the expression is as follows:

[0083]

[0084] In the formula, p(V|a,b) is the likelihood function of V under the condition of given parameters a and b. This likelihood function is obtained by integrating and dividing out the parameters C, H, W and in the likelihood function p(V,C,H,W,β|a,b).

[0085] S93. Based on the likelihood function obtained in step S92, use Jensen's inequality to derive the lower bound of the log-likelihood function as the new objective function. The expression of the new objective function is:

[0086] logp(V|a,b)≥∫ C,H,W,β q(C,H,W,β)

[0087]

[0088] logp(V,C,H,W,β|a,b)dCdHdWdβ

[0089] In the formula, q(C,H,W,β) is the variational distribution;

[0090] S94. According to the principle of the variational Bayesian EM algorithm, transform the variational distribution q(C,H,W,β) in step S93 into the following decomposable form:

[0091]

[0092] In the formula, q(C f,n ) is a matrix where the random variable follows a multinomial distribution, q(w(t) f,k ) and q(h k,n ) follow the Gamma distribution and share the same parameters, and q(β k ) is the inverse Gamma distribution;

[0093] Based on step S94, the optimal variational distribution is obtained, and then the lower bound is calculated. According to the lower bound, the update formulas for hyperparameters a and b are inferred to obtain the update formulas for hyperparameters a and b, as shown below:

[0094]

[0095] In the formula, a and b are parameters that control the relevant parameter β of the inverse Gamma distribution, K is the number of basis functions in the dictionary W, ψ is the digamma function defined as k and K is the number of basis functions in the dictionary W, ψ is the digamma function defined as , and respectively represent the shape and scale parameters of the inverse Gamma distribution q(h k,n ).

[0096] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0097] 1. The present invention solves the sub-optimal problem in estimating the auditory dictionary W in convolutional non-negative coefficient coding, and can directly obtain the three-dimensional auditory dictionary W with short-time spectral information from the observed data V without estimating the activation H.

[0098] 2. The present invention solves the problem that the number of convolutional non-negative basis functions cannot be automatically obtained, avoiding the underfitting problem caused by insufficient model complexity and the overfitting problem caused by too many model basis functions.

[0099] 3. The present invention can continuously iteratively estimate the important hyperparameters in the assumed model directly from the data, avoiding various deficiencies caused by the manual setting method. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 It is a diagram for approximating the cochlear spectrum using convolutional non-negative sparse coding;

[0101] Among them, the small blocks in W represent the learned auditory basis functions. The sparsity of the activation function can be observed in the activation function H;

[0102] Figure 2 It is an intuitive comparison diagram of minimal representation and redundant representation; (a) is the minimal representation, (b) is the redundant representation, and the white empty squares in the figure represent zero coefficients, while the colored squares represent non-zero elements;

[0103] Figure 3 It is an intuitive diagram of the generative model with the intermediate variable C;

[0104] Figure 4Probabilistic graphical models of the algorithm proposed in the present invention and the ARD NMF algorithm; (a) is the probabilistic graphical model of ARD NMF, and (b) is the probabilistic graphical model of the algorithm proposed in the present invention;

[0105] Figure 5 It is a graph of the synthetic dataset;

[0106] Figure 6 Results learned by the algorithm proposed in the present invention; (a) is the learned convolutional basis function, (b) is the change curve of the relevant parameters during 1000 iterations, and (c) is the value of the relevant parameters when the iteration terminates;

[0107] Figure 7 It is a visualization graph of the convolutional basis function learned by the baseline method; (a) is the basis function obtained by Algorithm 1, (b) is the basis function obtained by Algorithm 2, (c) is the basis function obtained by Algorithm 3, (d) is the basis function obtained by Algorithm 4, (e) is the basis function obtained by Algorithm 5 when λ = 0.005, and (f) is the basis function obtained by Algorithm 5 when λ = 0.003;

[0108] Figure 8 They are the notes of "Mary had a little lamb";

[0109] Among them, the activated notes are E 4 , D 4 , C 4 , D 4 , E 4 , E 4 , E 4 ;

[0110] Figure 9 It is a comparative study of the algorithm proposed in the present invention and DR-NMF of Algorithm 6 under the condition of K = 6;

[0111] Figure 10 It is a comparative study of the algorithm proposed in the present invention and the baseline algorithm with convolutional basis functions; (a) to (e) respectively represent the results of Algorithms 1 to 5, and (f) represents the result of the algorithm proposed in the present invention;

[0112] Figure 11 It is the result obtained by running the algorithm proposed in the present invention on the cochlear spectrum of the sound "Siren" signal;

[0113] Among them, the left sub-graph shows the cochlear spectrum of the "Siren" signal, the middle sub-graph shows the numerical change of the relevant parameters during 1000 iterations, and the right sub-graph shows the value of the finally stable relevant parameters after 1000 iterations;

[0114] Figure 12are basis functions learned from the cochlear spectrum of the sound signal “Siren”; (a) is the STK result learned by the algorithm proposed in the present invention, and (b)-(f) are the STK results learned by algorithms 1 to 5 respectively;

[0115] Figure 13 are basis functions learned from the cochlear spectrum of the sound signal “car alarm”; (a) is the STK result learned by the algorithm proposed in the present invention, and (b)-(f) are the STK results learned by algorithms 1 to 5 respectively;

[0116] Figure 14 are the results obtained from the unsupervised sound source separation task; (a) is the cochlear spectrum obtained by mixing “hammering” and “woman speaking”, (b) are the basis functions learned from the mixed signal, (c) is the change curve of relevant parameter values during the iteration process, (d) is the cochlear spectrum of the real “hammering” signal in the mixed sound, (e) is the cochlear spectrum reconstructed using the basis functions corresponding to the maximum correlation parameter, (f) is the cochlear spectrum corresponding to the real “woman speaking” signal in the mixed sound, and (g) is the cochlear spectrum reconstructed using the remaining basis functions;

[0117] Figure 15 are the results obtained from the supervised sound source separation task; (a) is the cochlear spectrum mixed with “hammering” and “woman speaking”, (b) are the basis functions obtained after removing the irrelevant basis functions from the model. The basis functions in the red box correspond to the basis functions learned from the “hammering” signal, and the basis functions in the orange box correspond to the basis functions learned from the “woman speaking” signal; (c) is the cochlear spectrum of the real “hammering” signal in the mixed signal; (d) is the cochlear spectrum reconstructed using the basis functions in the red box; (e) is the cochlear spectrum of the real “woman speaking” signal in the mixed signal; (f) is the cochlear spectrum reconstructed using the basis functions in the orange box. Detailed implementation manners

[0118] The technical solutions in the embodiments of the present invention will be described clearly and completely below. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0119] I. Specific construction process of the method provided by the embodiments of the present invention

[0120] In the embodiment of the present invention, joint sparsity and conjugate prior are introduced in the process of learning an auditory dictionary under the convolutional non-negative framework, and then an objective function with MMLE under the convolutional non-negative framework is constructed based on MMLE. Finally, the optimization problem is solved by variational inference based on expectation-maximization (EM). By Figure 1 It can be seen that the basis function has a clear time-frequency spectrum structure. In the convolution operation, it is shifted and scaled along the time axis of the i-th row in H.

[0121] The embodiment of the present invention provides a method for learning an extremely simple convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation, which specifically includes the following steps:

[0122] 1. Introduction of joint sparsity and conjugate prior

[0123] Figure 2 An intuitive comparison diagram of the extremely simple representation and the redundant representation is given. From Figure 2 it can be seen that, compared with the redundant representation, the extremely simple representation has a better effect in reflecting the internal structure of the observation. In Figure 2 (b), the activation functions that often appear simultaneously (such as the activation functions of #1 and #2, #3 and #4) will be combined into a basis function, similar to Figure 2 (a), which is consistent with the grouping principle in the auditory system. Generally, this kind of sparsity is called the collaborative sparse of the convolutional non-negative basis function.

[0124] Inspired by NMF with ARD, the activation function and the basis function are connected together using a Gamma prior with the same hyperparameters.

[0125]

[0126] In the formula, is the Gamma distribution, and β > 0, α > 0, β > 0, w(t) f,k and h k,n represent the elements in W and H respectively.

[0127] Because the Gamma distribution is conjugate to the Poisson distribution, the Gamma distribution is the preferred choice. If the value of α k is fixed to 1, the Gamma distribution can be transformed into an exponential distribution. In ARD NMF, all β k are called correlation parameters.

[0128] 2. Construction of a new objective function

[0129] From the perspective of the generative model, assuming that the observed data V is generated by a Poisson distribution with parameters W and H, the generative model of the observed data V is as follows:

[0130]

[0131] In the formula, when n - t < 1, h k,n-t is zero, and zero-padding is performed on the elements in the matrix shift operation. The elements in V are v f,n , and the symbol represents the Poisson distribution, and Γ(·) represents the Gamma distribution, and w(t) f,k represents the elements in W.

[0132] Taking the marginal likelihood as the new objective function and integrating out all possible activation functions, the expression of the marginal likelihood is as follows:

[0133] p(V|W) = ∫ H p(V,H|W)dH (14)

[0134] In the formula, p(V|W) is the probability of V given W, and p(V,C,H|W) is the joint probability distribution given W;

[0135] As Figure 3 shown, the matrix V is divided into several time slices along the dimension of t of C, and is used to approximate each time slice C(t).

[0136] Therefore, formula (14) can be transformed into the maximization problem of the following formula under the given dictionary W:

[0137] p(V|W) = ∫ C,H p(V,C,H|W)dCdH (18)

[0138] In the formula, p(V|W) is the probability of V given W, and p(V,C,H|W) is the joint probability distribution given W.

[0139] 3. Maximum Marginal Likelihood Estimation

[0140] To maximize formula (18), the lower bound to be maximized is derived based on variational inference. The Expectation-Maximization (EM) algorithm is used to derive the update algorithm for optimizing the lower bound.

[0141] (1) Variational Inference

[0142] To maximize the marginal likelihood of formula (18), Jensen's inequality is used to derive its lower bound:

[0143]

[0144] where ∝ denotes equality without considering the constant term, and q(C, H) is the variational distribution used to approximate the exact posterior distribution.

[0145] Therefore, the lower bound to be maximized has the following form:

[0146] Q = ∫ C,H q(C, H) log p(V, C, H|W) dC dH (20)

[0147] where q(C, H) is an intractable instrumental distribution, and log p(V, C, H|W) represents the log distribution to be solved;

[0148] When the variational distribution is equal to the posterior distribution, the lower bound Q is tight. To simplify the integral in Equation (20), mean-field variational inference assuming all variables in q(C, H) are independent is used to simplify the integral as follows:

[0149]

[0150] where q(C f,n ) is a multinomial distribution with probability p t,k,f,n , q(h k,n ) is a Gamma distribution, and its estimated shape and scale parameters are and C f,n represents a matrix of size T×K, and the expression of C f,n is as follows:

[0151]

[0152] where c T-1,K,f,n represents the element in the tensor C with four dimensions.

[0153] (2) E step

[0154] According to the principle of the variational Bayesian EM algorithm, the functional derivatives of q(H) and q(C) are calculated to derive the update formula for updating the variational distribution. The factors q(h k,n ) and q(C f,n ) should satisfy the following expectations:

[0155]

[0156] where <·> π denotes the expectation under the probability distribution π, and the symbol A -i,j denotes the elements in A excluding ai,j The remaining elements outside, logp(V, C, H|W) represents the logarithmic distribution to be solved. To calculate the expectations of equations (23) and (24), logp(V, C, H|W) is expanded in the following form:

[0157]

[0158] where δ(·) represents the Kronecker delta function. In equation (25), it can be seen that the first two terms in the second row control the likelihood of the data, which corresponds to the generalized Kullback-Leibler divergence (KLD), and the third term corresponds to the regularization term derived from the prior probability assumption of H. If we observe the prior probability logp(H), we can find that -h k,n / β k monotonically decreases as β k increases. And when α k is a positive number, -α k logβ k monotonically increases as β k increases. Therefore, when the algorithm converges, a part of β will remain at a large value, and the remaining part will become very small. This property can remove the basis functions that are irrelevant in representing the data from the model through pruning.

[0159] Based on equations (23), (24), and (25), the following optimal variational distribution can be derived:

[0160]

[0161] where ∝ means equality up to a constant factor; the symbol <·> represents the expectation under a given probability distribution, c t,k,f,n represents the elements in the tensor C, w(t) f,k and h k,n represent the elements in W and H respectively, Γ(·) represents the Gamma distribution, α and β are the parameters in the Gamma distribution, and α > 0, β > 0. When n + t > N, <c t,k,f,n+t > is zero, corresponding to the left shift operation of the matrix. h k,n-t is the element in the matrix obtained by shifting the matrix H t units to the right along the n dimension, and ρ(n) is defined as follows:

[0162]

[0163] Because in the convolutional structure, the activation function has a shared property, Equations (26), (27), and (28) need to be carefully considered during the derivation process. This problem can be considered by forwardly examining the relationship between each column in H and the four-dimensional tensor C. Specifically, the optimal variational distribution q(C f,n ) that satisfies Equation (27) can be obtained through the following probability parameters:

[0164]

[0165] where w(τ) f,l and w(t) f,k represent the elements in W. Different symbols are used to avoid possible confusion during summation. w(t) f,k represents the elements in W, and h k,n-t refers to the element in the matrix obtained by shifting matrix H t units to the right along the n dimension. h l,n-τ refers to the element in the matrix obtained by shifting matrix H τ units to the right along the n dimension.

[0166] Similarly, the optimal variational distribution q(h k,n ) that satisfies Equation (26) has the following sufficient statistics:

[0167]

[0168]

[0169] where and respectively represent the shape and scale parameters of the variational distribution q(h k,n ) that follows the Gamma distribution; c t,k,f,n represents the element in tensor C, and c t,k,f,n+t represents the element obtained by shifting tensor C t units to the left along the n dimension. w(t) f,k represents the element in W.

[0170] The expectations <c t,k,f,n > in Equation (30) and <logh k,n > in Equation (29) can be obtained through the following method:

[0171]

[0172] where ψ is the digamma function defined as .

[0173] From the above content, it can be found that in Equations (30) and (31), the challenging problem of inferring sufficient statistics in the convolutional non-negative framework is solved through forward updates.

[0174] (3) Step M

[0175] Based on Bayes' rule, the objective function can be rewritten in the following form:

[0176]

[0177] where Q represents the lower bound of the objective function, <logp(C|W,H)> q(C,H) represents the expectation of logp(C|W,H) under the given probability distribution q(C,H), <logp(H)> q(C,H) represents the expectation of logp(H) under the given probability distribution q(C,H), w(t) f,k and h k,n represent the elements in W and H respectively, and h k,n-t refers to the element in the matrix obtained by shifting matrix H t units to the right along the n dimension.

[0178] Derive the lower bound of the MMLE objective function as shown in Equation (34). Calculate the partial derivative of the lower bound with respect to the dictionary W, and the following relationship can be obtained:

[0179]

[0180] where w(t) f,k represents the element in W, c t,k,f,n represents the element in the tensor C, and h k,n-t refers to the element in the matrix obtained by shifting matrix H t units to the right along the n dimension;

[0181] Therefore, the update criterion for the elements in W is as follows:

[0182]

[0183] where the symbol <·> represents the expectation under a given probability distribution, w(t) f,k represents the element in W, c t,k,f,n represents the element in the tensor C, h k,n-t refers to the element in the matrix obtained by shifting matrix H t units to the right along the n dimension, and Equations (30) and (31) respectively give and the definitions of. If Equations (29) and (32) are substituted into Equation (36), it can be found that the update criterion formula has a multiplicative iterative form.

[0184] III. Inference of Hyperparameters

[0185] Note that in Equations (30) and (31), αk and β k are both hyperparameters, and there is currently no good way to select these parameters. In the setting of the present invention, by following the prior of ARD NMF, setting α k = 1, the Gamma distribution is transformed into an exponential distribution, which is consistent with l 1 -ARD NMF in the literature. In addition, an inverse Gamma distribution prior given in Equation (10) is added for β k . Therefore, adjusting the hyperparameters changes from the original adjustment of α k and β k to determining the parameters of the shape and scale in the inverse Gamma distribution, which are the hyperparameters a and b respectively.

[0186] The selection of these two hyperparameters is very important. Existing literature proposes a method for determining the parameters based on the moment method. However, when a relatively large a is selected, it is difficult for the algorithm to obtain good results. Existing research shows that variational inference can achieve comparable results to the MCMC method, but the latter MCMC method consumes more computing resources. Therefore, variational inference is selected to solve the hyperparameter estimation problem.

[0187] Comparison between the algorithm of the present invention and the algorithms in the literature is as Figure 4 shown, Figure 4 (a) represents the probability graph model of ARD NMF in the literature, Figure 4 (b) is the probability graph model of the algorithm proposed by the present invention, where each element in W(t) is denoted as w(t) f,k . As mentioned before, the four-dimensional tensor C (the elements of which are denoted as c t,k,f,n ) is used to partition the observed data V (each element of which is denoted as v f,n ).

[0188] 1. Problem description

[0189] Before solving for W and H, first estimate a and b in the probability graph model in Figure 4 (b). The problem is modeled as follows:

[0190]

[0191] By integrating out all the random variables W, H, C and the related parameter β, the following likelihood p(V|a, b) can be obtained, and the expression is as follows:

[0192]

[0193] The priors in W and H are given in Equations (11) and (12) respectively. By using Jensen's inequality, the lower bound of the log-likelihood function is derived as the new objective function, and the expression of the new objective function is:

[0194]

[0195] where q(C, H, W, β) is the factorized variational distribution given by Equation (40). The bound is tight when this distribution is equal to the posterior probability p(C, H, W, β|V, a, b).

[0196] 2. E-step

[0197] The variational distribution can be written in the following factorized form:

[0198]

[0199] where q(C f,n ) is a matrix where the random variable follows a multinomial distribution, q(w(t) f,k ) and q(h k,n ) follow Gamma distributions and share the same parameters. q(β k ) is an inverse Gamma distribution.

[0200] To calculate the above expectations, the following derivations are carried out:

[0201]

[0202] This will result in the following variational distribution:

[0203]

[0204]

[0205] where c t,k,f,n represents the element in tensor C, c t,k,f,n+t represents the element obtained by shifting tensor C t units to the left along the n dimension, w(t) f,k represents the element in W, h k,n-t is the element in the matrix obtained by shifting matrix H t units to the right along the n dimension, h k,n represents the element in H, and a and b are the parameters controlling the parameter β k that follows an inverse Gamma distribution.

[0206] 3. M-step

[0207] According to the variational distribution, the lower bound is calculated as follows:

[0208]

[0209] In the formula, <logβ k > is equal to ψ is defined as the digamma function, and a and b are parameters controlling the parameter β related to the inverse Gamma distribution k To maximize the lower bound, a heuristic algorithm commonly used in the NMF framework is resorted to. The gradient of Q is expressed as the difference between the positive part and the negative part as follows:

[0210]

[0211] Then the heuristic algorithm that maximizes the bound can be written as:

[0212]

[0213] The multiplicative iteration algorithm that ensures non-negativity of the parameters during the update process is as shown in Equation (65). According to this formula, the derivatives of Q with respect to a and b are calculated respectively, and the calculation formulas are as follows:

[0214]

[0215] In the formula, a and b are parameters controlling the parameter β related to the inverse Gamma distribution k , K is the number of basis functions in the dictionary W, ψ is defined as the digamma function, and respectively represent the shape and scale parameters of the inverse Gamma distribution q(h k,n ).

[0216] Therefore, the following update formula can be obtained by maximizing the lower bound:

[0217]

[0218] In the formula, a and b are parameters controlling the parameter β related to the inverse Gamma distribution k , K is the number of basis functions in the dictionary W, ψ is defined as the digamma function, and respectively represent the shape and scale parameters of the inverse Gamma distribution q(h k,n ).

[0219] At this time, the key hyperparameters a and b can be automatically updated. Using the estimated a and b, the dictionary W can be learned according to the Figure 4 probabilistic graphical model in (b).

[0220] IV. Experiments

[0221] 1. Experimental Setup

[0222] The method provided by the embodiment of the present invention is compared with six existing baseline algorithms, and the details of the six baseline algorithms are listed in Table 1 below.

[0223] Table 1 Baseline Algorithms for Comparison

[0224]

[0225]

[0226] Literature [1] W. J.H. McDermott, “Ecological origins of perceptual grouping principles in the auditory system,” Proceedings of the National Academy of Sciences, vol. 116, no. 50, pp. 25355 - 25364, Dec. 2019.

[0227] Literature [2] D. Fagot, H. Wendt, C. Févotte and P. Smaragdis, “Majorization - minimization algorithms for convolutive NMF with the beta - divergence,” In Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 8202 - 8206, Brighton, United Kingdom, 2019.

[0228] Literature [3] J.H. McDermott, E.P. Simoncelli. “Sound texture perception via statistics of the auditory periphery: evidence from sound synthesis,” Neuron, vol. 71, no. 5, pp. 926 - 940, 2011.

[0229] Reference [4] E.L. Mackevicius, A.H. Bahle, A.H. Williams, et al. “Unsupervised discovery of temporal sequences in high-dimensional datasets, with applications to neuroscience,” Elife, vol. 8, no. e38471, pp. 1-42, 2019.

[0230] Reference [5] V. Renkens and H. Van hamme, “Automatic relevance determination for nonnegative dictionary learning in the gamma-poisson model,” Signal Processing, vol. 132, pp. 121-133, 2017.

[0231] 2. Evaluation on the synthetic dataset

[0232] The performance of the algorithm proposed in the embodiment of the present invention is evaluated on Figure 5 the synthetic dataset shown in the figure. This dataset has an obvious structure. If the horizontal axis is regarded as the time axis and the vertical axis is regarded as the frequency axis, it can be regarded as an amplitude spectrum with an inherent time-frequency structure.

[0233] It can be seen through Figure 5 that the synthetic data has three repeated structural patterns. Each pattern spans 20 points in the horizontal direction. The values of each structured pattern vary from 7 to 13, and the values of the remaining elements are all zero. The one with a zigzag roof shape is called type I, and the monotonically increasing and monotonically decreasing ones are called type II and type III respectively. Type I and type II are repeated four times, while type III is only repeated twice.

[0234] Set T = 20 and K = 10, and run the algorithm proposed in the embodiment of the present invention on the above synthetic dataset to test whether the algorithm proposed by the present invention can capture short-time structures as a dictionary for efficient encoding of synthetic data. Initialize a and b with 1, and initialize W and H with random non-negative values. Iterate the proposed algorithm 1000 times. The obtained results are presented in Figure 6 the figure.

[0235] It can be seen through Figure 6It can be seen that after 1000 iterations, the estimated correlation parameters become stable, and most of the correlation coefficients therein become small values close to zero. In addition, through Figure 6 it can be found that the values of three correlation parameters are significantly larger than the remaining values, which is consistent with the actual model order. Further, it can be found that the correlation parameters of type I and type II are similar because both are repeated four times. The correlation parameters of type III are smaller than those of type II because type III is repeated twice while type II is repeated four times. Therefore, the proposed algorithm can automatically learn the repeated time-frequency basis functions and remove redundancy. The obtained basis functions are extremely simple and can avoid overfitting. In addition, the estimated correlation parameters can be used to indicate the contribution of each basis function in representing the observed data.

[0236] For the purpose of comparison, the baseline algorithms listed in Table 1 are also iterated 1000 times. For the fairness of comparison, like the algorithm proposed in the present invention, the parameters of the baseline algorithms are also set to K = 10 and T = 20. The basis functions learned by Algorithms 1 to 4 are respectively shown in Figure 7 (a) to Figure 7 (d). The hyperparameter λ in Algorithm 5 is used to control the balance between data fitting and overfitting. Two different values of λ are manually selected, that is, when λ = 0.005 and λ = 0.003 respectively, the obtained basis functions are respectively shown in Figure 7 (e) and Figure 7 (f).

[0237] As can be seen from Figure 7 , when K is greater than the true value, the baseline method fails to discover the potential hidden components, but instead finds some overlapping patterns of the true patterns. The baseline algorithms (Algorithms 1 to 4) tend to obtain redundant representations spanning all basis functions. Algorithm 5 has the opportunity to obtain an extremely simple representation. However, its performance is highly sensitive to the choice of hyperparameters. In contrast, the algorithm proposed in the present invention reflects the internal structure of the observed data. And the algorithm proposed in the present invention can learn highly interpretable repeated basis functions, which is consistent with the efficient coding hypothesis being able to achieve a good balance between data fitting and sparsity. The key parameters a and b can be automatically estimated, and all other variables can be integrated out by integration.

[0238] 3. Note discovery in piano sounds

[0239] This part of the content studies the data in the literature (N. Gillis, L. T. K. Hien, V. Leplat, and V. Y. F. Tan, “Distributionally robust and multi-objective nonnegative matrix factorization,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 8, pp. 4052-4064, 2022.). The signal is from the first 4.7 seconds of the music “Mary had a little lamb” and consists of three notes, namely E4, D4, and C4. First, the signal is downsampled to 16 kHz, and a 512-point Hamming window with 50% overlap is used to calculate the short time Fourier transformation (STFT). The piano_Mary dataset is as Figure 8 shown.

[0240] In this part of the content, the number of basis functions of all algorithms is set to K = 6, which is twice the number of real notes. The obtained results are shown in Figure 9 and Figure 10 . Considering that Algorithm 6 does not have a convolutional structure capable of capturing the short-time dynamic characteristics of sound, T in the proposed algorithm is set to 1 for fair comparison. In addition, the l 1 -norm value of the rows in H obtained by Algorithm 6 indicates the contribution of the basis functions in W. The l 1 -norm values are arranged in descending order to evaluate the contribution of each basis function in representing the input data, as shown in Figure 9 shown.

[0241] When T = 4, the results of the algorithm to be evaluated with convolutional basis functions are as shown in Figure 10 shown.

[0242] For the two choices of T, the algorithm proposed in the present invention can always obtain better results than the baseline algorithm. Obviously, the algorithm proposed in the present invention can always discover the correct basis functions consistent with the number of real notes by using relevant parameters. In addition, the relevant parameters corresponding to E 4 , D 4 and C 4 have good consistency with the number of times these notes repeat.

[0243] In contrast, the baseline algorithm tends to split a single note into multiple notes, which goes against the hypothesis of efficient coding. Algorithms such as Algorithm 4, Algorithm 5, and Algorithm 6 tend to regard the initial part of each note (i.e., the mechanical vibration in the piano when a specific note is triggered) as an independent basis function, but this is not the simplest mode of the dictionary.

[0244] 4. Evaluation on Real-World Audio Datasets

[0245] This section will present the experiments on the real-world audio dataset used in the literature (

[39] S. Norman-Haignere, N. G. Kanwisher, J. H. McDermott, “Distinct cortical pathways for music and speech revealed by hypothesis-free voxel decomposition,” Neuron, vol. 88, no. 6, pp. 1281-1296, 2015.). This dataset contains 165 different natural sounds that people often hear in daily life. Each sound lasts for 2 seconds. All signals are downsampled from 44.1 kHz to 16 kHz.

[0246] Learning an auditory dictionary from real-world datasets is more challenging than from synthetic data because the number of true basis functions is unknown. To address this issue, sounds with an obvious repetitive structure are selected to compare the proposed algorithm and the baseline algorithm.

[0247] The cochlear spectrogram, which can simulate the mammalian cochlea and provide a rough model of the auditory nerve input to the brain, is used for time-frequency analysis. A cochlear spectrogram with 128 channels is used, equally spaced on the equivalent rectangular bandwidth (ERB) scale. The center frequencies are distributed in the range of 80 Hz to 8 kHz. The window length is 320 points, equivalent to 20 ms at a sampling rate of 16 kHz. The frame shift is half of the window length, i.e., 10 ms.

[0248] First, the algorithm proposed in the present invention is run on the "siren", i.e., the "Siren" sound for testing, and the results are as Figure 11 shown. As can be seen from Figure 11, this sound has a quasi-periodic short-time structure. Set K = 5 and T = 36 (i.e., 370 ms) for the algorithm proposed in the present invention and the baseline algorithm, and all algorithms run 1000 iterations. The obtained STK results are as Figure 12 shown. The changes in the values of the relevant parameters during the 1000 iterations and their final results are respectively in Figure 11It is presented in the middle sub - figure and the right - hand sub - figure.

[0249] From Figure 12 (a), it can be seen that the algorithm proposed in the present invention only obtains one main basis function, and the relevant parameters corresponding to this basis function are significantly larger than other relevant parameters. When ignoring the subtle differences between quasi - periodic time structures, this is consistent with auditory perception. As a comparison, the baseline algorithms (Algorithms 1 to 4) obtain similar results, and the visualization results of the relevant learned basis functions are given in Figure 12 (b) - (e). From Figure 12 (f), it can be seen that the result obtained by Algorithm 5 is similar to the algorithm proposed in the present invention, but its effect depends to a large extent on the choice of hyperparameters.

[0250] Running the algorithm proposed in the present invention on "car alarm", that is, on the cochlear spectrum corresponding to the "car alarm" sound signal, from Figure 13 (a), it can be seen that this signal also has a quasi - periodic time - frequency structure. Set the parameters K = 12 and T = 16 (i.e., 170 ms) for the algorithm proposed in the present invention and Algorithm 1, and run 1000 iterations. Display the basis functions and estimated activation functions learned by the proposed algorithm and the baseline algorithm from the cochlear spectrum of the "car alarm" signal. In addition, according to the values of the relevant parameters after the algorithm converges, calculate the activation functions corresponding to the largest 4 relevant parameters, and calculate their cross - correlation functions. The obtained results are presented in Figure 13 as follows.

[0251] It can be seen that the basis function dictionary obtained by the proposed algorithm is more intuitive compared to the baseline algorithm. The probability that the basis functions generated by the proposed algorithm appear simultaneously is lower, and the representation form is also more concise. The cross - correlation function between the activation functions corresponding to the largest 4 relevant parameters reveals that the observed data is a sound with a periodic structure.

[0252] 5. Application in source separation

[0253] In this part, test whether the algorithm proposed in the present invention can extract meaningful time - frequency features from the mixed sound source signal, and demonstrate that the algorithm proposed in the present invention is meaningful in both unsupervised and supervised source separation tasks.

[0254] In the unsupervised separation task, two real - world sound signals: "hammering" (sound source 1) and "woman speaking" (sound source 2) are mixed to form an observed signal (i.e., Figure 14(as shown in (a)). The cochlear spectrum of the mixed signal has parameter settings consistent with previous experiments. Set K = 12 and T = 16 (i.e., corresponding to a sound duration of 170 ms), and iterate the algorithm proposed by the present invention 1000 times. The resulting results are given in Figure 14 as follows.

[0255] As can be seen from Figure 14 (b) and Figure 14 (c), after 1000 iterations, the algorithm proposed by the present invention can learn the basis functions of 6 STKs, with correlation parameters significantly greater than those of other basis functions. If the cochlear spectrum is reconstructed using the basis functions with large correlation parameters, the resulting results (i.e., Figure 14 (e)) are consistent with the actual "hammering" signal (i.e., Figure 14 (d)).

[0256] As a control, the cochlear spectrum reconstructed using the remaining basis functions (i.e., Figure 14 (g)) is consistent with the real "woman speaking" (i.e., Figure 14 (f)). Note that the above source separation is performed in an unsupervised manner without pre-given training data. The simultaneously occurring basis functions tend to be combined into a basis function that can better reflect the internal background noise. Through the correlation parameters, it is very convenient to locate and extract the "hammering" component, whose appearance frequency domain is much higher than the basis functions in the speech signal "woman speaking".

[0257] As a comparison, since the baseline algorithm does not have correlation parameters, its activation function is used as a substitute. Similar to the algorithm proposed by the present invention, assume that the most active basis function is the "hammering" (sound source 1) component, and the remaining part is regarded as the "woman speaking" (sound source 2) component. Calculate the output SNR of the reconstructed spectra of sound source 1 and sound source 2, and the corresponding results are given in Table 2 below.

[0258] Table 2 Output signal-to-noise ratio in unsupervised separation tasks

[0259]

[0260]

[0261] From the results in Table 2, it can be seen that the algorithm proposed by the present invention has better effects compared to the baseline algorithm. In addition, the algorithm proposed by the present invention only requires half of the memory overhead of Algorithms 1 to 4 because the remaining basis functions are not relevant to the observations, which can be seen from Figure 14 (b) and Figure 14Observed in (c). In contrast, Algorithm 5 only uses 3 basis functions, but as can be seen from Table 2, its effect is inferior to the algorithm proposed by the present invention.

[0262] In the unsupervised separation task, the same signals as in the previous experiment are used. First, 30 basis functions with T = 4 (i.e., 50 ms) are used, and the corresponding training signals are learned by the algorithm proposed by the present invention or the baseline algorithm. Then, the cochlear spectra of the mixed signals are independently projected onto the datasets of these two basis functions. Finally, the output SNR of the estimated "hammering" and "woman speaking" reconstructed spectra is calculated respectively. The obtained results are given in Table 3.

[0263] Table 3 Output Signal-to-Noise Ratio in the Supervised Separation Task

[0264]

[0265] By using relevant parameters to indicate the absorption of each basis function, after removing the irrelevant basis functions, the algorithm proposed by the present invention has learned a total of 13 basis functions, while Algorithms 1 to 4 have obtained approximately 60 basis functions. Therefore, by using only 21.67% of the memory of the baseline algorithm, the algorithm proposed by the present invention can achieve better results than the baseline algorithm. Algorithm 5 has only learned 9 basis functions (less than the algorithm proposed by the present invention), but its effect is far worse than that of the algorithm proposed by the present invention. Therefore, the performance of the algorithm proposed by the present invention is superior to all these baseline algorithms.

[0266] Although the embodiments of the present invention have been shown and described, those of ordinary skill in the art can understand that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the claims and their equivalents.

Claims

1. A method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation, characterized in that: The following steps are involved: S1. Based on joint sparsity and conjugate prior, the activation function and basis function are connected together, and the expressions of the activation function and basis function are as follows: In the formula, is a Gamma distribution, and x≥0,α>0,β>0,w(t) f,k and h k,n represent the elements in W and H respectively; S2. Assume that the observed data V is generated by a Poisson distribution with parameters W and H. Based on the expressions of the activation function and basis function obtained in step S1, the generative model of the observed data V is obtained as follows: In the formula, when nt<1, h k,n-t is zero, v f,n is an element in V, represents a Poisson distribution, and Γ(·) represents the Gamma distribution, w(t) f,k represents an element in W; S3. Based on the generation model obtained in step S2, an objective function is constructed and converted into the following equation; p(V|W)=∫ C,H p(V,C,H|W)dCdH Where p(V|W) is the probability of V given W, and p(V,C,H|W) is the joint probability distribution given W; S4. Use Jensen's inequality to derive the lower bound of the equation obtained in step S3 to obtain the derived lower bound. The expression of the lower bound is: Q=∫ C,H q(C,H)logp(V,C,H|W)dCdH Where q(C,H) is the tool distribution that is easy to handle, and logp(V,C,H|W) represents the logarithmic distribution to be solved; S5. Simplify the lower bound expression obtained in step S4 and approximate it using a decomposable and tractable variational distribution of the following form: In the formula, q(C f,n ) is the probability p t,k,f,n The multinomial distribution of q(h k,n ) is the Gamma distribution, C f,n Represents a matrix of size T×K; S6. According to the principle of variational Bayesian EM algorithm, the simplified parameter factor q(h k,n ) and q(C f,n ) is calculated and q(h k,n ) and q(C f,n ) is updated with the optimal variational distribution, the update formula is as follows: logq(h k,n ) logq(C f,n ) In the formula, ∝ means that the two are equal under the condition of ignoring the constant term; the symbol <·> represents the expectation under a given probability distribution condition, c t,k,f,n Represents the elements in the tensor C, w(t) f,k and h k,n denote the elements in W and H respectively, Γ(·) denotes the Gamma distribution, α and β are the parameters in the Gamma distribution, and α>0, β>0, When n+t>N, <c t,k,f,n+t > is zero, corresponding to the left shift operation of the matrix, h k,n-t Refers to the elements in the matrix obtained by shifting the matrix H right by t units along the dimension n; And satisfy the optimal variational distribution q(h k,n ) has the following sufficient statistics: In the formula, and They represent the variational distribution q(h) that follows the Gamma distribution. k,n )’s shape and scale parameters; c t,k,f,n Represents the elements in tensor C, c t,k,f,n+t represents the element obtained by shifting the tensor C left by t units along the dimension of n, w(t) f,k represents an element in W; S7. Based on the Bayesian rule, the lower bound of the equation obtained in step S3 is derived to obtain the lower bound of the derived objective function. The expression of the lower bound of the objective function is: Q∝∫ C,H q(C,H)log(p(C|W,H)p(H))pdCdH =<logp(C|W,H)> q(C,H) +<logp(H)> q(C,H) In the formula, Q represents the lower bound of the objective function.<logp(C|W,H)> q(C,H) represents the expectation of logp(C|W,H) under a given probability distribution q(C,H),<logp(H)> q(C,H) represents the expectation of logp(H) under a given probability distribution q(C,H), w(t) f,k and h k,n Represent the elements in W and H respectively, h k,n-t It means that the matrix H is the element in the matrix obtained by shifting t units right along the dimension of n; S8. According to the expression of the lower bound of the objective function obtained in step S7, an update criterion in element W is obtained. The update criterion is as follows: In the formula, the symbol <·> represents the expectation under a given probability distribution, w(t) f,k represents the elements in W, c t,k,f,n Represents the elements in the tensor C, h k,n represents the elements in H, h k,n-t Refers to the elements in the matrix obtained by shifting the matrix H right by t units along the dimension n; S9, the hyperparameter α in the sufficient statistics obtained in step S6 k and β k The adjustment is transformed into a β that follows the inverse Gamma distribution k The adjustment of the hyperparameters a and b of the shape and scale in the middle is based on the principle of the variational Bayesian EM algorithm. The update formulas of the hyperparameters a and b are inferred to obtain the update formulas of the hyperparameters a and b. The update formulas of the hyperparameters a and b are as follows: Where a and b are the parameters β related to the control obeying the inverse Gamma distribution k , K is the number of basis functions in the dictionary W, and ψ is defined as The digamma function, and They are respectively subject to Gamma distribution q(h k,n )’s shape and scale parameters; S10. Automatically update the hyperparameters a and b according to the update formula obtained in step S9. The dictionary W can be learned according to the updated hyperparameters a and b, the update formula obtained in step S6, and the update criterion obtained in step S8.

2. The method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation according to claim 1, characterized in that: Step S1 comprises the following steps: S11, combining the simultaneously occurring activation functions into a basis function; S12. Use the Gamma prior with the same hyperparameters to connect the activation function with the basis function obtained in step S11. The expressions of the activation function and the basis function are as follows: In the formula, is a Gamma distribution, and x≥0, and parameters α>0, β>0, w(t) f,k and h k,n Represent the elements in W and H respectively.

3. The method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation according to claim 1, characterized in that: Step S3 includes the following steps: S31, taking the likelihood function of the generative model obtained in step S2 as the objective function, and dividing the activation function by integration, to obtain the expression of marginal likelihood, as shown below: p(V|W)=∫ H p(V,H|W)dH Where p(V|W,H) is the likelihood function of the generative model, and p(H) is the prior of the matrix H; S32, using a tensor C with four dimensions as an intermediate variable to be introduced into the generative model obtained in step S2, segmenting the observed data V, thereby decoupling the coupling relationship between different basis functions and time slices, and converting the generative model obtained in step S2 into the following combination form: In the formula, c t,k,f,n is an implicit variable of the four-dimensional tensor C, w(t) f,k represents the elements in W, h k,n-t Refers to the elements in the matrix obtained by shifting the matrix H right by t units along the dimension n; S33. Based on step S32, the marginal likelihood expression obtained in step S31 is converted into the following equation: p(V|W)=∫ C,H p(V,C,H|W)dCdH Where C is an intermediate variable with four dimensions, and p(V,C,H|W) is the joint probability distribution under the given W condition.

4. The method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation according to claim 1, characterized in that: Step S6 comprises the following steps: S61, the simplified parameter factor q(h k,n ) and q(C f,n ) meets the following expectations: In the formula, <·> π It represents the expectation under the probability distribution π, symbol A -i,j It means that the elements in A are i,j The remaining elements, logp(V,C,H|W) are the logarithmic distributions to be solved; S62, expand logp(V,C,H|W) in step S61 to obtain the expanded formula: logp(V,C,H|W) =logp(V|C)+logp(C|W,H)+logp(H) Where δ(·) represents the Kronecker delta function; S63, according to the expected expression of step S61 and the expansion obtained in step S62, derive q(h k,n ) and q(C f,n ) is updated with the optimal variational distribution, the update formula is as follows: logq(h k,n ) logq(C f,n ) In the formula, ∝ means that the two are equal under the condition of ignoring the constant term; the symbol <·> represents the expectation under a given probability distribution condition, c t,k,f,n Represents the elements in the tensor C, w(t) f,k and h k,n denote the elements in W and H respectively, Γ(·) denotes the Gamma distribution, α and β are the parameters in the Gamma distribution, and α>0, β>0, When n+t>N, <c t,k,f,n+t > is zero, corresponding to the left shift operation of the matrix, h k,n-t It is the element in the matrix obtained by shifting the matrix H right by t units along the dimension n. And satisfy the optimal variational distribution q(h k,n ) has the following sufficient statistics: In the formula, and They represent the variational distribution q(h) that follows the Gamma distribution. k,n )’s shape and scale parameters; c t,k,f,n Represents the elements in tensor C, c t,k,f,n+t represents the element obtained by shifting the tensor C left by t units along the dimension of n, w(t) f,k Represents an element in W.

5. The method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation according to claim 1, characterized in that: Step S8 comprises the following steps: S81. According to the expression of the lower boundary of the objective function obtained in step S7, the partial derivative of the lower boundary with respect to the dictionary W is calculated to obtain the following relationship: Where w(t) f,k represents the elements in W, c t,k,f,n Represents the elements in the tensor C, h k,n-t It is the element in the matrix obtained by shifting the matrix H right by t units along the dimension n; S82. According to the relationship obtained in step S81, an update criterion in the element W is obtained. The update criterion is as follows: In the formula, the symbol <·> represents the expectation under a given probability distribution, w(t) f,k represents the elements in W, c t,k,f,n Represents the elements in the tensor C, h k,n-t It is the element in the matrix obtained by shifting the matrix H right by t units along the dimension n.

6. The method for learning a minimalist convolutional non-negative dictionary from a spectrogram based on maximum marginal likelihood estimation according to claim 1, characterized in that: Step S9 comprises the following steps: S91, follow the ARD NMF prior and set α k =1, converting the Gamma distribution into an exponential distribution, β k Adding the inverse Gamma distribution prior will adjust the hyperparameter α in the sufficient statistics obtained in step S6 k and β k The adjustment of is converted into the adjustment of the hyperparameters a and b of the shape and scale in the inverse Gamma distribution; S92, estimate the hyperparameters a and b in step S91 to obtain the likelihood function p(V|a,b), which is expressed as follows: p(V|a,b) =∫ C,H,W,β p(V,C,H,W,β|a,b)dCdHdWdβ Where p(V|a,b) is the likelihood function of V given the parameters a and b, which is obtained by integrating and dividing the parameters C, H, W and in the likelihood function p(V,C,H,W,β|a,b); S93. Based on the likelihood function obtained in step S92, the lower bound of the log-likelihood function is derived as a new objective function by using Jensen's inequality. The expression of the new objective function is: logp(V|a,b) ≥∫ C,H,W,β q(C,H,W,β) ∝∫ C,H,W,β q(C,H,W,β) logp(V,C,H,W,β|a,b)dCdHdWdβ Where q(C,H,W,β) is the variational distribution; S94. According to the principle of variational Bayesian EM algorithm, the variational distribution q(C, H, W, β) in step S93 is transformed into the following decomposable form: In the formula, q(C f,n ) is the matrix of random variables following multinomial distribution, q(w(t) f,k ) and q(h k,n ) follow the Gamma distribution and share the same parameter β k ,q(β k ) obeys the inverse Gamma distribution; S95. Based on step S94, the optimal variational distribution is obtained, and then the lower boundary is calculated. The update formulas of the hyperparameters a and b are inferred according to the lower boundary to obtain the update formulas of the hyperparameters a and b. The update formulas of the hyperparameters a and b are as follows: Where a and b are the parameters β related to the control obeying the inverse Gamma distribution k , K is the number of basis functions in the dictionary W, and ψ is defined as The digamma function, and They respectively represent the inverse Gamma distribution q(h k,n )’s shape and scale parameters.

Images