Continuous-Time and Discrete-Time Singular Value Decomposition of an Impulse Response Function

This paper proposes the continuous-time singular value decomposition (SVD) for the impulse response function, a special kind of Green’s functions, in order to find a set of singular functions and singular values so that the convolutions of such function with the set of singular functions on a specified domain are the solutions to the inhomogeneous differential equations for those singular functions. A numerical example was illustrated to verify the proposed method. Besides the continuous-time SVD, a discrete-time SVD is also presented for the impulse response function, which is modeled using a Toeplitz matrix in the discrete system. The proposed method has broad ap-plications in signal processing, dynamic system analysis, acoustic analysis, thermal analysis, as well as macroeconomic modeling.


Introduction
Impulse response function (IRF) has been considerably used for signal processing, dynamic system analysis, acoustic analysis, thermal analysis, and macroeconomic modeling. The significance of such function is that for any input, the output can be calculated in terms of the input and the impulse response. To determine an output directly in the time domain requires the convolution of the input with the impulse response, which is very complicated. Therefore, it is usually to calculate the Laplace transform of a system's output by the multiplication of the Laplace transform of the IRF with the input's Laplace transform in the frequency domain at first; and then find the output in the time domain through an inverse Laplace transform.
This study proposes a direct method to calculate the system output through the singular value decomposition (SVD) algorithm. As a popular technique of factorizing a matrix, a variety of forms of SVD have been developed recently and applied for solving a wide range of scientific and engineering problems. For exam, Ramirez-Velarde et al. used SVD on the video correlation matrix to perform the principal component analysis [1]. Liu et al. developed an SVD-based minutia matching method for finger vein recognition, which performs better than other methods [2]. Wan et al. proposed a quasi SVD method to extract algebraic features from human face with enhanced face recognition rate and speed [3]. Zhang et al. applied a higher-order SVD method for denoising magnetic resonance data to improve the inspection quality and reliability of quantitative image analysis [4]. Shih et al. proposed an adaptive parametrized block-based SVD for preserving the edge structure and avoiding blurred image after the compressing process [5]. Ghazy et al. presented a block based watermarking scheme using the SVD algorithm to embed encrypted watermarks into digital images. Experimental results showed that the SVD-based method is superior to the traditional method for embedding encrypted watermarks and extracting them efficiently under attacks [6]. Cai et al. used SVD for model predictive control design of distribution systems and demonstrated the advantages of their method through two case studies [7]. Katsaounis et al. derived a new way to detect combinatorial equivalence of symmetric factorial designs based on the SVD of design matrices to provide a fast screen for non-equivalence [8]. Berenguer et al. used SVD of successive iterated solutions on subdomains interfaces to computing a low-rank approximation of the Aitken's formula at a high computational efficiency. The results were then applied for solving heterogeneous 3D groundwater flow problems [9]. di Dio applied SVD for analyzing temperature dependent radial distribution functions and showed that the generality of the SVD approach allows investigating various temperature dependent structures in a unified way [10]. Abuturab proposed a new color image security system based on SVD in gyrator transform domains and proved robustness and other advantages of the developed method in enhancing the security [11]. Kavaklioglu applied SVD for modeling Turkey's electricity consumption, in which the SVD was used to reduce the dimensionality of the problem and to provide robustness to the estimations [12]. However, in the previous studies, SVD was only applied for discrete systems (matrices) and this paper aims to present a continuous-time SVD approach for the first time to facilitate the calculation of convolution of a set of singular functions and the IRF  The analysis below shows how to find a set of "singular" functions ( ) The left hand side of Equation (1) is Duhamel's integral [13] for an input The complete solution process for Equation (2) includes three steps: 1) simplification of left hand side of Equation (2); 2) simplification of right hand side of Equation (2); 3) solve the simplified Equation (2).
Step 1: Simplifying The integral on left hand side of Equation (2) can be solved as: The first term on the right hand side of Equation (3) can be simplified as: sin atan e sin e cos e 1 1 In order to find j φ to satisfy Equation (3), the second term on the right hand side of that equation can be set to zero: Solutions of the above equation are: Substituting Equations (4) and (5) into (3), a more compact form can be obtained as: Replace j ω with j φ using Equation (6) and the left hand side of Equation (2) can be eventually simplified as: Step 2: Simplifying (6) on the right hand side of Equation (2) Step 3: Solve Equation (2) In order to satisfy Equation (2), j φ , j ω and j D must meet some requirements. This step completely solves Equation (2) and discusses the characteristics of j φ , j ω and j D from the final solution. Apply the obtained simplified terms to the original problem by substituting Equations (8) and (9) into Equation (2), we can have: The term ( ) e 1 n t − can be canceled from both sides, which implies that the integer n will not affect the final results. Equation (10) becomes (11) can be expressed expanded as: Because the left hand side of Equation (11) Equation (13) is a transcendental equation with infinite roots so it is reasonable to label the roots in ascending order such that: 0 1 2 ω ω ω < < < etc. Thus, Equation (11) can be further reduced to: which yields For an integer j.
As can be seen from Equations (13) and (15), the value n does not affect the results for either j D or j ω , for convenience we set n j = to further simplify the solution (Equation (16)) as: Thus, Equation (2) has been completely solved with solutions provided through Equations (13), (15), and (17). Nevertheless, 0 j ω = is excluded even it is a solution of Equation (13) because it will lead to a trivial solution for Equations (2) and (17). Thus, the solved singular functions can be represented as and conjecture 1 is proved.

