Signal Processing Interpolation and Problem Conditioning Educational Workbench

This work presents a newer and more complete version of an educational tool to be used in signal processing interpolation-related subjects. Besides the consolidation of acquired theoretical knowledge, the tool allows now its users to apply three error geometry patterns, test minimum or maximum dimension signal reconstruction algorithms and problem conditioning through the analysis of the matrix spectral radius or the condition number: these new features gives the possibility to alter the problem definitions to the desired goal before the reconstruction begins. The time unit that measures the algorithms performance is (nlogn) thus independent from the machine’s architecture. A video of the developed tool can be seen in: https://www.dropbox.com/s/t6yiiuy31ramxse/FILME_1.avi.


Introduction
Nowadays, namely in engineering areas, there are a great amount of subjects that require a large and deeper knowledge of theoretical aspects related to a specific scientific area.Sometimes students don't really understand the need to have skills on important concepts and lack the perception to understand its utility or even if they really work.During the last decades, with the quick spread of information technologies, some new tools arose whose purpose is to ensure and consolidate the theoretical concepts taught in regular teaching classes.The main idea is to give students a way to experiment concepts taught in classes, as well as giving them the opportunity to test them in different scenarios.With these tools students and other kind of users can practice and enhance skills by themselves as self-taught persons.
In this paper a new version of the "Signal Processing Interpolation Educational Tool" (SPIEW) is presented (Costa et al., 2012).In the first version we propose a new and easy-to-use tool which will help students consolidate concepts of signal processing and complement their knowledge by enabling graphical results after using complex signal reconstruction algorithms.The tool permits to achieve calculation in just one hit or step by step through the iterative reconstruction refinement process and the assessment of the reconstructed signals.
Now we present an improved version of SPIEW.Its main new features are: the possibility to apply three error geometry patterns (interleaved, random and burst); control the redundancy of the signals by applying a proper oversampling factor, r, in the band limiting operation; test two different classes of signal reconstruction algorithms and the usage of the iterative and direct calculation methods; the previously set up the problem conditioning through the analysis of the spectral radius, (S), of the iteration matrix (in the case of iterative methods are used) and the computation of the condition number, k(A), in the case of the direct methods are used.In this way, its user can employ these parameters to know in advance if the problem will converge to an accurate solution, even before running one of the reconstruction methods available.A previous analysis of key parameters permit to chose the most suitable method to be used in the signal reconstruction test scenarios.
Since all the chosen methods use the FFT computation, (nlogn), the time necessary to compute a FFT is used as reference.So, the performance is measured through a relative Time Unit (TU) that is the time to compute a FFT in a particular machine.

Reconstruction Algorithms
The signal reconstruction algorithm used in SPIEW was the discrete version of the Papoulis algorithm (Papoulis, 1975): the Papoulis Gerchberg algorithm (Ferreira, 1994a;Jones, 1986).A thorough description of this maximum dimension algorithm can be found in (Costa et al., 2012).The Papoulis algorithm is still available to run but a new minimum dimension reconstruction algorithm was included like the one used in (Neves et al., 2008).Two different methods can be used to implement this algorithm: iterative and non-iterative calculation.We also call the non-iterative method the direct calculation method.This way students can test and compare maximum and minimum dimension algorithms with direct and iterative formulations.In the next subsections the minimum dimension algorithm is presented.

Minimum Dimention Algorithm
This subsection describes a minimum dimension algorithm which also requires band-limited signals of finite-dimension similarly to the Papoulis Gerchberg algorithm.
The minimum dimension algorithm is a system of only l equations corresponding to the l unknown samples when the total number of samples is n and n> l.
To establish the basic concepts, the specific case of an original signal x [n] For this signal, Equation (1) becomes

… (2)
For reconstruction purposes let us assume that the 2 nd and 4 th samples of x[n] are lost.Then, because only the unknown samples are to recover, the set of equations, Equation (2), will be limited to those including the unknown samples.In each of these equations, we are interested in separating the right side terms of Equation (2) containing the unknown samples (x 2 , x 4 ) from those containing the known ones.These yields, , which is equivalent to Let us denote by u the subset of the original signal x[n] which contains the unknown values.In this case, u={x 2 ,x 4 } is of cardinality l=2.Also, let us define U = {i 1 , ..., i l } as the set of subscripts of l unknown samples in x [n].In the present case, U = {2, 4}.Therefore, Equation (4) can be written as or, in matrix form where S is a ll principal submatrix of B, as defined in Equation ( 4), and h is the (n-l) dimensional vector in the second sum of Equation ( 5), which is a linear combination of the known samples of x [n].S is the system matrix (Ferreira, 1994b).Theoretically, Equation ( 6) has a unique solution regardless the number and distribution of the lost samples.
On the one hand, by successive operations, Equation ( 6), becomes equivalent to , and a basis to implement the non-iterative method is so established.On the other hand Equation ( 6) suggests the needed formula to the iteration process for the case where a non-relaxation method is used.Then u (i) is obtained at iteration i and the solution is given by the limit regardless of u (0) .
Direct calculation of u as given in Equation ( 7) has the advantage of being done in one single step and gives an exact solution, providing that (I-S) -1 exists.The drawback is that in case (I-S) -1 does not exist or the calculation is not possible and no reliable approximation is possible.In practice, there are several factors which may lead to serious difficulties in calculating the inverse of I-S.For example, if one of the eigenvalues of S,  i (S), is close enough to the unity, then the computation of (I-S) -1 may become very difficult, even impossible, which leads to an ill conditioned problem.In such cases an iterative method is advantageous because it may be used to circumvent this difficulty and to find a relatively accurate approximation for the solution, u.Despite the fact of having an ill conditioned problem, in direct calculation method such problem is impossible to solve whereas in iterative methods an approximation is always possible to be found, though its accuracy may not be very high.
The conditions under which these equations provide a solution for u can be found in (Ferreira, 1994a) and will be discussed in the next subsection.

