Signal estimation device and program

Abstract

The original signal is estimated with high accuracy using a sample signal based on the observed signal and information representing its missing values. [Solution] The estimated original signal is optimized according to an objective function based on the estimated original signal, a sample signal based on observation signals obtained by one or more sensors, missing information representing missing values ​​in the sample signal, a first graph representing the position of the sensors and the relationship between the sensors, a second graph representing the discrete frequencies of the sample signal and the relationship between the discrete frequencies, and a third graph representing the relationship between discrete times and adjacent discrete times.

Description

[Technical Field]

[0001] This invention relates to a technique for estimating an original signal from an observed signal, and more particularly to a technique using graph signal processing. [Background technology]

[0002] Assume a scenario where one or more sensors acquire observation signals, and sample signals are acquired at regular time intervals based on these observation signals. However, not all observation signals are necessarily obtained, and the sample signals have gaps. In this situation, the original signal is estimated using the sample signals and information representing the gaps. However, this is an underdetermined problem, so other conditions need to be introduced.

[0003] Non-patent document 1 discloses a technique for transforming an underdetermined problem into a decision problem by introducing a sensor graph representing the position of sensors and the relationships between sensors, and a known time graph (path graph) representing the relationships between discrete times and adjacent discrete times. Since the time graph is known, the original signal and the sensor graph become the targets of estimation. [Prior art documents] [Non-patent literature]

[0004] [Non-Patent Document 1] X. Jiang, Z. Tian, ​​and K. Li, “A graph-based approach for missing sensor data imputation,” IEEE Sens. J., vol. 21, no. 20, pp. 23133-23144, Oct. 2021. [Overview of the project] [Problems that the invention aims to solve]

[0005] However, simply introducing sensor graphs and time graphs may not be sufficient to estimate the original signal with high accuracy.

[0006] In this invention, the original signal is estimated with high accuracy using a sample signal based on the observed signal and information representing its missing values. [Means for solving the problem]

[0007] In order to solve the above problems, the present invention optimizes the estimated original signal according to an objective function based on the estimated original signal, a sample signal based on observation signals obtained by one or more sensors, missing information representing missing values ​​in the sample signal, a first graph representing the position of the sensors and the relationship between the sensors, a second graph representing the discrete frequencies of the sample signal and the relationship between the discrete frequencies, and a third graph representing the relationship between discrete times and adjacent discrete times. [Effects of the Invention]

[0008] This makes it possible to estimate the original signal with high accuracy using the sample signal based on the observed signal and information representing its missing values. [Brief explanation of the drawing]

[0009] [Figure 1] Figure 1 is a block diagram illustrating the configuration of the signal estimation system according to the embodiment. [Figure 2] Figure 2 is a flowchart illustrating the signal estimation method of the embodiment. [Figure 3] Figure 3 is a conceptual diagram illustrating the product graph in the embodiment. [Figure 4] Figure 4 is a block diagram illustrating the hardware configuration of the embodiment. [Modes for carrying out the invention]

[0010] Embodiments of the present invention will be described below with reference to the drawings. [Definition] First, let's define the terms. N:N represents the number of sensors. Examples of sensors include microphones, ultrasonic sensors, and biosignal sensors. N is an integer greater than or equal to 1; for example, N is an integer greater than or equal to 2. T:T represents the number of time intervals. Examples of time intervals include time frames and time windows. T is an integer greater than or equal to 1; for example, T is an integer greater than or equal to 2. F: F represents the number of discrete frequencies. Examples of discrete frequencies include frequency bins. F is an integer greater than or equal to 1; for example, F is an integer greater than or equal to 2.

[0011] R: R represents the set of real numbers. C:C represents the set of complex numbers. ∈:η∈Γ indicates that η is an element of Γ. T:η T This represents the transpose of η. ◎:◎ represents the Hadamard product operator. (+):(+) represents the Cartesian Product operator. (×):(×) represents the Kronecker product operator. [×]:[×] represents the Strong Product operator. Scalar product: The product of a scalar value ζ1 and a scalar value ζ2 is denoted as ζ1ζ2. 1 Q :1 Q represents a Q-dimensional column vector where all elements are 1. Q is an integer greater than or equal to 1. 1 T :1 T This represents a row vector where all elements are 1. I Q :I Q This represents the identity matrix of a Q×Q matrix. diag(χ): diag(χ) represents the operation of converting χ into a diagonal matrix. When χ is a square matrix, diag(χ) represents the operation of converting the square matrix χ into a diagonal matrix where the components on the main diagonal of χ (diagonal components) are the components on the main diagonal and the other components are 0. When χ is a vector, diag(χ) represents the operation of converting the vector χ into a diagonal matrix where χ is the component on the main diagonal and the other components are 0. log(η): log(η) represents a matrix whose elements are the natural logarithms of the elements of matrix η. The size of log(η) is the same as the size of matrix η. tr(χ): tr(χ) represents the trace of the square matrix χ (the sum of the elements on the main diagonal of the square matrix χ). vec(η): vec(η) represents a column vector obtained by arranging all the elements of matrix η or tensor η in a column according to a predefined rule. For example, vec(η) may be a column vector obtained by directly flattening all the elements of η (row-first or column-first), or a column vector obtained by expanding η into a matrix according to a predefined rule and then arranging the matrix in a column. ||η||: ||η||1 represents the norm of η. ||η||1: ||η||1 represents the L1 norm of η. ||η|| 1,off : ||η|| 1,off represents the L1 norm of the off-diagonal elements of η. ||η||2: ||η||2 represents the L2 norm of η. ||η|| * : ||η|| * represents the nuclear norm of η. For the nuclear norm, see, for example, Reference 1, etc. Reference 1: C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan, “Tensor Robust Principal Component Analysis with a new tensor nuclear norm,” IEEE Trans. Pattern Anal. Mach. Intel., vol. 42, no. 4, pp. 925-938, Apr. 2020. Note that in Reference 1, the Fast Fourier Transform is performed on the third axis and matrix formation is performed on the first and second axes. However, it is also possible to perform the Fast Fourier Transform on the second axis and matrix formation on the third and first axes, or to perform the Fast Fourier Transform on the first axis and matrix formation on the second and third axes.

