Signal detection by means of orthogonal decomposition

Matched filtering is a well-known method frequently used in digital signal processing to detect the presence of a pattern in a signal. In this paper, we suggest a time variant matched filter, which, unlike a regular matched filter, maintains a given alignment between the input signal and the template carrying the pattern, and can be realized recursively. We introduce a method to synchronize the two signals for presence detection, usable in case direct synchronization between the signal generator and the receiver is not possible or not practical. We then propose a way of realizing and extending the same filter by modifying a recursive spectral observer, which gives rise to orthogonal filter channels and also leads to another way to synchronize the two signals.


Introduction
This paper presents methods whose development was inspired by attempts [1,2] to automatically and continuously monitor the integrity of the signal acquisition chain of the new beam loss monitoring (BLM) system [3,4] currently being deployed at CERN.1 A similar solution, not capable of continuous operation, is deployed at CERN's flagship particle accelerator, the Large Hadron Collider (LHC) [5]. The BLM systems are complex measurement systems with many channels, containing a multitude of detectors, cables and connectors. In such measurement systems, one may wish to automatically detect problems such as manufacturing defects, installation issues or component degradation: cut or disconnected cables, bad contacts at connectors or terminating impedances, aging or radiation-related effects in components.
The integrity of such a signal chain may be assessed by applying a known excitation to its input(s) and measuring the resulting response at its output(s). In particular, this excitation could be a periodic test signal designed and generated by the user -this is the case we consider in this paper. One can then detect whether the expected response is present at the output of the signal chain. Nowadays, because of the ubiquity of microprocessors and reprogrammable logic devices, this is usually done in the digital domain.
Matched filtering [6] is frequently used for detecting the presence of patterns due to its simplicity. It involves calculating the cross-correlation of the input signal with a template carrying the known pattern and identifying the maxima of the resulting function. The matched filter has been shown to be the optimal linear time invariant filter for maximizing the signal-to-noise ratio in signals contaminated by additive stationary white noise. Application examples include image processing in transmission radiography [7] and in the measurement of swirling fluidized beds [8], spectroscopic measurement of electron source stability [9] and plasma rotation velocity measurement at a tokamak [10].
A pattern match is indicated by the characteristic cross-correlation peaks, as illustrated in figure 1. This means that a definitive decision concerning the presence of the pattern is only possible at the occurrence of the peak, when the template is synchronous with the pattern occurring in the input signal.
Since we assume a periodic excitation, this result can be improved by circularly rotating the samples of the template. A rotation of one template sample at the arrival of each input sample maintains a given alignment between the input signal and the template. If this alignment synchronizes the two signals, the value of the cross-correlation peak can be obtained in each time step. The resulting filter can be interpreted as a time variant (LTV) matched filter, whose coefficients are modified in each time step.
Let us consider the periodic pattern as the linear combination of cosine components with fixed frequencies and prescribed amplitudes, periodically obeying a given phase relation at certain instants during their time evolution. If the input signal deviates significantly from any of these parameters, no pattern match should be reported. The orthogonal filter channels presented in this paper provide a way of identifying such situations.
After briefly recalling the principle of matched filtering, section 2 describes a recursive method for the implementation of a time variant matched filter. A solution is then suggested to synchronize the template to the input signal of the receiver in order to obtain the cross-correlation peak in every time step, if using a synchronization signal is not possible or the related costs are too high. The rest of the document shows how the same result can be obtained and extended by incorporating the test pattern into a recursive observer [11]. The observer itself is introduced briefly in section 3. The suggested modification, the resulting orthogonal filter channels and another possibility for synchronization are presented in section 4. The operation of the filters is demonstrated in section 5, then section 6 concludes the paper.