Convergence Analysis
It is necessary that (I-S) -1 exists or the iterative processes converge to find a problem solution.It is demonstrated in (Ferreira, 1994a) that for low pass signals of dimension n with s known samples and q nonzero harmonics, reconstruction iterative algorithm converges if: which is equal to It is important to refer that q=2m + 1 where m is the sample bandwidth of the Low Pass Filter (LPF) (Costa et al., 2012).So Equation ( 11) can be rewritten as follows: which establishes a sufficient condition on the density of the known samples.More precisely, the number of known samples s must exceed twice the number of nonzero harmonics.For low pass signals, the condition is consistent with the well-known classical sampling theorem.
The sufficient condition expressed by ( 12) is obviously necessary too.After all, low-pass signals with q=2m+1 nonzero harmonics are described by 2m+1 parameters, and therefore one cannot expect to reconstruct them from less than 2m+1 samples.
In order to use the algorithms available in the new tool it is necessary that all signals are band limited.This means that Equation (12) establishes a necessary condition in our tool.Because only low pass signals are used, the necessary condition is met.
For the specific case of the Papoulis-Gerchberg algorithm, a relaxation constant, , is introduced to improve the convergence rate, as stated in the previous tool version.In (Ferreira, 1994a) it is demonstrated that the algorithm converges if both conditions are satisfied: 0<<2 and 2 1 Equation ( 12).Because we have considered =1 and used low pass signals the iterative process converges.
Regarding minimum dimension method, the characterization of the system matrix, S, is fundamental to infer about the existence of solution to the problem formulated in Equation ( 2).Namely its eigenvalues,  i (S), and in particular its spectral radius, (S), condition the problem.(S) values close to 1 (the unity) make the computation of (I-S) -1 very difficult, even impossible.Theoretically, the condition (S) < 1 would always guarantee the solution possible.However, due to the some computational representation constraints, small variations on input values can cause big variations on the output and cause the problem to be ill conditioned.In these cases a new implemented parameter must to be taken into account: the condition number of (I-S) -1 , presented latter below.
According to these conditions, it is possible to put a reconstruction problem into a point such that it is well conditioned.It is known that the eigenvalues of the system's matrix S depend on the distribution of the missing samples.In particular, its spectral radius is more likely to be unitary for burst distributions rather than for equidistant missing samples (Ferreira, 1994c).In the case of signal reconstruction it is interesting to note that, if the distribution of the missing samples U = {i 0 , i 1 , …, i l-1 } is equidistant by some fixed integer p1, that is, U = {i 0 p, i 1 p, …, i l-1 p}, then the eigenvalues  i of S are given by rp /p, i.e., where rp denotes the greatest integer less than or equal to rp and rp denotes the smallest integer equal or greater than rp.In the particular case of r=rp /p, the eigenvalues of S are all the same  i (S)=r, i.
Given the above analysis, it is possible to put the problem into a well-conditioning point by properly selecting the gap between the missing samples.Then it is possible to put  i (S) close to either r or its multiples, regardless of the number of missing samples.By using an appropriate choice of the oversampling and interleaving factors, r and p, respectively, it is possible to put  i (S) less enough than the unity in order to control the reconstruction accuracy and processing speed (i.e., such that mr is an integer).
Figure 1 shows the maximum and the minimum eigenvalues of S as a function of the interleaving factor, p, for a given bandwidth, defined by r=0.6.As the figure shows, greater values of p lead to better-conditioning problems because  max decreases as p increases.Also, when the product mr is an integer, all eigenvalues are equal since they are  i (S) = r, as stated before.In Figure 1, this occurs for p=5 and p=10.To better understand and control the problem conditioning, in the presented tool, some new parameters are calculated: the spectral radius of the S matrix, delimited in Equation ( 13) and the Condition Number of the (I-S) matrix (Equation ( 14)) (Golub & Van Loan, 1996;Kincaid & Cheney, 2002).Through the analysis of both values, and also the convergence condition, Equation ( 12), the user can know in advance if the reconstruction algorithm will converge to an accurate solution or not.Therefore through the properly chosen condition point the signal reconstruction is guaranteed.
In practice, when (S) tends to be very near or even equal to the unity, the problem may become ill-conditioned, even if it is theoretical well-conditioned due to the fact that (S)<1.Condition numbers near or above 10 15 do not generally lead to well-posed problems, even with double precision arithmetic (64-bit floating point, in many machines).The condition number of matrix A, k(A), is calculated according the equation presented in (Kincaid & Cheney, 2002)   Iteration Residual Error can also act as comparative values between the iterative methods performance.In the direct method, instead of the iteration error evolution graphic, both the reconstruct and corrupted signal versions are shown.
This application is to be used in Digital Signal and Image Processing related subjects.So far, the users' receptivity has been encouraging.Currently we have used in this tool high frame rate sequences and still images obtained with the MVD method (Rusanovskyy et al., 2013), still we believe it can be used in images from different areas, such as medical, e.g.
As future work we aim to add new features to SPIEW (Karthik & Prabhu, 2010) and conduct more tests with graduated students (MSc & PhD).The aim is finding the critical points that establish the appropriate method according to a previous analysis of the problem conditioning.