[0012] Graph G=(ν, ε): A graph G = (ν, ε) represents a pair of vertices ν and edges ε that connect them. Vertex ν is a fundamental element that constitutes a graph G, and is also called a node or point. Vertex ν represents information about an individual object. Information about an individual object includes, for example, information about the object itself, its position, its state, its properties, and its performance. In this form of graph G, the vertices are, for example, spatial positions, discrete frequencies, discrete times, etc. An edge ε is a line segment connecting two vertices ν, representing the relationship between those vertices ν. If two vertices ν are related, they are connected by an edge ε; if they are not related, they are not connected by an edge ε. Graphs G can be classified into (1) undirected graphs, where edges have no direction and represent undirected relationships between vertices; (2) directed graphs (or digraphs), where edges have direction and represent unidirectional relationships from one vertex to another; (3) weighted graphs, where weights representing the relationships between vertices are assigned to the edges; and (4) unweighted graphs, where weights representing the relationships between vertices are not assigned to the edges.

[0013] Adjacency Matrix W: The adjacency matrix W of graph G=(ν, ε) is the number of vertices v in graph G. i ∈ν and vertex v j ∈ν is edge e i,j w is a value that indicates whether or not they are connected by ∈ε. i,j ∈R is a Θ×Θ matrix W∈R whose (i,j) components are ∈R. Θ×Θ Here, Θ represents the number of vertices ν in graph G. Θ is an integer greater than or equal to 1, for example, Θ is an integer greater than or equal to 2. Also, i=1,...,Θ and j=1,...,Θ. If graph G is an undirected graph, then vertex v i ∈ν and vertex v j ∈ν is edge e i,j If they are connected by ∈ε, then (e i,j (If ∈ε then) w i,j >0, otherwise w i,j = 0. If graph G is a directed graph, then vertex v i From ∈ν to vertex v j Edge e toward ∈ν i,j If they are connected by ∈ε, then w i,j >0, otherwise w i,j = 0. If graph G is a weighted graph, w i,j is vertex v i ∈ν and vertex v j This represents the weight corresponding to the strength of the relationship with ∈ν. If the graph G is an unweighted graph, w i,j is a constant (for example, w i,j =1). For example, if Θ=4 and graph G is an undirected and unweighted graph, an example of the adjacency matrix W is given by equation (1) below.

number

[0014] Degree Matrix: The order matrix D is the matrix of the vertices v of the graph G. i Degree d of ∈ν i Θ×Θ matrix W∈R whose (i,i) component (diagonal element) is ∈R Θ×Θ Therefore, if graph G is an undirected graph, the degree d iis vertex v i This represents the number of edges ε connected to a given point. If the graph G is a directed graph, then the degree d is also represented. i The vertex is v i From ∈ν to vertex v j This can represent the number of edges ε connected to ∈ν, or the degree d. i The vertex is v j From ∈ν to vertex v i This can be divided into two cases: one representing the number of edges ε connected to ∈ν, and the other representing the degree d of the former. i A matrix D whose elements are d is called an out-degree matrix, and the latter degree d i A matrix D whose elements are 1 is called an in-degree matrix. The matrix D has an adjacency matrix W and a column vector 1. Θ Using D=diag(W1 Θ ) can be expressed as follows. For example, in the case of the adjacency matrix W in equation (1), the order matrix D is as shown in equation (2) below.

number

[0015] Laplacian matrix: The Laplacian matrix is ​​defined as L = DW ∈ R for an adjacency matrix W and a degree matrix D. Θ×Θ This is a Θ×Θ matrix that satisfies the following condition. For the adjacency matrix W in equation (1) and the order matrix D in equation (2), the Laplacian matrix L is given by equation (3) below.

number

[0016] Product Graphs: A product graph is a graph constructed by combining multiple graphs. In other words, a product graph is a graph obtained by multiplying graphs. Representative product methods for graphs include the direct product, the Kronecker product, and the strong product. A product graph obtained by the direct product is called a direct product graph, a product graph obtained by the Kronecker product is called a Kronecker product graph, and a product graph obtained by the strong product is called a strong product graph.

[0017] Cartesian product graph: Cartesian product of graph G1 and graph G2: graph G c =G1(+)G2 is a graph that satisfies the following relationships (4) and (5). W (+) =W1(+)W2(4) L (+) =L1(+)L2(5) Here, W (+) This is a Cartesian product graph G c W1 represents the adjacency matrix of graph G1, W2 represents the adjacency matrix of graph G2, and L (+) This is a Cartesian product graph G c L1 represents the Laplacian matrix of graph G1, and L2 represents the Laplacian matrix of graph G2.

[0018] Kronecker product graph: The Kronecker product of graph G1 and graph G2 is graph G k =G1(×)G2 is a graph that satisfies the following relationships (6) and (7). W (×) =W1(×)W2(6) L (×) =D1(×)D2-W1(×)W2(7) Here, W (×) This is the Kronecker product graph G c W1 represents the adjacency matrix of graph G1, W2 represents the adjacency matrix of graph G2, and L (×) This is the Kronecker product graph G k D1 represents the Laplacian matrix of graph G1, and D2 represents the order matrix of graph G2.

[0019] Square graph: The strong product of graph G1 and graph G2 is graph G s =G1[×]G2 is a graph that satisfies the following relationships (8) and (9). W [×] =W (×) +W (+) (8) L [×] =L (×) +L(+) (9) Here, W [×] This is a strongly multiplicative graph G s This represents the adjacency matrix of L [×] This is a strongly multiplicative graph G s This represents the Laplacian matrix.

[0020] [principle] Next, the principle of this embodiment will be explained. As mentioned above, Non-Patent Document 1 discloses a technique that transforms an underdetermined problem into a decision problem by introducing a sensor graph representing the position of sensors and the relationships between sensors, and a time graph representing the relationships between discrete times and known adjacent discrete times. This makes it possible to estimate the original signal while considering the positional relationships of the sensors and the time. However, Non-Patent Document 1 does not consider the periodicity of the signal. That is, if the original signal has periodicity, there is a relationship between observed signals at time intervals that are integer multiples of the period. For example, acoustic signals, ultrasonic signals, biological signals, seismic waveform signals, electrical signals, optical signals, mechanical vibration signals, electromagnetic wave signals, etc., are periodic signals, and there is a relationship between observed signals that are integer multiples of the period. However, the method in Non-Patent Document 1 cannot utilize this feature. In this embodiment, the estimation accuracy of the original signal is improved by further introducing a frequency graph that represents this feature.

[0021] <Question setting> A noise signal is superimposed on the original signal emitted from one or more signal sources, and N sensors 121, ..., 12 N Assume the situation observed by sensors 121, ..., 12 N The observed signals are digitized, and sample signals based on the observed signals are obtained at regular time intervals. However, not all observed signals are obtained, and there are gaps in the sample signals. Here, sensor 12 n Based on the observed signals (n=1,...,N), discrete time 14 t The sample signal obtained at (t=1,...,T) z n,t Let ∈R represent the sample signal z n,t An N×T matrix with (n,t) elements is used for the sample signal Z∈R. N×Tis expressed as. Also, the original signal corresponding to the sample signal z n,t is expressed as x n,t ∈R, and the original signal x n,t is expressed as an estimated original signal X ∈ R of an N×T matrix with (n, t) elements N×T . Further, the noise signal corresponding to the sample signal z n,t is expressed as e n,t ∈R, and the noise signal e n,t is expressed as an estimated noise signal E ∈ R of an N×T matrix with (n, t) elements N×T . Also, the mask representing the loss of the sample signal z n,t is expressed as m n,t ∈R, and the mask m n,t is expressed as a mask (loss information) M ∈ R of an N×T matrix with (n, t) elements N×T . Note that the sample signal z n,t , the original signal x n,t , the noise signal e n,t , and the mask m n,t are values in the time domain. Here, the sample signal Z can be represented by the following observation model shown in Equation (9a) using the mask M, the estimated original signal X, and the estimated noise signal E Z = M ◎ (X + E) (9a) In this embodiment, it is an object to obtain the estimated original signal X based on the determined mask M and the acquired sample signal Z

[0022] <Introduction of Graphs> However, as it is in Equation (9a), it is an ill - determined problem. Therefore, in this embodiment, a sensor graph G1 (first graph), a frequency graph G2 (second graph), and a time graph G3 (third graph) are introduced The sensor graph G1 = (ν1, ε1) is a graph representing the positions of sensors 121,..., 12 N and the relationships between sensors 121,..., 12 N . That is, the vertex ν1 of the sensor graph G1 represents the positions of sensors 121,..., 12 N , and the edge ε1 represents the relationships between sensors 121,..., 12 N . The adjacency matrix of the sensor graph G1 is expressed as W1 ∈ R N×N , and the degree matrix is expressed as D1 ∈ R N×NThis is expressed as follows, and the Laplacian matrix is ​​L1∈R N×N This is how it is expressed. The frequency graph G2=(ν2, ε2) shows the discrete frequencies of the sample signal Z, 131, ..., 13 F and discrete frequencies 131, ..., 13 F This is a graph that shows the relationship between them. That is, the vertex ν2 of the frequency graph G2 is the discrete frequency 131, ..., 13 F This represents the discrete frequencies 131, ..., 13 F This represents the relationship between elements. The adjacency matrix of the frequency graph G2 is W2∈R F×F This is expressed as follows, and the order matrix D2∈R F×F This is expressed as follows, and the Laplacian matrix is ​​L2∈R F×F This is how it is expressed. The time graph G3=(ν3, ε3) is a discrete time graph of 141, ..., 14 F and adjacent discrete time 14 t ,14 t+1 This is a graph that shows the relationship between two points. Here, t=1,...,T-1. That is, the vertices ν3 of the time graph G3 are discrete time points 141,...,14 F This represents the case where edge ε3 is adjacent to discrete time 14 t ,14 t+1 This represents the relationship between elements. The time graph G3 is a known path graph (see, for example, Non-Patent Document 1). The adjacency matrix of the time graph G3 is W3∈R. T×T This is expressed as follows, and the order matrix D3∈R T×T This is expressed as follows, and the Laplacian matrix is ​​L3∈R T×T This is how it is expressed. The sensor graph G1 (first graph), frequency graph G2 (second graph), and time graph G3 (third graph) may be undirected graphs, directed graphs, weighted graphs, or unweighted graphs.

[0023] <This methodology> This signal estimation device uses an estimated source signal X and one or more sensors 121, ..., 12 N Sample signal Z based on the observed signal obtained by, mask (missing information) M representing missing values ​​in sample signal Z, sensors 121, ..., 12 NPosition and sensors 121, ..., 12 N Sensor graph G1 (first graph) showing the relationship between the two, discrete frequencies of sample signal Z 131, ..., 13 F and discrete frequencies 131, ..., 13 F Frequency graph G2 (second graph) showing the relationship between the two, and discrete time 141, ..., 14 F and adjacent discrete time 14 t ,14 t+1 The estimated original signal X is determined according to an objective function S based on a time graph G3 (third graph) representing the relationship (t=1,...,T-1). The signal estimation device uses, for example, information representing the sample signal Z to determine the estimated original signal X according to the objective function S. Here, since the time graph G3 is known, the targets for estimation are the fixed original signal X, the sensor graph G1, and the frequency graph G2. For example, the objective function S is a convex function with respect to each of the fixed original signal X, the sensor graph G1, and the frequency graph G2, and a non-convex function with respect to the set of fixed original signal X, the sensor graph G1, and the frequency graph G2.

[0024] Preferably, the signal estimation device of this embodiment performs the following: an estimated original signal optimization process that optimizes the estimated original signal X while fixing the sensor graph G1 (first graph) and frequency graph G2 (second graph) according to an objective function S; a first graph optimization process that optimizes the sensor graph G1 (first graph) while fixing the estimated original signal X and frequency graph G2 (second graph) according to an objective function S; and a second graph optimization process that optimizes the frequency graph G2 (second graph) while fixing the estimated original signal X and sensor graph G1 (first graph) according to an objective function S. For example, the signal estimation device of this embodiment repeats the estimated original signal optimization process, the first graph optimization process, and the second graph optimization process until a predetermined termination condition is met. This makes it possible to estimate the estimated original signal X even when it is not possible to optimize the estimated original signal X, sensor graph G1 (first graph), and frequency graph G2 (second graph) simultaneously.

[0025] Preferably, the objective function S in this embodiment is a function that represents the fidelity of the estimated original signal X and mask (missing information) M to the sample signal Z, the smoothness of the estimated original signal X with respect to the product graph of the sensor graph G1 (first graph), frequency graph G2 (second graph), and time graph G3 (third graph), the low rank of the estimated original signal X, the connectivity of the sensor graph G1 (first graph), the connectivity of the frequency graph G2 (second graph), the sparsity of the sensor graph G1 (first graph), and the sparsity of the frequency graph G2 (second graph). Here, the fidelity of the estimated original signal X and mask M to the sample signal Z represents the degree to which the observation model in equation (1) based on the estimated original signal X and mask M matches the sample signal Z. The higher the fidelity, the better the observation model matches the sample signal Z. The smoothness of the estimated original signal X with respect to the product graph represents the tendency for the estimated original signal X to take similar values ​​for vertices that have a strong relationship on the product graph. The higher the smoothness, the stronger this tendency. The low-rank nature of the estimated original signal X indicates whether the estimated original signal X can potentially be represented with a low rank. The lower the rank of the estimated original signal X, the lower the rank it can be represented with. The connectivity of a graph indicates how much the graph is grouped into a single connected component. The higher the connectivity of a graph, the more likely it is to be grouped into a single connected component. The sparsity of a graph indicates how little the graph has unnecessary edges. The higher the sparsity of a graph, the fewer unnecessary edges it has. For more information on these meanings, please refer to, for example, Non-Patent Document 1. Note that the product graph may be, for example, a Cartesian product graph, a Kronecker product graph, or a strong product graph.

[0026] For example, the objective function S in this embodiment includes the sum of terms that decrease as the fidelity of the estimated original signal X and mask M to the sample signal Z increases, terms that decrease as the smoothness of the estimated original signal X with respect to the product graph increases, terms that decrease as the estimated original signal X can be represented at a lower rank, terms that decrease as the connectivity of the sensor graph G1 (first graph) increases, terms that decrease as the connectivity of the frequency graph G2 (second graph) increases, terms that decrease as the sparsity of the sensor graph G1 (first graph) increases, and terms that decrease as the sparsity of the frequency graph G2 (second graph) increases. For example, the estimated original signal optimization process in this embodiment obtains information representing the estimated original signal X that minimizes the objective function S while the sensor graph G1 (first graph) and frequency graph G2 (second graph) are fixed; the first graph optimization process obtains information representing the sensor graph G1 (first graph) that minimizes the objective function S while the estimated original signal X and frequency graph G2 (second graph) are fixed; and the second graph optimization process obtains information representing the frequency graph G2 (second graph) that minimizes the objective function S while the estimated original signal X and sensor graph G1 (first graph) are fixed.

[0027] Equation (10) below shows a specific example of the objective function S. S=(1 / 2)||M◎(XZ)|| 2 +(α / 2)y T Ly+μ||Y|| * +Σ i=1,2 (-β1 (i) tr(log(diag(L i )))+β2 (i) ||L i || 1,off ) (10) Here, α,μ,β1 (1) ,β1 (2) ,β2 (1) ,β2 (2) ∈R is a non-negative real constant, for example, a positive real constant. Y represents the estimated original signal in the time-frequency domain, where Y=Fourier(X)∈C N×F×V The following conditions are met: Y=Fourier(X) represents the operation of converting the estimated original signal X in the time domain to the estimated original signal Y in the time-frequency domain. That is, the estimated original signal Y ∈ C in the time-frequency domain.N×F×T The estimated original signal X∈R N×T The original signal x is an element of n,1 ,...,x n,T ∈R in the time interval 14' τ The original signal y in the time-frequency domain is obtained by converting each step to the time-frequency domain. n,f,τ This is a third-order tensor whose elements are . Here, n=1,...,N is sensor 12 n This represents f=1,...,F is a discrete frequency 13 f This represents τ=1,...,V is in the time interval 14' τ This represents the signal. V is an integer greater than or equal to 1, and V may or may not be T. For the conversion from the time domain to the time-frequency signal, for example, the Short-Time Fourier Transform can be used. y also represents the estimated original signal in the time-frequency domain, and y = vec(Y) ∈ C NFT The following conditions are met. L is the Laplacian matrix of the product graph of the sensor graph G1, the frequency graph G2, and the time graph G3. If the product graph is a Cartesian product, then L = L1(+)L2(+)L3 (see Figure 3). If the product graph is a Kronecker product, then L = L1(×)L2(×)L3. If the product graph is a strong product, then L = L1[×]L2[×]L3. tr(log(diag(L i )) is log(diag(L i Since it is the trace of the Laplacian matrix L i This represents the sum of the diagonal elements. Here, (1 / 2)||M◎(XZ)|| 2 This term becomes smaller as the fidelity of the estimated original signal X and the mask M to the sampled signal Z increases. (α / 2)y T Ly is a term that becomes smaller as the smoothness of the estimated original signal X increases with respect to the product graph of the sensor graph G1, frequency graph G2, and time graph G3. μ||Y|| * This is the term that becomes small enough for the estimated original signal X to be represented by a low rank. -β1 (1) tr(log(diag(L1)) is a term that becomes smaller as the connectivity of the sensor graph G1 increases, and -β1 (2) tr(log(diag(L2))) is a term that becomes smaller as the connectivity of the frequency graph G2 increases. β2 (1)||L1|| 1,off This term becomes smaller as the sparsity of the sensor graph G1 increases, and β2 (2) ||L2|| 1,off This term becomes smaller as the sparsity of the frequency graph G2 increases. The signal estimation device of this form, for example, finds the solution Y,L1,L2 to the optimization problem of equation (11) below.

number

[0028] <Solution to the optimization problem in equation (11)> An example of how to solve the optimization problem in equation (11) is given. First, the objective function S in equation (10) is (1 / 2)||M◎(XZ)|| 2 Mask M∈R N×T The sample signal z n,t The missing part of sensor 12 n We transform this into a time-frequency domain mask based solely on the absence of (n∈{1,...,N}). That is, a time-domain mask M∈R N×T The time-frequency mask M shown in equation (12) below tf ∈{0,1} N×FT It transforms into this. M tf =m(×)1 FT T ∈{0,1} N×FT (m∈{0,1} N ) (12) Here, the mask M∈R N×T The time domain mask m included n,1 ,...,m n,T time interval 14' τ The mask obtained by converting each step to the time-frequency domain is m n,f,τ This is expressed as, mask m n,f,τ An N×FT matrix with elements M tf ∈R N×FV This is how it is expressed. n=1,...,N is sensor 12 n This represents f=1,...,F is a discrete frequency 13 f This represents τ=1,...,V is in the time interval 14' τ This represents the conversion from the time domain to a time-frequency signal, for example, the short-time Fourier transform can be used. Also, the objective function S in equation (10) is (1 / 2)||M◎(XZ)|| 2 Estimated original signal X∈R N×T and sample signal Z∈R N×T The estimated original signal X in the time-frequency domain tf ∈{0,1} N×FT and sample signal Z tf ∈R N×FV It is transformed into this. Here, the estimated original signal X∈R N×T The original signal x in the time domain included n,1 ,...,x n,T ∈R in the time interval 14' τThe estimated original signal in the time-frequency domain obtained by converting each to the time-frequency domain is x n,f,τ This is expressed as the estimated original signal x in the time-frequency domain. n,f,τ An N×FT matrix with elements is used to estimate the original signal X in the time-frequency domain. tf ∈R N×FV Let's assume that the sample signal Z∈R N×T The time-domain sample signal z included n,1 ,...,z n,T ∈R in the time interval 14' τ The time-frequency domain sample signal obtained by converting each one to the time-frequency domain is z n,f,τ This is expressed as the sample signal z in the time-frequency domain. n,f,τ An N×FT matrix with elements is used to represent the sample signal Z in the time-frequency domain. tf ∈R N×FV This is how it is expressed. Also, M tf ,X tf ,Z tf The order of elements in M ​​is predetermined. tf ,X tf ,Z tf The same position corresponds to the same (n,f,τ) m n,f,τ ,x n,f,τ ,z n,f,τ The following is arranged. Therefore, the objective function S in equation (10) is (1 / 2)||M◎(XZ)|| 2 This can be transformed into equation (13) below. (1 / 2)||M tf ◎(X tf -Z tf )|| 2 (13) Furthermore, equation (13) can be transformed into equation (14) below. (1 / 2)||M tf ◎(X tf -Z tf )|| 2 =(1 / 2)||M tf ◎X tf -Z tf || 2 =(1 / 2)||vec(M tf ◎X tf )-vec(Z tf )|| 2 (14)

[0029] Here, Fourier(X) processes sensors 121, ..., 12 N Since they are identical in each case, vec(Y) can be expressed as shown in equation (15) below. vec(Y)=(FOURIER(×)I N )vex(X tf ) (15) Here, FOURIER is an N×N matrix. From equation (15), the following equation (16) holds. vec(X tf )=(FOURIER -1 (×)I N )vec(Y) (16) Here, FOURIER -1 This is the inverse matrix of FOURIER. From equations (12) and (16), we have vec(M) in equation (14). tf ◎X tf ) can be transformed as shown in equation (16) below. vec(M tf ◎X tf ) =diag(vec(M tf ))vec(X tf ) =diag(vec(m(×)1 FT T ))(FOURIER -1 (×)I N )vec(Y) =(I FT (×)diag(m))(FOURIER -1 (×)I N )vec(Y) =(FOURIER -1 (×)I N )(I FT (×)diag(m))vec(Y) =(FOURIER -1 (×)I N )diag(vec(M tf ))vec(Y) (17) Therefore, equation (14) can be transformed using equation (17) to become equation (18). (1 / 2)||M tf ◎(Xtf -Z tf )|| 2 =(1 / 2)||vec(M tf ◎X tf )-vec(Z tf )|| 2 =(1 / 2)||(FOURIER -1 (×)I N )diag(vec(M tf ))vec(Y)-vec(Z tf )|| 2 =(1 / 2)||M tens ◎Y-Fourier(Z)|| 2 (18) Here, M tens ∈R N×F×T is vec(M tens )=vec(M tf This is a third-order tensor that satisfies ). Therefore, the objective function S in equation (10) can be transformed into equation (19) below. S=(1 / 2)||M tens ◎Y-Fourier(Z)|| 2 +(α / 2)y T Ly+μ||Y|| * +Σ i=1,2 (-β1 (i) tr(log(diag(L i )))+β2 (i) ||L i || 1,off ) (19) Here, Fourier(Z)∈C N×F×V The sample signal Z∈R N×T The sample signal z is an element of n,1 ,..., z n,T ∈R in the time interval 14' τ The time-frequency domain sample signal z obtained by converting each step to the time-frequency domain n,f,τ This is a third-order tensor whose elements are .

[0030] Therefore, the signal estimation device of this form can be implemented by, for example, repeatedly and alternately optimizing Y, L1, and L2 independently according to equation (11) for the objective function S of equation (19). Thus, the optimization problem of equation (11) is divided into the subproblems of equations (20), (21), and (22) below.

number

number

number

[0031] <Solution to the subproblem of equation (20)> An example of how to solve the subproblem of equation (20) is given. Here, we demonstrate how to solve the subproblem of equation (20) using the Alternating Direction Method of Multipliers (ADMN). Variable U∈C N×F×T and the dual variable P∈C N×F×T By introducing the Lagrange multiplier (real number) ρ>0, the extended Lagrange function L(Y,U,P) in equation (20) becomes equation (23) below. L(Y,U,P)=(1 / 2)||M tens ◎Y-Fourier(Z)|| 2 +(α / 2)y T Ly+μ||U|| * +(ρ / 2)||Y-U+P / ρ|| 2 (twenty three) In contrast, the subproblem in equation (18) can be solved by repeatedly applying the following update equations (24), (25), and (26) to the variables Y, U, and P.

number

number

[0032] <Solution to the update equation (24)> The optimization of the update equation in equation (24) is the minimization of a convex quadratic function. Therefore, L(Y,U (κ) ,P (κ) The gradient ∇ of Y ) Y L(Y,U (κ) ,P (κ) Regarding the equation ∇ Y L(Y,U (κ) ,P (κ) Solve )=0 and take the solution Y as Y (κ+1)This would suffice. However, the size of the matrix included in this equation is enormous (NFT × NFT), making this method impractical. Therefore, the most reasonable way to solve the update equation in equation (24) is with a gradient-based iterative method, such as the conjugate gradient method. Note that the gradient used in the iterative method is ∇ Y L(Y,U (κ) ,P (κ) ) can be calculated as shown in equation (27) below. ∇ Y L(Y,U (κ) ,P (κ) )=M tf ◎Y-Fourier(Z)+αvec T (Ly)+ρ(Y-U+P / ρ) (27)

[0033] <Solution to the update equation (25)> The optimization of the update equation in equation (25) is a proximity mapping calculation. That is, by calculating equation (28) below, U (κ+1) You just need to calculate this.

number

[0034] <Solution to the subproblem of equation (21)> Searching for a solution L1 ∈ Lap to the subproblem of equation (21) from among matrices Laps that satisfy the constraints of the Laplacian matrix is ​​computationally intensive. Therefore, we will search for an L1 that satisfies the constraint that the upper triangular components of the adjacency matrix are non-negative. To do this, we will transform equation (21) and solve the optimization problem expressed by this transformed equation using PDS (Primal-Dual Splitting). Here, the first and second terms of equation (21) are differentiable, and the third term and the constraint that the upper triangular components of the adjacency matrix are non-negative, when expressed as indicator functions, are proximity-mappable.

[0035] <Term 1 of Equation (21)> The first term of equation (21) is (α / 2)y TLy is differentiable with respect to the adjacency matrix W1 of the sensor graph G1. Therefore, using the adjacency matrix W1, (α / 2)y T Let Ly be expressed. Here, L is the Laplacian matrix of the product graph of the sensor graph G1, frequency graph G2, and time graph G3, but the result will differ depending on the type of product. Below, Y1∈C N×FT is a tensor Y∈C N×F×T This represents the matrix obtained by expanding Y1. - ∈C N×FT is vec(Y1 - This represents a matrix that satisfies ) = diag(y)·y. W2 is the adjacency matrix of the frequency graph G2, and W3 is the adjacency matrix of the time graph G3. C is a constant. (i) If the product graph is a direct product graph and L = L1(+)L2(+)L3, y T Ly=tr(S1 (×) W1) However, S1 (×) =Y1 - (W2(×)W3)1 FT 1 N T -Y1(W2(×)W3)Y1 T That is the case. (ii) If the product graph is a Kronecker product graph and L = L1(×)L2(×)L3, y T Ly=tr(S1 (+) W1) + C. However, S1 (+) =Y1 - 1 FT 1 N T -Y1Y1 T That is the case. (iii) If the product graph is a strong product graph and L = L1[×]L2[×]L3, y T Ly=tr(S1 [×] W1) + C. However, S1 [×] =S1 (×) +S1 (+) That is the case. From now on, we will not distinguish between products. T Let Ly = tr(S1W1) + C. Next, define the replication matrix B1 of equation (29) below.

number

number

[0036] <Terms 2 and 3 of Equation (21)> The second term of equation (21) can be transformed into equation (31) below. -β1 (1) tr(log(diag(L1)))=-β1 (1) 1 T log(R1w1) (31) Here,

number

[0037] <Transformation and solution of the subproblem (21)> From equations (30), (31), and (32) above, the subproblem of equation (21) can be transformed into equation (33) below.

number

number

number

number

[0038] <Solution to the subproblem in equation (22)> Similar to equation (21), the subproblem of equation (22) can be transformed into equation (37) below.

number

number

number

number

number

[0039] [First Embodiment] Next, a first embodiment of the present invention will be described. <Structure> The configuration of the signal estimation system 1 in this embodiment is illustrated below. As illustrated in Figure 1, the signal estimation system 1 of this embodiment includes a signal estimation device 11 and a sensor group 12. The signal estimation device 11 includes a memory 110, a signal processing unit 111, a time-frequency domain conversion unit 112, an estimated original signal optimization unit 113, a first graph optimization unit 114, a second graph optimization unit 115, a control unit 116, and a time-domain conversion unit 117. The signal estimation device 11 performs each process under the control of the control unit 116. Information obtained from each part of the signal estimation device 11 is stored in the memory 110 sequentially and used when called upon as needed. The sensor group 12 consists of N sensors 121, ..., 12 N It has sensors 121, ..., 12 N This is a device for observing signals. Examples of signals include acoustic signals, ultrasonic signals, biological signals, seismic waveform signals, electrical signals, optical signals, mechanical vibration signals, electromagnetic wave signals, etc. Sensors 121, ..., 12 N Specific examples include microphones, ultrasonic transducers, piezoelectric sensors, piezoresistive sensors, electrocardiogram sensors, electroencephalogram sensors, electromyogram sensors, seismometers, accelerometers, voltage sensors, and photodiodes.

[0040] <Processing> Next, the processing of the signal estimation device 11 in this embodiment will be explained. A noise signal is superimposed on the original signal emitted from one or more signal sources, and N sensors 121, ..., 12 N Assume the situation observed by sensors 121, ..., 12 N The observed signals are sent one by one to the signal processing unit 111. The signal processing unit 111 processes the sensor 12 n The observed signal (n=1,...,N) is digitized, and the sample signal z n,t Obtain ∈R. Each sensor 12 n Based on the observed signals (n=1,...,N), each discrete time 14 t The sample signal z obtained at (t=1,...,T) n,t ∈R is buffered in the signal processing unit 111, and the sample signal Z∈R N×T This is sent to the time-frequency domain conversion unit 112 (step S111).

[0041] The time-frequency domain conversion unit 112 converts the sample signal Z into the time-frequency domain, and the sample signal Fourier(Z)∈C in the time-frequency domain N×F×V The result is obtained and sent to the estimated original signal optimization unit 113 (step S112).

[0042] The Estimated Original Signal Optimization Unit 113 uses the time-frequency domain sample signal Fourier(Z) (information representing the sample signal Z) to optimize the estimated original signal X according to the objective function S, while fixing the sensor graph G1 and frequency graph G2. For example, according to equations (11) and (19), the Estimated Original Signal Optimization Unit 113 obtains and outputs Y=Y^ (information representing the estimated original signal X) which minimizes the objective function S with L1 and L2 fixed and Y as the variable. The initial values ​​of L1 and L2 are predetermined. When L1^ is sent to the Estimated Original Signal Optimization Unit 113 as described later, L1 is fixed to the latest L1^ sent, and when L2^ is sent to the Estimated Original Signal Optimization Unit 113 as described later, L2 is fixed to the latest L2^ sent. The obtained Y^ (information representing the estimated original signal X) is sent to the first graph optimization unit 114 and the second graph optimization unit 115 (step 113).

[0043] The first graph optimization unit 114 optimizes the sensor graph G1 according to the objective function S, while fixing the estimated original signal X and the frequency graph G2. For example, according to equations (11) and (19), the first graph optimization unit 114 obtains and outputs L1 = L1^(information representing the first graph) which minimizes the objective function S with Y and L2 fixed and L1 as the variable. Note that Y is fixed to the latest Y^ sent from the estimated original signal optimization unit 113. If L2^ is sent to the first graph optimization unit 114 as described later, L2 is fixed to the latest L2^ that was sent. The obtained L1^(information representing the first graph) is sent to the estimated original signal optimization unit 113 and the second graph optimization unit 115 (step 114).

[0044] The second graph optimization unit 115 optimizes the frequency graph G2 according to the objective function S, while fixing the estimated original signal X and the sensor graph G1. For example, according to equations (11) and (19), the second graph optimization unit 115 obtains and outputs L2 = L2^(information representing the second graph) which minimizes the objective function S while fixing Y and L1 and making L2 a variable. Note that Y is fixed to the latest Y^ sent from the estimated original signal optimization unit 113, and L1 is fixed to L1^ sent from the first graph optimization unit 114. The obtained L2^(information representing the second graph) is sent to the estimated original signal optimization unit 113 and the first graph optimization unit 114 (step 115).

[0045] The control unit 116 determines whether or not the termination condition is met (step S116). The termination condition can be anything. Examples of termination conditions include the update width of Y being less than or equal to a predetermined width, or the number of updates reaching a predetermined number. If the termination condition is not met, the process returns to step S113. On the other hand, if the termination condition is met, the estimated original signal optimization unit 113 sends the latest Y^ (information representing the estimated original signal X) as the finally optimized Y~ to the time-domain conversion unit 117, and the process proceeds to step S117.

[0046] When the time-domain conversion unit 117 receives Y~, it converts X~ = Fourier -1The estimated original signal X~ in the time domain is obtained by (Y~) and output (step S118).

[0047] <Features of this form> In this embodiment, the estimated original signal was optimized according to an objective function based on the estimated original signal, a sample signal based on observation signals obtained by one or more sensors, missing information representing missing values ​​in the sample signal, a first graph representing the sensor positions and relationships between sensors, a second graph representing the discrete frequencies of the sample signal and the relationships between discrete frequencies, and a third graph representing the relationships between discrete times and adjacent discrete times. This allows for the acquisition of an estimated original signal that considers not only the relationships between sensors and adjacent discrete times, but also the relationships between discrete frequencies. As a result, the original signal can be estimated with high accuracy using the sample signal based on the observation signal and the information representing its missing values.

[0048] Furthermore, in this embodiment, preferably, the estimated original signal optimization unit 113 optimizes the estimated original signal according to the objective function S while fixing the first graph and the second graph; the first graph optimization unit 114 optimizes the first graph according to the objective function while fixing the estimated original signal and the second graph; and the second graph optimization unit 115 optimizes the second graph according to the objective function while fixing the estimated original signal and the first graph. This makes it possible to estimate the estimated original signal even when it is not possible to optimize the estimated original signal, the first graph, and the second graph simultaneously.

[0049] Furthermore, preferably in this embodiment, the objective function is a function that represents the fidelity of the estimated original signal and the sampled signal of missing information, the smoothness of the estimated original signal with respect to the product graph of the first graph, the second graph, and the third graph, the low rank of the estimated original signal, the connectivity of the first graph, the connectivity of the second graph, the sparsity of the first graph, and the sparsity of the second graph. This makes it possible to estimate an estimated original signal with optimized fidelity, smoothness, low rank, connectivity, and sparsity.

[0050] Furthermore, in this embodiment, the objective function S of equation (10) is preferably (1 / 2)||M◎(XZ)||2 Mask M∈R N×T The sample signal z n,t The missing part of sensor 12 n We transform the mask into a time-frequency domain mask based solely on the absence of (n∈{1,...,N}). That is, we transform it into the mask of equation (12). This reduces the computational complexity.

[0051] [Second Embodiment] The objective function S exemplified in equations (10) and (19) was a function representing the fidelity of the estimated original signal X and the mask (missing information) M to the sample signal Z, the smoothness of the estimated original signal X with respect to the product graph of the sensor graph G1 (first graph), frequency graph G2 (second graph), and time graph G3 (third graph), the low rank of the estimated original signal X, the connectivity of the sensor graph G1 (first graph), the connectivity of the frequency graph G2 (second graph), the sparsity of the sensor graph G1 (first graph), and the sparsity of the frequency graph G2 (second graph). However, the objective function may also be a function representing the fidelity of the estimated original signal and missing information to the sample signal, the smoothness of the estimated original signal with respect to the product graph of the first, second, and third graphs, the low rank of the estimated original signal, the connectivity of the first graph, the connectivity of the second graph, the sparsity of the first graph, and the sparsity of the second graph. Alternatively, the objective function may represent other properties in addition to these. For example, the objective function may represent robustness to outliers and noise, stronger physical and time-series assumptions, physical and practical constraints, etc. Alternatively, the objective function does not have to represent the smoothness of the estimated original signal with respect to the product graph of the first, second, and third graphs, the low-rank nature of the estimated original signal, the connectivity of the first graph, the connectivity of the second graph, the sparsity of the first graph, and at least some of the sparsity of the second graph. For example, the real constants μ,β1 mentioned above. (1) ,β1 (2) ,β2 (1) ,β2 (2) Some of the values ​​may be zero.

[0052] [Hardware configuration] The functions realized by the components described herein may be implemented in a circuitry or processing circuitry, including general-purpose processors, application-specific processors, integrated circuits, ASICs (Application Specific Integrated Circuits), CPUs (a Central Processing Unit), conventional circuits, and / or combinations thereof, programmed to realize the functions described herein. A processor includes transistors and other circuits and is considered a circuitry or processing circuitry. A processor may be a programmed processor that executes a program stored in memory.

[0053] In this specification, circuitry, unit, and means are hardware programmed to perform or execute the functions described herein. Such hardware may be any hardware disclosed herein, or any hardware known to be programmed to perform or execute the functions described herein.

[0054] If the hardware is a processor that is considered to be a type of circuitry, then the circuitry, means, or unit is a combination of hardware and software used to constitute the hardware and / or processor.

[0055] For example, the signal estimation device 11 in each embodiment is a device configured by a general-purpose or dedicated computer, which is equipped with a processor (hardware processor) such as a CPU (central processing unit) and memory such as RAM (random-access memory) and ROM (read-only memory), executing a predetermined program. That is, the signal estimation device 11 in each embodiment has, for example, a processing circuitry configured to implement each of its respective parts. This computer may have one processor and memory, or it may have multiple processors and memories. This program may be installed on the computer, or it may be pre-recorded in ROM, etc. Furthermore, some or all of the processing units may be configured using electronic circuits that realize processing functions independently, rather than electronic circuits that realize the functional configuration by loading a program, such as a CPU. Also, the electronic circuits that constitute one device may include multiple CPUs.

[0056] Figure 4 is a block diagram illustrating the hardware configuration of the signal estimation device 11 in each embodiment. As illustrated in Figure 4, the signal estimation device 11 in this example includes a CPU (Central Processing Unit) 10a, an input unit 10b, an output unit 10c, a RAM (Random Access Memory) 10d, a ROM (Read Only Memory) 10e, an auxiliary storage device 10f, a communication unit 10h, and a bus 10g. The CPU 10a in this example includes a control unit 10aa, an arithmetic unit 10ab, and a register 10ac, and performs various arithmetic processing according to various programs loaded into the register 10ac. The input unit 10b is an input terminal, keyboard, mouse, touch panel, etc., to which data is input. The output unit 10c is an output terminal, display, etc., to which data is output. The communication unit 10h is a LAN card, etc., controlled by the CPU 10a which has loaded a predetermined program. Furthermore, RAM 10d is an SRAM (Static Random Access Memory), DRAM (Dynamic Random Access Memory), etc., and has a program area 10da where a predetermined program is stored and a data area 10db where various data is stored. Furthermore, auxiliary storage device 10f is, for example, a hard disk, MO (Magneto-Optical disc), semiconductor memory, etc., and has a program area 10fa where a predetermined program is stored and a data area 10fb where various data is stored. Furthermore, bus 10g connects CPU 10a, input unit 10b, output unit 10c, RAM 10d, ROM 10e, communication unit 10h, and auxiliary storage device 10f so that information can be exchanged. CPU 10a writes the program stored in the program area 10fa of auxiliary storage device 10f to the program area 10da of RAM 10d according to the loaded OS (Operating System) program. Similarly, CPU 10a writes various data stored in the data area 10fb of auxiliary storage device 10f to the data area 10db of RAM 10d. Then, the address on RAM10d where this program and data are written is stored in register 10ac of CPU10a.The control unit 10aa of the CPU 10a sequentially reads these addresses stored in register 10ac, reads programs and data from the area on RAM 10d indicated by the read addresses, sequentially has the arithmetic unit 10ab execute the calculations indicated by the programs, and stores the calculation results in register 10ac. This configuration realizes the functional configuration of the signal estimation device 11.

[0057] The program describing this process can be recorded on a computer-readable recording medium. Examples of computer-readable recording media are non-transitory recording media. Examples of such recording media include magnetic recording devices, optical discs, magneto-optical recording media, and semiconductor memory.

[0058] Furthermore, this program may be distributed, for example, by selling, transferring, or lending portable recording media such as DVDs or CD-ROMs on which the program is recorded. Alternatively, the program may be stored in the storage device of a server computer and distributed by transferring it from the server computer to other computers via a network. The program describing this processing (computer program) may be included in a computer program product.

[0059] A computer executing such a program may, for example, first store the program recorded on a portable storage medium or a program transferred from a server computer in its own storage device. Then, when processing is to be executed, the computer reads the program stored on its own storage medium and executes the processing according to the read program. Alternatively, the computer may directly read the program from the portable storage medium and execute the processing according to that program, or it may sequentially execute the processing according to the received program each time a program is transferred to it from a server computer. Furthermore, the processing may be executed by a so-called ASP (Application Service Provider) type service, where the processing function is realized only by execution instructions and result acquisition, without transferring the program from the server computer to this computer. Furthermore, the processing may be executed using a so-called SaaS (Software as a Service) type service, where a part of the server computer is made available to the user along with the program. In this form, the program includes information used for processing by an electronic computer that is equivalent to a program (data that is not a direct instruction to the computer but has the property of defining the computer's processing).

[0060] Furthermore, in this configuration, the device is configured by executing a predetermined program on a computer, but at least a part of these processes may be implemented in hardware.

[0061] [Other variations] It should be noted that the present invention is not limited to the embodiments described above. For example, in the first embodiment, an example was shown in which the process of step S113 is performed, followed by the process of step S114, followed by the process of step S115, and followed by the process of step S116. However, the order of at least some of the processes of step S113, step S114, and step S115 may be changed. In addition, the various processes described above may not only be executed in chronological order as described, but may also be executed in parallel or individually as needed, depending on the processing capacity of the device performing the processing. Furthermore, the Frobenius norm may be used instead of the L1 norm described above. Needless to say, other modifications can be made as appropriate without departing from the spirit of the present invention. [Explanation of symbols]

[0062] 11. Signal Estimation Device 113 Estimated Original Signal Optimization Unit 114 Graph Optimization Section 115 Graph Optimization Section 117 Time Domain Conversion Unit

Claims

1. A signal estimation device having an estimated original signal optimization unit that optimizes the estimated original signal according to an objective function based on an estimated original signal, a sample signal based on observation signals obtained by one or more sensors, missing information representing missing values ​​in the sample signal, a first graph representing the positions of the sensors and the relationships between the sensors, a second graph representing the discrete frequencies of the sample signal and the relationships between the discrete frequencies, and a third graph representing the relationship between discrete times and adjacent discrete times.

2. A signal estimation device according to claim 1, The estimated original signal optimization unit optimizes the estimated original signal while fixing the first graph and the second graph according to the objective function, A first graph optimization unit optimizes the first graph while fixing the estimated original signal and the second graph according to the objective function, A second graph optimization unit optimizes the second graph while fixing the estimated original signal and the first graph, according to the objective function, A signal estimation device having the following features.

3. A signal estimation device according to claim 1 or 2, A signal estimation device in which the objective function is a function that represents the fidelity of the estimated original signal and the missing information to the sample signal, the smoothness of the estimated original signal with respect to the product graph of the first graph, the second graph and the third graph, the low rank of the estimated original signal, the connectivity of the first graph, the connectivity of the second graph, the sparsity of the first graph and the sparsity of the second graph.

4. A program for causing a computer to function as a signal estimation device according to claim 1.

Images