Time variant matched filtering
A matched filter [6] is a finite impulse response (FIR) filter suitable for detecting the presence of a pattern in a noisy signal. Its impulse response consists of the samples of the pattern in reverse  The first row of figures shows the cross-correlation between 6 periods of a sine function of 3 Hz as pattern, with the same signal repeated twice and contaminated with white noise at SNR = 0 dB. The second row of figures is obtained in a similar manner but using a linear swept-frequency sine (chirp) from 0 Hz to 5 Hz over 2 s. The third row is a combination of the previous two: it shows the cross-correlation of the ideal cosine as pattern with the noisy chirp as input signal. The characteristic cross-correlation peaks are clearly visible in the first two output signals. Note that the crosscorrelation function of the two non-corresponding signals in row 3 has a considerably lower amplitude than the other two, indicating the absence of the pattern in the input signal.
order, thus the output of the filter is the cross-correlation of the input signal and the template. In the following, this principle is recalled briefly, then the suggested LTV filter is introduced along with its properties.
The indexes in some of the following equations may appear excessively complicated. This ensures their consistency with the subsequent formulas.

Matched filtering
If the input signal and the template are denoted by y n and h n , respectively, the samples of the corresponding cross-correlation function are expressed as where the limits for l are given by the finite length N of the template h n : h n ≡ 0 if n [0, N). Figure 1 shows a few examples of cross-correlation waveforms. The peaks thus obtained can be used for detecting the presence of the pattern and, if present, identifying its location in time.
A circular buffer of depth N is required for storing the preceding samples of y n and N multiplyand-accumulate operations need to be performed at the reception of each new y n sample to obtain the next value of v n .

Using a time variant pattern
If the template contains a finite number of periods of a periodic signal, a constant alignment between the input of the matched filter and the template can be maintained by using a time varying pattern. This requires each y p value to be multiplied by the same h q coefficient in all sums it is present in. The following filter responds to this requirement: With a periodic input signal, this filter produces the same point of the cross-correlation function in every time step, as illustrated in figure 2. The values yielded by eq. (2.1) and eq. (2.2) correspond if n = m × N, m ∈ Z. Thus, eq. (2.2) can be interpreted as a linear time variant matched filter.
Without loss of generality, we can assume y n = 0 for n < 0, which also implies x n = 0 for n < 0. This means that eq. (2.2) can be expressed recursively as 3) see section 2.3. The block diagram of the corresponding recursive LTV matched filter is shown in figure 3. The evaluation of eq. (2.3) also requires a circular sample buffer of depth N so that y n−N is available for the calculation of x n . Nevertheless, the calculation of each new x n value only requires one multiplication and two additions, thus the corresponding computational burden is reduced considerably with respect to eq. (2.1) and eq. (2.2), which are identical in this respect.

Properties of the recursive LTV matched filter
Equation (2.3) describes a moving window filter with a window length of N samples. We assumed that y n = 0 for n < 0. This means that y n−N ≡ 0 while n ≤ N − 1, thus during this time, x n accumulates the samples of y n multiplied by the corresponding samples of h n . Afterwards, a new y n gets added to the buffer in each time step, while the contribution from y n−N at the other end of the -4 -  window is cancelled since h n mod N = h (n−N ) mod N . Accordingly, the recursive and non-recursive formulations of the LTV matched filter, along with x n in eq. (2.2) and x n in eq. (2.3), are equivalent. As a result, the values yielded by eq. (2.1) and eq. (2.3) also correspond if n = m × N, m ∈ Z. For other values of n, however, the LTV matched filter behaves differently. This is due to the fact that the matched filter yields values of v n for all possible alignments between h n and y n successively, which is responsible for the typical peaks characterizing the cross-correlation function in the case of a match, as shown in figure 1.
In contrast, the recursive LTV matched filter described in section 2.2 keeps the alignment between y n and h n constant. Therefore after convergence, a single point of the cross-correlation function, corresponding to this fixed alignment, is obtained for the N most recent samples of y n whenever a new sample is received. As illustrated in figure 2, this means that we get the value from the same single position along the cross-correlation function in every time step. This may be an advantage or a disadvantage depending on the corresponding alignment between y n and h n : if the two are synchronized, the filter is going to yield the maximum cross-correlation peak value at every new input sample, but if they are not, it will consistently return a value under the maximum.