Illustrative Example
An illustrative example is presented here to validate the developed method.
In Equation (13), we assume 10 T = and the roots of that equation can be found using Mathcad, which are the intersections between function curves ( ) From Figure 1, the roots for Equation (13)  Equations (15) and (17) provide values for j D and j ϕ . Let 0,1, 2,3 j = , substitute j and j ω solved from Figure 1 into Equation (15) and (17)    . Note that T is not the period of the sine waves. It needs to be noted that the numerical example in this paper was modeled and solved using Mathcad, the math software for engineering calculations.

Discrete-Time SVD of Transformable to Symmetric Matrices
As a counter part of the continuous-time SVD, a discrete-time SVD is also conducted for the Green's function ( ) e t τ − − . In discrete system, such function is represented as a Toeplitz matrix. The entire SVD process is presented as follows.
Given a matrix A with the following property: where R and L are suitable orthogonal matrices. It is obvious that the product

SVD of a Lower Triangular Toeplitz Matrix
Let P be a permutation matrix used to reverse the order of the rows. Then P must be an orthogonal matrix with ones along the main antidiagonal line and zeros elsewhere. Thus if the matrix A is a Toeplitz matrix, then the product P A * must be a Hankel matrix. In Equation (22), if we assume R P = , A is a Toeplitz matrix, and L is the identity matrix L I = , then according to SVD theory, the matrix S should be a matrix obtained by assembling orthonormal eigenvectors from the Hankel matrix P A * and Λ is a diagonal matrix with the corresponding eigenvalues of P A * along the diagonal line. This shows that the singular values and the right singular vectors of a lower triangular Toeplitz matrix are just the eigenvalues and eigenvectors of the Hankel matrix obtained by flipping upside-down the lower triangular Toeplitz matrix. However, it needs to mention that even the eigenvalues of P A * can be positive or negative, the singular values of A must be positive.
The following process explains how to find U and V in Equation (24) for the lower triangular Toeplitz matrix A.
The eigenvalue and eigenvector for the Hankel matrix P A * is: Since i λ can be either positive or negative but the singular value i s must positive, we have:

( )
, and For example, if matrix A is:

Resolve the Example in Section 3 Using Discrete-Time SVD
In this section, the illustrative example solved in Section 3 will be resolved by the discrete-time SVD algorithm. For a discrete decomposition, the original continuous function is first discretized as 100 discrete values, from which a Toeplitz matrix is formed and the validated SVD process presented in Equations (22) to (26) is applied to find finding U, V and D for such Toeplitz matrix. Since the SVD algorithm has been implemented into Mathcad, so the entire process is solved using that software and the detailed codes are presented below, where the embedded figure shows the discretized function e t y − = , which is displayed  Figure 4 compares the first column of U with the discretized normalized function 1 f (Equation (18)). From that figure and Figure 2, it is confirmed that the discrete-time SVD can correctly find the output for the IRF and the first column of the SVD matrix U is the discrete solution for such problem, which can be considered as a counterpart of Equation (19) in discrete domain.
In Figure 3 and Figure 4, the is the ith column extracted from the matrix U, which is obtained following the discrete-time SVD process presented in Figure 3 and from Equations (22) to (26).

Discussion
In this section, the discrete-time SVD is performed to solve the problem raised in Equation (2) in a discrete system. From the solution it can be seen that the presented method can accurately yield the results, which is the ith column of the SVD matrix U. The critical technique employed here is to flip the lower triangle Toeplitz matrix upside down to convert it to a Hankel matrix through a permutation matrix. Comparing to the Toeplitz matrix, which is usually used to model the IRF, the Hankel matrix is easier to be decoupled therefore making computing effort easier and more efficient. The same "flipping upside-down" technique is also applied for the continuous-time SVD (the function term (1)) to simplify the solution process.    (2) and then it can be solved following the same approach discussed in Section 2 (Equations (3) to (17)).

Conclusion
This study presented and validated continuous-time and discrete-time SVD for the impulse response function, a special kind of Green's functions. For both the continuous-time and discrete-time SVD, the "flipping upside-down" strategy is employed to simplify the solution process. Illustrative example confirms the accuracy and efficiency of the present methods. Even though this paper mainly focuses on the SVD of the Green's functions, it is obvious that the solution can follow a similar approach. Comparing to existing methods, one advantage of the presented method is that by using such method, the time-domain output can be directly determined from the convolution of the input with the IRF without going through Laplace and inverse Laplace transforming. Considering the diversity of impulse response functions involved in multidisciplinary engineering problems, the method presented in this study can be extensively developed for broad applications. In the future, the SVD approach can be further modified to absorb some features of Adomian decomposition method [14] [15] [16] for being used to solve more engineering problems such as Cauchy type singular integral equation [17] [18], Bagley-Torvik equation [19], Fredholm integro-differential equation [20] [21] [22].