Synchronizing the filter
In order to find the best alignment between the input signal and the template, the recursive filter in eq. (2.3) can be extended to a filter bank. The filter bank containing N instances of the filter is described as follows: with 0 ≤ k < N. The filter bank realizes all possible alignments between h n and y n .
After convergence at n = N − 1, the highest output value can be used to identify channel l, the channel best synchronized to the input signal. Naturally, this is only reasonable if this channel indicates the presence of the pattern. At this point, all other channels of the filter bank can be disabled and channel l can continue operation as synchronized recursive LTV matched filter.

On the computational complexity of synchronization
Depending on the computing resources available and the requirements of the application, the implementation of the filter bank can be scaled broadly. All filters in the bank could have their own dedicated computing resources and operate in parallel. Alternatively, if the input sample rate and processing frequency permit it, the same computation resources could be shared among several channels and the computations might be carried out sequentially in a time-multiplexed manner. We might also consider using only one set of computing resources and operating one filter from the bank at a time, thereby obtaining one filter value per N input samples, increasing the time required for synchronization considerably.
After a successful synchronization, if the alignment between h n and y n can be assumed to be constant, the superfluous channels in the filter bank can be disabled permanently. If this is not the case, a periodic total recalibration may be desirable. In case the alignment is expected to drift slowly, some of the filter bank channels adjacent to the operational channel -such as l − 2, l − 1, l, l + 1 and l + 2 -may be kept active. If the highest filter output value shifts systematically to one of these additional channels, the channel in question can take over the role of the operational channel, thereby allowing the processing to maintain synchronicity.

The conceptual signal model
The recursive observer [11] is based on a conceptual signal model, shown in figure 4. The signal under measurement is assumed to be generated by this hypothetical dynamic system, and the observer estimates the parameter vector of the measured signal by estimating the state vector of the conceptual signal model.
Equations The LTV conceptual signal model is described by the following system equations: where x n is the state vector of the conceptual signal model, while y n is its output, that is, the signal under measurement. f 0 is the relative frequency of the fundamental harmonic with respect to the sampling frequency: where f 1 is the fundamental harmonic frequency and f s is the sampling frequency. The number of harmonic components K is such that K × f 1 < f s 2 ≤ (K + 1) × f 1 holds. In case of equality, i = −K, . . . , K + 1 and N = 2K + 2 in the equations above. This corresponds to a single additional channel at f s 2 .
-7 -2018 JINST 13 P03002 The coupling coefficient vectors c n determine the basis functions generated by the channels of the conceptual signal model.
In the case of the Fourier analyzer, the conceptual signal model is effectively a complex multisine generator reconstructing the input signal of the observer from its DFT. The state vector x n consists of the N = 2K + 2 complex Fourier coefficients, including DC and f s 2 . If y n is real-valued, the Fourier coefficients form complex conjugate pairs: x i,n = x * −i,n , where · * is the complex conjugate operator. The coupling coefficients are set to unit complex numbers according to eq. (3.5), thus each complex conjugate channel pair generates the samples of a single sinusoid.

The observer
The state variables of the conceptual signal model introduced in section 3.1 can be estimated by an appropriately designed observer [13]. Thus, in the case of the Fourier analyzer, this structure can estimate the complex Fourier coefficients of its input signal in its state variables. Figure 5 shows the block diagram of the observer designed to match the LTV conceptual signal model.
The system equations governing the operation of the LTV observer are: x n+1 =x n + g n e n =x n + g n (y n −ŷ n ) , (3.6) y n = c T nx n , (3.8) wherex n is the estimate of the state vector,ŷ n is the estimated signal value and e n is the estimation error.

JINST 13 P03002
The c n and g n coupling vectors determine the behavior of the observer. In the case of the Fourier analyzer, the coupling vector g n is set according to eq. (3.10): it is the product of the observer gain α N , a tunable parameter setting the dynamical behavior of the observer, and the coupling vector c n .
If the coupling vectors are adjusted as in eq. (3.5) and eq. (3.10), with α = 1 and f 1 = f s N , the result is a dead-beat observer. This means that the transients of the observer decay in at most N steps and afterwards, its state vector contains the recursive DFT of the N most recent samples of the input signal.

Dead-beat observers
Dead-beat settling means that the state variables of the observer converge to those of the conceptual signal model, and the estimation error e n reaches 0 in N (or fewer) time steps. As discussed in section 3.2, with appropriate settings, the Fourier analyzer provides dead-beat convergence. In a more general case, the sets of coupling vectors c T i and g i need to form a biorthogonal system for the recursive observer to exhibit dead-beat settling behavior [11].
This requires C T = G −1 with The vectors c T i and g i contain the values of the coupling coefficients corresponding to all channels in time step i. Conversely, the column vectors of C T and the row vectors of G consist of the coupling coefficient values corresponding to a given channel over an entire time period.

Incorporating the test pattern into the observer
If the recursive observer is operated in a configuration with dead-beat settling (see section 3.3), the overall transfer function of the structure in figure 5 from y n toŷ n is H P = z −N [14]. This means that any individual channel of the structure can be realized on its own without implementing the parallel channels, by using a structure akin to the recursive LTV matched filter shown in figure 3. Similar structures implementing a sliding DFT have been suggested in the literature [15][16][17].
In the following, we suggest a way of modifying the C T and G matrices of the Fourier analyzer so that one channel performs recursive matched filtering on its input signal with an arbitrary periodic template. This modification also yields orthogonal matched filter channels producing a measure of the "error" of the matched filter. Moreover, the dead-beat properties of the structure are also preserved, thus implementing an arbitrarily low number of channels instead of the entire structure is sufficient.

The modified Fourier analyzer
Let us start with a Fourier analyzer with a number of channels N matching the length of the desired template. Ideally, the harmonic components of the template should be harmonics of the resulting fundamental frequency f 1 = f s N of the observer (see section 3).
-9 - In this configuration, let T be the set of positive harmonic indexes contributing to the template: T = {t 1 , t 2 , . . . t m }. Let us denote their respective amplitudes as w i = A i e ϕ i , i ∈ T. All modifications described in the following need to be extended to the complex conjugate counterparts of the channels concerned.
We can modify the c coupling coefficients corresponding to channel t 1 so that this channel, combined with its complex conjugate pair, generates the template at the output. In order to do this, let so that 2 × Re c t 1 ,n = h n . As a result, the columns of the C T matrix are no longer orthogonal. In order to restore their orthogonality, one may resort to one of the methods used for QR decomposition such as the Gram-Schmidt process [18] or Householder reflections [19]. The method of choice needs to be applied to the columns of C T corresponding to the harmonics of the template. We used MATLAB's built-in qr function, which appears to rely on the modified Gram-Schmidt process, although this is not explicitly documented. The leftmost columns of the Q matrix thus obtained contain a set of orthonormal vectors that span the corresponding subspace of C N , and the first one is parallel to the initial modified basis vector calculated in eq. (4.1). It should be noted that the two vectors may be antiparallel, which should be corrected by a multiplication by −1 as required for best results later on. We scaled these vectors by i ∈T |w i | 2 so that they match the norm of the basis vector obtained in eq. (4.1), then copied them back into C T .
The resulting modified C T matrix is nonsingular, thus the corresponding g coupling coefficients yielding a dead-beat observer can be calculated as G = C T −1 .
It can be shown that the g coupling coefficients corresponding to channel t 1 in the resulting G matrix can be expressed using those of the initial G matrix as Since in the case of the Fourier analyzer, g i,n = 1 N c * i,n (see section 3.2), it follows from eq. (4.1) and eq. (4.2) that Re g t 1 ,n ∝ h n . This means that in the Fourier analyzer modified according to the above, Re x t 1 ,n ∝ x n in eq. (2.3), thus the same recursive LTV matched filter output can be obtained by this method.
If the g t 1 ,n coefficients and the input signal y n are not synchronized, the resulting crosscorrelation value does not correspond to the maximum of the function, as described in section 2.3. Naturally, the synchronization method based on a filter bank presented in section 2.4 is still applicable. Another method is proposed in section 4.2.
In this scheme, due to the orthogonality of the row vectors of G, the other channels in T yield x values whose real parts are the cross-correlation of the input signal with a signal consisting of the same harmonics as the template, but orthogonal to it. Along with Im x t 1 ,n , these may be used to supplement the use of Re x t 1 ,n as matched filter output value to better judge the "goodness" of the cross-correlation fit: higher amplitudes in these channels indicate more prominent components in the frequency bands of interest not matching the structure of the template.
-10 -Since the dead-beat behavior of the observer is preserved, any channel of the structure may be realized on its own without the need to realize the entire structure. One might be interested in realizing only the channels in T, for instance. Unless one wants to restore the signal component(s) generated by the channel(s) in question, the output coupling using the c t i ,n coefficients may also be omitted.

Synchronization in the modified Fourier analyzer
We propose a way of synchronizing the coefficient sequence g t 1 ,n with the input signal y n by estimating the time lag ∆t n between the two. The suggested method relies on the state variables of certain regular Fourier analyzer channels: those corresponding to the harmonics of the template. As mentioned earlier in section 4, these FA channels may also be realized on their own without implementing the entire feedback loop.
The phase measured in the ith regular Fourier analyzer channel (i ∈ T) in the nth time step can be expressed as: where ϕ i = arg w i is the phase of the corresponding harmonic component of the template, f i is the (absolute) frequency of the harmonic and 2πn i,n accounts for the fact that we only get phase information mod 2π from the Fourier analyzer. The estimate for ∆t n that can be calculated from the phase of channel i is then: However, we do not know the value of n i,n . Nevertheless, we are expecting to get the same value for ∆t n in all channels in T, thus ∆t n = ∆t t 1 ,n = ∆t t 2 ,n = . . . = ∆t t m ,n . Going from ∆t t 1 ,n = ∆t t 2 ,n to ∆t t m−1 ,n = ∆t t m ,n , we obtain a system of m − 1 Diophantine equations with m integer variables in n i,n : 1 2π where f T is the reciprocal of the period of the template: f T = 1 T T , thus each f t i f T is an integer. Since the right-hand side of each Diophantine equation is an integer, it is desirable to round the left-hand side to the nearest multiple of gcd to reduce the contribution from noise in the measured phases. We can then solve the system of equations using, for instance, the method presented in [20]. This procedure does not require inverting matrices, and interestingly, it yields a solution even if some of the equations are inconsistent. In this case, some of the solutions are not going to be integers.
Substituting the values obtained for n i,n into eq. (4.4), we obtain m estimates for ∆t n . If one of these ∆t n estimates were chosen according to some criterion, it could be used to synchronize the g i,n (and eventually c i,n ) coefficients to the input signal by shifting them accordingly. Since the resolution of the ∆t n estimates is not limited by the sampling grid, doing this exactly would require interpolation between the precalculated coefficient samples.
We decided to choose the best sample number offset o n for the coefficient vectors based on the ∆t i,n estimates, thus avoiding interpolation. This is equivalent to choosing the nearest channel of -11 -the filter bank in eq. (2.4). The best sample number offset for channel i can be expressed as where f S is the sampling frequency and · is the nearest integer operator. The sample offset values should be considered as o i,n mod N, where N is the number of channels in the original observer. In order to choose the best o i,n , we propose the following algorithm: • Create channel clusters: using a greedy algorithm, form all possible clusters of channels.
• Assign a score to each cluster: the number of integer n i,n values corresponding to the cluster if all o i,n values are identical, −1 otherwise.
• Find the best estimate: select all clusters having the highest occurring score. For each of these clusters, select the member with the highest amplitude as given by the corresponding w i coefficient. Of these, select the one with the highest w i. The corresponding o i,n becomes o n . Thus, using the notation of figure 5, the state variable update equation takes the form which differs from eq. (2.3) structurally since it accounts for the eventual variation of the g i,n coefficients over time due to the synchronization. The determination of o n through ∆t above causes a significant computational overhead. This overhead can be reduced using the following considerations: • If the value of ∆t is expected to be stable, it is sufficient to calculate it once at a suitably chosen time instant.
• If the value of ∆t is expected to drift slowly, it may be sufficient to calculate it more seldom than the sampling frequency.
• The value for o n obtained with the full calculation may be used as a baseline for subsequent calculations. We evaluate eq. (4.4) for the channel previously selected, reusing the n i value previously obtained. If eq. (4.6) yields a value different from o n , we repeat the whole calculation.

Response to a transient excitation
The transient input signal used for this study, shown in figure 6a, is   with w t 1 = A t 1 e ϕ t 1 = 3 e 1 , w t 2 = A t 2 e ϕ t 2 = 1, w t 3 = A t 3 e ϕ t 3 = 0.7 e 2 , f t 1 = 6 Hz, f t 2 = 3 Hz, f t 3 = 13 Hz and t = t + 0.419 s. The sampling rate is 1 kHz.
The template (see figure 6b) matches the input signal: Figure 6c shows the cross-correlation between the two signals.
The response of the modified Fourier analyzer introduced in section 4.1 to this transient excitation is demonstrated in figure 7. The template in eq. (5.2) is incorporated into the observer. Synchronization between the input signal and the template is ensured by the method based on Diophantine equations presented in section 4.2.  Since the fundamental frequency of the observer is 1 Hz, its dead-beat transients subside in 1 s. This applies to the observer channels used for synchronization and also those used for measurement. As a result, as seen in figure 7, the latter exhibit a coarse transient phase while the former are in transient, then the latter settle in a fine transient phase of identical duration.

Response to a non-ideal excitation
For a series of further tests, we used the following input signal: .237 s and all other parameters are the same as in section 5.1. Once more, the sampling rate is 1 kHz.
The w t i components are intended to represent the useful signal, while the d t i components correspond to an additive in-band disturbance switched on at t = 3 s. As a result, this part of the signal still contains the same frequency components, but with different amplitude and phase relations than those prescribed by the template. This means that ideally, the detection algorithm should be able to identify this situation as problematic. Moreover, the test signal is contaminated with additive white noise, SNR = 20 dB. The resulting signal is shown in figure 8. Figure 9 shows the output of an ordinary matched filter with a template corresponding to eq. (5.2) applied to the test signal. After the initial transient, numerous positive peaks are present in the signal. The in-band disturbance does not significantly change the amplitudes of the peaks, although at this point, the signal is no longer a good match for the template.
The operation of the modified Fourier analyzer described in section 4.1, complemented by the synchronization method suggested in section 2.4, is demonstrated in figure 10. By operating the recursive LTV matched filter presented in section 2.2 with the same synchronization method, we would obtain the behavior seen in Re {Ch1}. After the transients caused by turning on the disturbance have settled, this channel shows a slight increase, however, the additional channels indicate that the input signal no longer matches the expected signal. Figure 11 shows the behavior of the same structure, but using the synchronization approach described in section 4.2. In this case, the disturbance has practically no influence on the value of the channel used for detection, but it is clearly revealed in the orthogonal channels.

Conclusions
In this paper, we presented a time variant filter for signal presence detection. As a result, its response can exhibit more advantageous characteristics than that of the initial matched filter. We described ways to implement the structure, along with two solutions to perform the synchronization required for its operation. One of the suggested methods gives rise to orthogonal cross-correlation channels, which allow improving the confidence of signal presence detection based on a matched filter.

Future work
It might be interesting to analyze the information content of the orthogonal channels more systematically. They may contain additional quantitative information related to the structure of the input signal.   (b) Time lag between the input signal and the template, obtained using the method relying on Diophantine equations. This value is used to choose the g i,n coefficient values to use in the next iteration. It is interesting to note that after the disturbance is switched on, this method converges back to the correct time lag value. This is due to the fact that there are no disturbing components in the channel with the highest design amplitude (see the description of the algorithm in section 4.2). Figure 11. Response of the modified Fourier analyzer from section 4.1 to the input signal described in section 5.2. The observer is synchronized using the solution suggested in section 4.2.