A Modulated Wideband Converter Model Based on Linear Algebra and Its Application to Fast Calibration

Abstract

In the context of cognitive radio, smart cities and Internet-of-Things, the need for advanced radio spectrum monitoring becomes crucial. However, surveillance of a wide frequency band without using extremely expensive high sampling rate devices is a challenging task. The recent development of compressed sampling approaches offers a promising solution to these problems. In this context, the Modulated Wideband Converter (MWC), a blind sub-Nyquist sampling system, is probably the most realistic approach and was successfully validated in real-world conditions. The MWC can be realized with existing analog components, and there exist calibration methods that are able to integrate the imperfections of the mixers, filters and ADCs, hence allowing its use in the real world. The MWC underlying model is based on signal processing concepts such as filtering, modulation, Fourier series decomposition, oversampling and undersampling, spectrum aliasing, and so on, as well as in-flow data processing. In this paper, we develop an MWC model that is entirely based on linear algebra, matrix theory and block processing. We show that this approach has many interests: straightforward translation of mathematical equations into simple and efficient software programming, suppression of some constraints of the initial model, and providing a basis for the development of an extremely fast system calibration method. With a typical MWC acquisition device, we obtained a speed-up of the calibration computation time by a factor greater than 20 compared with a previous implementation.

1. Introduction

Digital wireless radio signals are often composed of a small number of narrow-band transmissions spread across a wide spectrum range. For example, the Internet-of-Things (IoT) communications have recently emerged in contexts such as smart cities. Cognitive radios are able to manage the spectrum dynamically but require advanced sensing techniques for spectrum monitoring.

Basically, spectrum monitoring is based on the Shannon–Nyquist sampling theorem [1,2]. This theorem states that the signal must be sampled at a rate greater than its Nyquist frequency, which is twice its frequency band. However, when we have to monitor a large frequency band, this requirement can exceed the capabilities of existing Analog to Digital Converters (ADC) or require very expensive components. Furthermore, sampling at a very high rate may require huge storage capacities to store the digital samples.

Recently, new approaches have been proposed allowing sampling at sub-Nyquist rates. These approaches, known as compressed sensing, or compressive sampling [3], have emerged as a promising framework for signal acquisition in difficult applications, such as monitoring a wideband spectrum [4]. The basic idea of compressed sampling is to take advantage of the fact that a signal that has a sparse representation on a given basis can theoretically be recovered from a small set of linear measurements [5,6]. The price to pay is the need to develop sophisticated signal processing algorithms to reconstruct the signal from this small set of measurements, these algorithms being much more complex than the usual demodulators.

A great deal of the theoretical aspects of compressed sampling has been addressed in the literature. For example, many studies have been proposed in relation to the design of the measurement scheme as in [7,8]. However, few studies have considered the practical limitations of compressed acquisition. In fact, designing measurement schemes and applying them to practical acquisition systems remains a significant challenge.

In this context, the Modulated Wideband Converter (MWC) has been proposed as an efficient system for real-world compressed sampling [9]. The MWC does sub-Nyquist sampling without prior information about the spectral support of the transmitters present in the monitored wideband. It can be realized with existing devices [10] and has been successfully tested on real-world problems, including surveillance of a wideband spectrum [11].

A few real-world compressed sensing acquisition systems have already been proposed [11,12,13,14,15,16,17,18]. An analog circuit board with discrete components was designed by us to prototype compressed sensing based on the MWC principle [19].

A necessary step when using MWC in real-world conditions is the calibration of the acquisition system. Indeed, analog components are never ideal, especially when fed with wideband signals. Then, using a purely theoretical model leads to extremely poor performances of signal reconstruction. An efficient calibration method, which is considered a reference, has been proposed in [13]. It consists of estimating the sensing matrix column after column by injecting sine waves at a specific frequency and recording the corresponding output signals. The procedure is repeated by changing the input frequency until all columns are estimated. Some researchers have exploited this work to calibrate their systems or to propose variants of the calibration algorithms [18,20,21,22,23]. While this procedure is efficient, it can be time-consuming because the number of columns to estimate is usually at least a few dozen. That is why a few authors [24,25], including us [19], have recently proposed alternative calibration algorithms requiring only one input signal.

The MWC theoretical background is signal processing theory (filters, modulation, Fourier series, sampling theory, spectrum aliasing, etc.). Most signal processing theoretical tools are asymptotic. However, when signals are processed in real-world conditions, they are always finite; thus, block processing and purely matrix-based algorithms may be more natural and efficient.

Moreover, a quick look at the literature shows that most people use Matrix-based programming tools, such as Octave or Matlab, for signal processing in this context, but without really exploiting the power of these tools. To take full advantage of the power of Matrix-oriented software, it would then be preferable to process data by blocks instead of in-flow. It is, therefore, interesting to view the whole MWC acquisition and reconstruction in terms of block processing. The most natural framework to achieve this objective is matrix theory and linear algebra.

In this paper, we elaborate on an MWC model using linear algebra only (without any signal processing theory). While this approach will probably appear less intuitive than the approach based on signal processing, because most people are less familiar with linear algebra than with signal processing, it has strong advantages:

• Once the model is established, programming it becomes extremely simple, straightforward and efficient.
• Furthermore, computational complexity is significantly reduced.
• The border effects are implicitly taken into account in the model. Indeed, using a signal processing model, people have to deal with the fact that the signals processed in the real world are not infinite, while when using a linear algebra model, the finite nature of the data is implicitly taken into account and the mathematical equations are exact and not approximate.
• In the original version of the MWC, the number of physical branches is increased by a factor q, which must be an odd integer due to signal processing considerations. An interesting aspect of the linear algebra model is that it allows even integers for q.

The main contributions of this paper are:

• The development of a pure linear-algebra model of the MWC. Despite the fact that establishing this model is rather hard because it requires non-trivial matrix manipulations, once established, it is extremely simple and allows programming MWC-related software, such as calibration, in a very fast, compact and efficient way.
• Its application to the development of a very fast calibration method. With typical choices of parameters, the calibration is more than 20 times faster than our previous method (this previous method being itself very fast compared to a reference method because it required only one calibration signal instead of dozens of sinusoidal signals in the reference method).

The remainder of this paper is organized as follows: Section 2 provides the main mathematical tools used in the paper. Then, Section 3 presents an overview of our hardware acquisition board and the MWC principle. Section 4 establishes a system model based on linear algebra, and an equivalent model, useful for signal reconstruction and system calibration, is then derived in Section 5. In Section 6, we show how this model allows us to considerably improve a calibration method that we proposed previously, leading to speeding up the process by a factor greater than 20. Then, some experimental results are shown in Section 7.

2. Mathematical Background

2.1. Notations

Unless otherwise stated, lowercase symbols denote row vectors (e.g., x, p), uppercase symbols denote matrices (e.g., C, Z), $\overline{x}$ stands for the DFT (Discrete Fourier Transform) of x. The symbols $N,K,L,a,b$ will be used to denote the size of vectors or matrices.

We will denote $D_x$ as the square diagonal matrix whose diagonal is vector x.

The vectorization of a $K\times L$ matrix Q, denoted $vec\left(Q\right)$, is the $1\times KL$ row vector obtained by reading the matrix row after row, from top to bottom:

$$vec\left(Q\right)=\left(\begin{array}{ccc}{q}_{11}\cdots q_{1L}& q_{21}\cdots q_{2L}& \cdots q_{K1}\cdots q_{KL}\end{array}\right)$$

(1)

$M^{*}$ stands for the Hermitian transpose of matrix $M.$

$I_K$ stands for the $K\times K$ identity matrix.

The nearest lower or equal integer will be noted $\lfloor \rfloor$ and the nearest greater or equal integer $\lceil \rceil$.

2.2. Circulant Matrices

Let x be a $1\times N$ row vector. A circulant matrix $C_x$ is a square matrix whose first row is x and each next row is a circular shift one element to the right of the preceding row. That is:

$$C_x=\left(\begin{array}{ccccc}{x}_o& x_1& x_2& \cdots & x_{N-1}\\ x_{N-1}& x_0& x_1& \cdots & x_{N-2}\\ x_{N-2}& x_{N-1}& x_0& \cdots & x_{N-3}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_1& x_2& x_3& \cdots & x_0\end{array}\right)$$

(2)

It is convenient to define the cyclic permutation matrix as the $N\times N$ matrix below:

$$J_N=\left(\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& & 0\\ ⋮& & & \ddots \\ 0& 0& 0& & 1\\ 1& 0& 0& \cdots & 0\end{array}\right)$$

(3)

Then, $C_x$ is a polynomial in $J_N$:

$$C_x=\sum _{n=0}^{N-1}{x}_n{J}_N^n$$

(4)

The effect of multiplication of a matrix M by $J_N$ is as follows. The rows of $M{J}_N$ are the rows of M circularly shifted one element to the right. The columns of $J_{N}M$ are the columns of M circularly shifted one element to the top.

Matrices $J_N$ and $C_x$ commute:

$$J_N{C}_x=C_x{J}_N$$

(5)

because

$$\begin{array}{cc}\hfill J_N{C}_x& =J_N\left(\sum _{n=0}^{N-1}{x}_n{J}_N^n\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & =\sum _{n=0}^{N-1}{x}_n{J}_N^{n+1}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & =\left(\sum _{n=0}^{N-1}{x}_n{J}_N^n\right)J_N\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & =C_x{J}_N\hfill \end{array}$$

(9)

2.3. Discrete Fourier Transform (DFT)

Let us note $\omega$ the $N^{th}$ square root of unity below:

$$\omega =exp\left(-i\frac{2\pi }{N}\right)$$

(10)

The DFT matrix $F_N$ is an $N\times N$ square symmetric matrix whose element at row l column k is ${\omega }^{lk}$ (assuming row 0 is the first row, and column 0 the first column):

$$F_N=\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{N-1}\\ 1& {\omega }^2& {\omega }^4& \cdots & {\omega }^{2(N-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{N-1}& {\omega }^{2(N-1)}& \cdots & {\omega }^{{(N-1)}^2}\end{array}\right)$$

(11)

The inverse DFT matrix is

$$F_N^{-1}=\frac{1}{N}{F}_N^{*}$$

(12)

The DFT of a vector x is

$$\overline{x}=x{F}_N$$

(13)

and the inverse DFT (IDFT) is given by $\overline{x}{F}_N^{-1}$.

A circulant matrix is diagonalized by the DFT matrix. That is

$$C_x=F_N{D}_{\overline{x}}{F}_N^{-1}$$

(14)

It follows that the elements of $\overline{x}$ are the eigenvalues of $C_x$ and the columns of $F_N^{-1}$ are the eigenvectors.

We also have

$$D_{\overline{x}}=F_N^{-1}{C}_x{F}_N$$

(15)

and

$$\frac{1}{N}{C}_{\overline{x}}=F_N^{-1}{D}_x{F}_N.$$

(16)

2.4. Kronecker Product

The Kronecker product is a bilinear matrix operation, denoted by ⊗. If A is a $K\times L$ matrix and B is a $M\times N$ matrix, it produces the $KM\times LN$ block matrix C below:

$$C=A\otimes B=\left(\begin{array}{ccc}{a}_{11}B& \cdots & a_{1L}B\\ ⋮& \ddots & ⋮\\ a_{K1}B& \cdots & a_{KL}B\end{array}\right)$$

(17)

A useful property about the inverse is:

$${\left(A\otimes B\right)}^{-1}=A^{-1}\otimes B^{-1}$$

(18)

Assuming the sizes are such that one can form the matrix products $AC$ and $BD$, an interesting property, known as the mixed-product property, is:

$$(A\otimes B)(C\otimes D)=\left(AC\right)\otimes \left(BD\right)$$

(19)

The product is not commutative, but there exist permutation matrices (shuffle matrices) such that if A is an $a\times a$ square matrix and B a $b\times b$ square matrix, then [26]:

$$A\otimes B=P_{a,b}(B\otimes A)P_{b,a}$$

(20)

Matrix $P_{a,b}$ represents the permutation obtained when elements, written row by row in an $a\times b$ matrix, are read column by column. For instance, set $a=2$ and $b=3$, and write the elements $1,2,3,4,5,6$ row by row in a $a\times b$ matrix

$$\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\end{array}\right)$$

(21)

Reading column by column, we obtain $1,4,2,5,3,6$. The permutation matrix is then the matrix such that:

$$\left(\begin{array}{cccccc}1& 4& 2& 5& 3& 6\end{array}\right)=\left(\begin{array}{cccccc}1& 2& 3& 4& 5& 6\end{array}\right)P_{2,3}$$

(22)

That is:

$$P_{2,3}=\left(\begin{array}{c}1\\ & & 1\\ & & & & 1\\ & 1\\ & & & 1\\ & & & & & 1\end{array}\right)$$

(23)

If $N=KL$, an interesting property with the permutation matrix defined in (3) is

$$J_N^K=J_L\otimes I_K.$$

(24)

2.5. General Radix Identity

If N is a composite number, i.e., $N=KL$, then [26]:

$$F_N=(F_K\otimes I_L)T_{K,L}(I_K\otimes F_L)P_{K,L}$$

(25)

where $T_{K,L}$ is a diagonal matrix (twiddle matrix) and $P_{K,L}$ a permutation matrix (shuffle matrix defined in Section 2.4).

The twiddle matrix $T_{K,L}$ is an $N\times N$ diagonal matrix, the diagonal of which is $vec\left(Q_{K,L}\right)$ with $\omega$ defined in Equation (10) and

$$Q_{K,L}=\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{L-1}\\ 1& {\omega }^2& {\omega }^4& \cdots & {\omega }^{2(L-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{K-1}& {\omega }^{2(K-1)}& \cdots & {\omega }^{(K-1)(L-1)}\end{array}\right)$$

(26)

For instance, with $K=2$ and $L=3$, we have

$$Q_{2,3}=\left(\begin{array}{ccc}1& 1& 1\\ 1& e^{-i\pi /3}& e^{-2i\pi /3}\end{array}\right)$$

(27)

and the diagonal of $T_{2,3}$ is

$$diag\left(T_{2,3}\right)=\left(\begin{array}{cccccc}1& 1& 1& 1& e^{-i\pi /3}& e^{-2i\pi /3}\end{array}\right)$$

(28)

Let us note ${\theta }_K$ and $1_K$ as the $(1\times K)$ vectors below

$${\theta }_K=\left[\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right]$$

(29)

$$1_K=\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]$$

(30)

Note that for any $1\times L$ vector p, we have

$$(1_K\otimes p)T_{K,L}=(1_K\otimes p)$$

(31)

because only the first L elements of $1_K\otimes p$ are non null, and the L first elements of $T_{K,L}$ are ones.

Note also that

$$(1_K\otimes p)P_{K,L}=p\otimes 1_K$$

(32)

when elements of $1_K\otimes p$ are written row by row in a $K\times L$ matrix, the elements of p go on the first row and the $K-1$ next rows are null. Then, when this matrix is read column by column, we obtain elements of p separated by $K-1$ zeroes, that is, $p\otimes 1_K$

Similarly, it is easy to check that

$$T_{a,b}^{-1}(I_a\otimes 1_b^T)=I_a\otimes 1_b^T$$

(33)

and

$$P_{a,b}^{-1}(I_a\otimes {\theta }_b^T)={\theta }_b^T\otimes I_a.$$

(34)

2.6. Selection Matrix

Let us define the selection matrix $S_{N,K}^{\left(r\right)}$ as the $N\times K$ matrix below:

$$S_{N,K}^{\left(r\right)}=\left(\begin{array}{c}{0}_{r\times K}\\ I_K\\ 0_{(N-K-r)\times K}\end{array}\right)$$

(35)

If $x=\left(x_0\cdots x_{N-1}\right)$ is a $1\times N$ vector, then $y=x{S}_{N,K}^{\left(r\right)}$ is the $1\times K$ vector below:

$$y=\left(x_r\cdots x_{r+K-1}\right)$$

(36)

We consider the indices modulo N, so r may be negative.

2.7. Moore–Penrose Pseudo-Inverse

Let us consider a rectangular matrix Z whose size is $L\times K$ with $L\le K$. The Moore–Penrose pseudo-inverse [27] of Z, denoted $Z^{+}$, is a $K\times L$ matrix, which generalizes the concept of inverse and, among other interesting properties, provides a mean to compute a least squares solution to a system of linear equations that lacks an exact solution. The pseudo-inverse is defined and unique for all complex matrices. It is usually computed using the singular value decomposition (SVD).

Let us note the SVD of Z as [28]:

$$Z=US{V}^{*}$$

(37)

where U is a $L\times L$ unitary matrix (i.e., $U{U}^{*}=U^{*}U=I$), V is a $K\times L$ matrix with orthonormal columns (i.e., $V^{*}V=I$) and S is a diagonal matrix whose diagonal elements are the singular values (non-negative real numbers, ranked by decreasing order). The SVD exists for all complex matrices.

Here we consider a version of the $SVD$ usually called “thin-SVD”, which is a compact version of the more general $SVD$ decomposition (in which matrices S and V are larger), because this compact version is sufficient for the purpose of computing the pseudo-inverse. The computational cost of computing the thin-SVD is $6K{L}^2+20L^3$ flops ([28] p. 254). Note that for complex matrices, it is usual to redefine the floating point operation (flop) in order to count only one flop for the product of two complex numbers, while in reality it requires four real multiplications. Since only ratios between the computational costs of algorithms are of interest, applying this does not change the result.

The pseudo-inverse is given by:

$$Z^{+}=V{S}^{+}{U}^{*}$$

(38)

where $S^{+}$ is the pseudo-inverse of S. It is a diagonal matrix in which the diagonal contains the inverses of the singular values of S, which are above a small tolerance value, and 0 elsewhere.

The cost of the inversion plus the computation of the matrix product is $2L+K{L}^2\simeq K{L}^2$.

Overall, the cost of computing the pseudo-inverse is $7K{L}^2+20L^3$.

3. Acquisition Device and System Parameters

The MWC is a compressed sampling device that samples a signal $x\left(t\right)$ at a sampling frequency $F_s$ significantly lower than its Nyquist frequency $F_{nyq}$. The input signal is assumed sparse in the frequency domain. From the outputs of this acquisition device, one can reconstruct the input signal using a compressed sensing algorithm, such as Orthogonal Matching Pursuit (OMP) [29].

The principle of the MWC is shown in Figure 1:

• The input signal $x\left(t\right)$ is multiplied (using a mixer) by a scrambling signal $s\left(t\right)$.
• The resulting signal $v\left(t\right)$ goes through a low-pass filter whose impulse response is $h\left(t\right)$.
• Then, the filter output $w\left(t\right)$ is sampled by an Analog to Digital Converter (ADC), providing the output samples $y\left[n\right]$.

The scrambler $s\left(t\right)$ is a periodic signal: it is a basic waveform $p\left(t\right)$ repeated $F_p$ times per second. The analog waveform $p\left(t\right)$ itself is generated at sampling frequency $F_{nyq}$ from an L samples digital sequence, which is usually a pseudo-random sequence. Consequently, we have $F_p=F_{nyq}/L$.

The performance of the system can be enhanced by using M parallel branches with different scrambling signals. However, since generalization to M branches is trivial, we will restrict the discussions below to one branch.

The digital outputs $y\left[n\right]$ are provided at $F_s$ samples per second.

In previous practical realizations, in order to reduce aliasing, the ADC output samples go through a digital filter, which provides properly filtered samples at a frequency $F_{ss}$ lower than $F_s$. In the original MWC model, $F_{ss}$ is an odd multiple of $F_p$, that is $F_{ss}=q{F}_p$ with q an odd integer. In this paper, since the linear algebra model allows a less constrained post-processing, this digital filter is not required and q is not necessarily odd. Indeed, we will see that the linear algebra model also allows even values of q.

When designing an actual acquisition device, we have to choose some parameters:

• The sampling frequency $F_{nyq}$ of the scramblers, which will impact the Nyquist frequency of acceptable input signals (i.e., input signal maximum frequency must remain under $F_{nyq}/2$).
• The sampling frequency $F_s$ of the ADC, which should be significantly lower than $F_{nyq}$ (otherwise the system would have no interest compared to direct sampling at $F_{nyq}$). This frequency determines the subsampling factor $b=F_{nyq}/F_s$.
• The length L of the scrambler periodic pattern. This parameter determines the frequency of repetition $F_p=F_{nyq}/L$ of the scrambling pattern.

Figure 2 illustrates, in a very simplified case ($L=5$), an example of scrambling signal. It is formed from a length-L binary pseudo-random sequence, which is repeated. In the time domain, the duration of a binary symbol is $1/F_{nyq}$, and the period of the signal is $1/F_p$ (where $F_{nyq}$ and $F_p$ are the frequencies defined above).

Figure 3 illustrates the low-pass filter response. The scrambled signal, which occupies a frequency band of width $F_{nyq}$, goes through a low-pass filter. The filter output is sampled at a frequency $F_s$, high enough to avoid aliasing.

The scrambler and the ADC are controlled by a common central clock to avoid synchronization problems.

Reconstruction of the input signal, and calibration of the system, are based on the information provided by a block of a output samples. In order to avoid unnecessary mathematical complications, the value of a is chosen such that it corresponds to an integer number K of scrambling patterns, then $a=KL/b$. This output block then corresponds to $N=KL$ scrambler samples (and also to N input samples if the input signal were sampled at $F_{nyq}$). The size of the block determines the frequency resolution $F_s/a=F_{nyq}/N$.

For our experiments on real-world data, we designed a 4-channel MWC analog board (Figure 4), which was described in more detail in a previous paper [19]. The scramblers are sampled at $F_{nyq}=1$ GHz and their length is $L=96$. Therefore, their repetition frequency is $F_p=F_{nyq}/L=10.41667$ MHz. The device is then able to monitor a wideband spectrum of 1 GHz.

The prototype includes $M=4$ physical channels. Each channel features an M1-0008 mixer from MArki^©, and an SXLP-36+ low-pass filter from Mini-Circuits^©. The filter cut-off frequency is 40 MHz (at -3 dB). The SXLP-36+ filter was chosen because it is quite close to the ideal low-pass filter. Indeed, it has a sharp cut-off, linear phase and flat band (attenuation < 0.5 dB) in the frequency range (0–36) MHz. The ADC sampling frequency is $F_s=10F_p=104.1667MHz$ (at $F_s/2$ the filter attenuation is more than 30 dB); therefore, the downsampling factor is $b=9.6$.

Table 1. Parameters of our MWC prototype.

Symbol	Meaning	Value
M	Number of physical channels	4
L	Length of scramblers	96
$F_{nyq}$	Sampling frequency of scramblers = bandwidth to monitor	1 GHz
$F_s$	Sampling frequency of physical ADC	104.1667 MHz
$b=F_{nyq}/F_s$	Physical subsampling factor	9.6
$F_p=F_{nyq}/L$	Repetition frequency of scramblers	10.41667 MHz

sums up the main parameters.

The frequency response of the low-pass filter implemented on our acquisition board is shown in Figure 5 and its phase i n Figure 6. Details on filter calibration can be found in one of our previous papers [30].

$F_{nyq}$ being the Nyquist frequency of the input signal, we can consider a digital equivalent model at $F_{nyq}$ without loss of information. Furthermore, as previously mentioned, since calibration and signal processing are always performed on a limited amount of data in real-world applications, we can consider an input block of N samples (at $F_{nyq}$).

Modern implementations of the FFT [31] contain a special code to handle splittings not only of size 2 but also of sizes 3 (and sometimes 5 and 7). Therefore, for the efficiency of the FFT, we will preferably choose block sizes whose prime factors belong to $\{2,3,5,7\}$. In our experiments, we have taken $K=448$, $N=KL=43,008=2^{11}\times 3\times 7$ and $a=4480=2^9\times 5\times 7$. The frequency resolution is then $F_{nyq}/N\simeq 23$ kHz, which is far sufficient unless we would like to detect extreme narrow-band transmitters.

4. System Model and Matrix Representation

4.1. System Equations in the Time Domain

Let us note:

• x the vector representing the input signal.
• s the vector representing the scrambling signal.
• v the vector representing the scrambler output.
• h the vector representing the low-pass filter impulse response.
• w the vector representing the low-pass filter output.
• y the $1\times a$ vector containing the digital output samples (at $F_s$).

All of these vectors, except for y, are $(1\times N)$ vectors and represent the signals at $F_{nyq}$ samples per second. In Figure 7, we show the links between these vectors. In the time domain (top of the figure), the signals are represented by vectors. Symbol * stands for cyclic convolution. These vectors can be transposed in the frequency domain using a multiplication by matrix $F_N$ or $F_a$. A post-processing, described later, is then performed in the frequency domain. The post-processing outputs q vectors ${\tilde{y}}_n$ of size $1\times K$.

In the figure, we have used different symbols for down-sampling because the operation in the time and frequency domains are different. For instance, when b is an integer, down-sampling in the time domain consists of picking one sample out of b while its equivalent in the frequency domain is a multiplication of the down-sampling matrix $\mathsf{\Xi }$, which will be defined later.

Notations used below have already been defined in Section 2. Since the system is linear, in the time domain, we have

$$y=xB$$

(39)

where B is an $N\times a$ matrix. The structure of B can be easily computed from the system model (Figure 7):

$$B=D_s{C}_h\left(I_a\otimes 1_b^T\right)$$

(40)

Indeed, the scrambler output is given by:

$$v=x{D}_s$$

(41)

The filter output is:

$$w=v{C}_h$$

(42)

For the moment, let us consider that b is an integer (we will see later that this is not a requirement). In that case, down-sampling consists of picking one sample out of b in w. Mathematically, that is:

$$y=w\left(I_a\otimes 1_b^T\right)$$

(43)

Otherwise, down-sampling can be modeled using an interpolation matrix. However, we will not detail this because only the equations in the frequency domain will be useful for our purpose. We will see later that in the frequency domain, thanks to the presence of a low pass filter, b being an integer is not a requirement anymore.

4.2. System Equations in the Frequency Domain

Multiplying Equation (39) by $F_a$ and inserting the identity $F_N{F}_N^{-1}$ where appropriate, we obtain:

$$y{F}_a=\left(x{F}_N\right)\left(F_N^{-1}B{F}_a\right)$$

(44)

That is

$$\overline{y}=\overline{x}A$$

(45)

where A is an $N\times a$ matrix.

$$A=F_N^{-1}B{F}_a$$

(46)

The structure of A can be detailed further. Replacing B with its expression (Equation (40)) and inserting the identity $F_N{F}_N^{-1}$ where it is appropriate, we obtain:

$$\begin{array}{cc}\hfill A& =\left(F_N^{-1}{D}_s{F}_N\right)\left(F_N^{-1}{C}_h{F}_N\right)\left(F_N^{-1}\left(I_a\otimes 1_b^T\right)F_a\right)\hfill \end{array}$$

(47)

Then, using Equations (15) and (16) we obtain:

$$A=\frac{1}{N}{C}_{\overline{s}}{D}_{\overline{h}}\mathsf{\Xi }$$

(48)

As proved in the appendix (see Appendix A.1), the frequency-domain down-sampling matrix $\mathsf{\Xi }$ is:

$$\mathsf{\Xi }=\frac{1}{b}{\theta }_b^T\otimes I_a$$

(49)

That is:

$$\mathsf{\Xi }=\frac{1}{b}\left(\begin{array}{c}{I}_a\\ ⋮\\ I_a\end{array}\right)$$

(50)

where sub-matrix $I_a$ is repeated b times. Here we remind that h is a low-pass filter. Since the ADC sampling frequency is $F_s$, we assume that the elements of $\overline{h}$, which correspond to frequencies outside the interval $]-F_s/2,F_s/2$ [are almost null. Since $\overline{h}$ contains N elements, the frequency resolution is $F_{nyq}/N$, so $F_s/2$ corresponds to index $(F_s/2)/(F_{nyq}/N)=N/\left(2b\right)$, that is, $a/2$. Let us note

$$c=⌊a/2⌋$$

(51)

and

$$\delta =amod2$$

(52)

Therefore, the elements of $\overline{h}$ are almost null for indices outside the interval $[-c,c+\delta ]$ (the indices are considered modulo N). Hence, we can redefine $\mathsf{\Xi }$ as the $N\times a$ matrix below:

$$\mathsf{\Xi }=\frac{1}{b}\left(\begin{array}{cc}{I}_{c+\delta }& 0\\ 0& 0\\ 0& I_c\end{array}\right)$$

(53)

without changing the product $D_{\overline{h}}\mathsf{\Xi }$. Here, the zeros stand for null sub-matrices. We see that, thanks to the low-pass pass filter which leads to this structure of $\mathsf{\Xi }$, it is not required any more that b is an integer (this requirement was only due to the need for an integer number of occurrences of $I_a$ in Equation (50)).

Finally, let us define the ($1\times a)$ vector $\tilde{h}$:

$$\tilde{h}=\left[\begin{array}{cccccc}{\overline{h}}_0& \cdots & {\overline{h}}_{c+\delta -1}& {\overline{h}}_{-c}& \cdots & {\overline{h}}_{-1}\end{array}\right]$$

(54)

where the indices are modulo N. We then have:

$$D_{\overline{h}}\mathsf{\Xi }=\mathsf{\Xi }{D}_{\tilde{h}}$$

(55)

and the expression of matrix A becomes:

$$A=\frac{1}{N}{C}_{\overline{s}}\mathsf{\Xi }{D}_{\tilde{h}}$$

(56)

4.3. Unconstrained System Equations in the Time Domain

Now we can go back to the time domain to obtain a matrix B, which does not require b being an integer. We have:

$$\begin{array}{cc}\hfill y& =\overline{y}{F}_a^{-1}\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill & =\overline{x}A{F}_a^{-1}\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & =\left(\overline{x}{F}_N^{-1}\right)\left(F_{N}A{F}_a^{-1}\right)\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill & =xB\hfill \end{array}$$

(60)

where

$$B=F_{N}A{F}_a^{-1}$$

(61)

4.4. Fast Simulation of the Acquisition System

A first interest of the linear algebra model is that it makes the design of a fast simulator obvious. Indeed, multiplication by a diagonal matrix D is efficiently implemented as an element-by-element vectors product, and multiplication by a Fourier matrix F (or its inverse) is efficiently implemented by Fast Fourier Transform (FFT). On the contrary, multiplications by circulant matrices C should be avoided because of their computational cost. Then, the method to design a fast simulator is to insert identities $F{F}^{-1}$ or $F^{-1}F$ where it is appropriate in order to suppress the circulant matrices. For instance, we have:

$$\begin{array}{cc}\hfill y& =xB\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill & =x{F}_N\frac{1}{N}{C}_{\overline{s}}\mathsf{\Xi }{D}_{\tilde{h}}{F}_a^{-1}\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill & =x\left(F_N\frac{1}{N}{C}_{\overline{s}}{F}_N^{-1}\right)F_N\mathsf{\Xi }{D}_{\tilde{h}}{F}_a^{-1}\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & =x{D}_s{F}_N\mathsf{\Xi }{D}_{\tilde{h}}{F}_a^{-1}\hfill \end{array}$$

(65)

using Equation (16). Here we have only fast operations, as shown in Figure 8.

5. Equivalent Model and Post-Processing

5.1. Equivalent Model

Until now, we have not taken advantage of the periodicity of the scrambler. This opens the way to an equivalent model with interesting properties.

The scrambler is a $(1\times N)$ vector s, which contains K replica of a basic waveform represented by a $(1\times L)$ vector p. Then, the scrambler can be written:

$$s={\theta }_K\otimes p$$

(66)

and we have (see proof in Appendix A.2)

$$\overline{s}=K\overline{p}\otimes 1_K$$

(67)

It follows that $\overline{s}$ is sparse (only one element out of K is non-zero). It will be easier to take benefit of the sparsity of $\overline{s}$ if we permute $\overline{x}$ and $\overline{s}$ in the expression of $\overline{y}$:

$$\begin{array}{cc}\hfill \overline{y}& =\frac{1}{N}\overline{x}{C}_{\overline{s}}\mathsf{\Xi }{D}_{\tilde{h}}\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill & =\frac{1}{N}\overline{s}{C}_{\overline{x}}\mathsf{\Xi }{D}_{\tilde{h}}\hfill \end{array}$$

(69)

The proof is trivial: since the multiplication is commutative, we can permute x and s (see Figure 7); therefore, we can also permute $\overline{x}$ and $\overline{s}$.

Let us define the $L\times N$ matrix $C_{\overline{x}}^{\left(K\right)}$ obtained by picking one row out of K in $C_{\overline{x}}$. That is:

$$C_{\overline{x}}^{\left(K\right)}=\left(I_L\otimes 1_K\right)C_{\overline{x}}$$

(70)

More explicitly, that is:

$$C_{\overline{x}}^{\left(K\right)}=\left(\begin{array}{ccc}{\overline{x}}_0& \cdots & {\overline{x}}_{N-1}\\ {\overline{x}}_{-K}& \cdots & {\overline{x}}_{N-K-1}\\ ⋮& & ⋮\\ {\overline{x}}_{-(L-1)K}& \cdots & {\overline{x}}_{K-1}\end{array}\right)$$

(71)

where the indices are considered modulo N.

Let us denote

$$\tilde{p}=\frac{1}{L}\overline{p}$$

(72)

Using Equation (67), the mixer output becomes:

$$\begin{array}{cc}\hfill \frac{1}{N}\overline{s}{C}_{\overline{x}}& =\frac{1}{N}\left(K\overline{p}\otimes 1_K\right)C_{\overline{x}}\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill & =\frac{1}{L}\left(\left(\overline{p}{I}_L\right)\otimes 1_K\right)C_{\overline{x}}\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill & =\frac{1}{L}\overline{p}\left(I_L\otimes 1_K\right)C_{\overline{x}}\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill & =\tilde{p}{C}_{\overline{x}}^{\left(K\right)}\hfill \end{array}$$

(76)

An interesting property of matrix $C_{\overline{x}}^{\left(K\right)}$, which will be exploited later, is (see proof in Appendix A.3):

$$C_{\overline{x}}^{\left(K\right)}{J}_N^K=J_L{C}_{\overline{x}}^{\left(K\right)}$$

(77)

Finally, let us define

$$\tilde{y}=\overline{y}{D}_{1/\tilde{h}}$$

(78)

We then have:

$$\begin{array}{cc}\hfill \tilde{y}& =\frac{1}{N}\overline{s}{C}_{\overline{x}}\mathsf{\Xi }\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill & =\tilde{p}{C}_{\overline{x}}^{\left(K\right)}\mathsf{\Xi }\hfill \end{array}$$

(80)

5.2. Post-Processing

The post-processing extracts frequency blocks of K samples from $\tilde{y}$.

Using definition (35), let us note $S_n$ the $a\times K$ selection matrix below

$$S_n=b{S}_{a,K}^{(r+nK)}$$

(81)

and $R_n$ the $N\times K$ selection matrix below

$$R_n=S_{N,K}^{(r+nK)}$$

(82)

Denoting ${\tilde{y}}_n$ a $1\times K$ vector representing the selected data, we have:

$${\tilde{y}}_n=\tilde{y}{S}_n$$

(83)

${\tilde{y}}_n$ contains the elements of $\tilde{y}$ whose indices (modulo a) are in the interval $\Phi =[r+nK,r+nK+K-1]$.

The indices are considered modulo a, so r may be negative. We will consider that $\Phi \subset [-c,c+\delta ]$, so

$$\begin{array}{cc}\hfill \mathsf{\Xi }{S}_n& =\mathsf{\Xi }b{S}_{a,K}^{(r+nK)}\hfill \end{array}$$

(84)

$$\begin{array}{cc}\hfill & =S_{N,K}^{(r+nK)}\hfill \end{array}$$

(85)

$$\begin{array}{cc}\hfill & =R_n\hfill \end{array}$$

(86)

We can note that:

$$R_n=J_N^{-nK}{R}_0$$

(87)

This is a matrix similar to $R_0$ but with sub-matrix $I_K$ circularly shifted $nK$ positions downwards. We can note that we have also:

$$R_{n+1}=J_N^{-K}{R}_n$$

(88)

Eventually, using Equations (73) and (77) we have:

$$\begin{array}{cc}\hfill {\tilde{y}}_n& =\tilde{p}{C}_{\overline{x}}^{\left(K\right)}\mathsf{\Xi }{S}_n\hfill \end{array}$$

(89)

$$\begin{array}{cc}\hfill & =\tilde{p}{C}_{\overline{x}}^{\left(K\right)}{R}_n\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill & =\tilde{p}{C}_{\overline{x}}^{\left(K\right)}{J}_N^{-nK}{R}_0\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill & =\left(\tilde{p}{J}_L^{-n}\right)\left(C_{\overline{x}}^{\left(K\right)}{R}_0\right)\hfill \end{array}$$

(92)

$$\begin{array}{cc}\hfill & ={\tilde{p}}_n{Z}_{\overline{x}}\hfill \end{array}$$

(93)

where

$${\tilde{p}}_n=\tilde{p}{J}_L^{-n}$$

(94)

and

$$Z_{\overline{x}}=C_{\overline{x}}^{\left(K\right)}{R}_0$$

(95)

More explicitly, $Z_{\overline{x}}$ is the $L\times K$ matrix below:

$$Z_{\overline{x}}=\left(\begin{array}{ccc}{\overline{x}}_r& \cdots & {\overline{x}}_{r+K-1}\\ {\overline{x}}_{r-K}& \cdots & {\overline{x}}_{r-1}\\ ⋮& & ⋮\\ {\overline{x}}_{r-(L-1)K}& \cdots & {\overline{x}}_{r-(L-2)K-1}\end{array}\right)$$

(96)

were the indices are considered modulo N. The interesting feature in Equation (89) is that $Z_{\overline{x}}$ does not depend on n. Hence, we can take q different values of n and write

$$\left(\begin{array}{c}⋮\\ {\tilde{y}}_n\\ ⋮\end{array}\right)=\left(\begin{array}{c}⋮\\ {\tilde{p}}_n\\ ⋮\end{array}\right)Z_{\overline{x}}$$

(97)

assuming we know the filter frequency response (which should be the case because the filter is part of the acquisition system). More compactly, this fundamental equation can be noted:

$$Y=PZ$$

(98)

were Y is the $(q\times K)$ matrix below:

$$Y=\left(\begin{array}{c}⋮\\ {\tilde{y}}_n\\ ⋮\end{array}\right)$$

(99)

P is the $(q\times L)$ matrix below:

$$P=\left(\begin{array}{c}⋮\\ {\tilde{p}}_n\\ ⋮\end{array}\right)$$

(100)

and Z is the $(L\times K)$ matrix below:

$$Z=Z_{\overline{x}}$$

(101)

Therefore, the sizes of the matrices appearing in Equation (98) are $(q\times K),(q\times L),(L\times K).$ Then, if the number of non-zero rows in Z is less than q the matrix Z can be reconstructed from this equation using an algorithm such as OMP [29]. Eventually, from Z we can retrieve $\overline{x}$ as shown below. Indeed, it is easy to see that $\overline{x}$ can be rebuilt from Z with

$$\overline{x}=vec\left(A_{L}Z\right)J_N^{r+K}$$

(102)

where $A_L$ is the $K\times K$ anti-diagonal matrix:

$$A_L=\left(\begin{array}{cccc}0& \cdots & 0& 1\\ ⋮& & 1& 0\\ 0& & & ⋮\\ 1& 0& \cdots & 0\end{array}\right)$$

(103)

The effect of the multiplication of a matrix by $A_L$ on the left is to reverse the order of its rows.

If we have M channels instead of one in the physical system, the number of rows of Y becomes $qM$; hence, we can theoretically rebuild the signal if the number of non-zero rows in Z is less than $qM$.

Let us note $q=2\rho +\tau$ the Euclidean division of q by 2. In our experiments, we have set $r=0$ for q even and $r=-⌊K/2⌋$ for q odd. For n we take the integers in the interval $-\rho$ to $\rho +\tau -1$. This choice, while not compulsory, is designed to take into account equally distributed values around frequency 0 in $\tilde{y}$, which is a priori the best choice. Indeed, the indices of the samples taken into account go from $r-\rho K$ to $r+(\rho +\tau )K-1$, that is:

• For q odd: from $-\rho K-⌊K/2⌋$ to $\rho K+⌈K/2⌉-1$
• For q even: from $-\rho K$ to $\rho K-1$

For this choice, Figure 9 illustrates how the elements of $\overline{x}$ are arranged into matrix $Z_{\overline{x}}$ and Figure 10 illustrates how the elements of $\tilde{y}$ are arranged into matrix Y.

5.3. Application of the Equivalent Model to Reconstruction

The input of the reconstruction algorithm is the vector y provided by the acquisition device. The output is an estimation of $\overline{x}$.

We assume that:

• The frequency response $\overline{h}$ of the low-pass filter is known (or has been estimated). Then $\tilde{h}$ can be precomputed using Equation (54).
• Matrix P has been precomputed, using Equations (72), (94) and (100), or (better) has been estimated by the calibration process (see Section 6).

The procedure is as follows:

• Using an a-points FFT compute $\overline{y}$
• Compensate the filter by computing $\tilde{y}$ (Equation (78))
• Extract q sub-vectors ${\tilde{y}}_n$ from $\tilde{y}$ (Equation (83))
• Compute matrix Y (Equation (99))
• Use a compressed sensing algorithm, such as OMP [29], to estimate matrix Z from Equation (98)
• Reconstruct $\overline{x}$ from Z using Equation (102)

As an illustration of how the linear algebra model makes things simple from the programming point of view, this is the Octave program, which builds matrix Y from y and then obtains x from the reconstructed matrix Z:

ytilde = fft(y,[],2)./htilde;

ind = 1+mod(r+[-rho*K:(rho+tau)*K-1], a);

Y = reshape(ytilde(ind),K,q).’;

% Insert here estimation of Z from Y and P (using OMP, for instance)

xbs = reshape(Z(L:-1:1,:).’,1,N);

xb = circshift(xbs,[0 r+K]);

x = ifft(xb,[],2);

If there are $M>1$ physical channels, the $q\times K$ matrices Y corresponding to each channel are stacked vertically, leading to a $qM\times K$ matrix Y.

5.4. Interest of Q Even

The linear algebra model allows even values of q, instead of previous models that required q to be odd. The main interest is that it puts lower constraints on the design of the acquisition board. If the acquisition board is already available, it may also allow a better use of the MWC output data if the acquisition board is not perfectly optimized.

Let us consider our own acquisition board, which was designed before we established the linear-algebra model. We remind that the ADC sampling frequency is $F_s\simeq 104.2$ MHz and the scramblers repetition frequency is $F_p=F_{nyq}/L\simeq 10.42$ MHz. Using $F_s=F_{nyq}/b$ and $N=KL=ab$, it is easy to see that $F_p=K{F}_s/a$. In the frequency domain, the a output samples correspond to $F_s$, then $qK$ output samples correspond to $qK{F}_s/a=q{F}_p$.

• With $q=7$, we put into Y a total of $qK$ samples corresponding to a frequency half-band $q{F}_p/2=36.47$ MHz. This choice perfectly fits the frequency response of the low-pass filter, which is an almost perfectly flat and linear phase in [DC-36MHz] (see Figure 5 and Figure 6).
• With $q=6,$ we would put into Y samples corresponding to a frequency half-band $q{F}_p/2=31.26$ MHz. This corresponds to an even better area of the filter response, but in doing this, we would not use all the available information.
• On the contrary, with $q=8$, we would put it into Y samples corresponding to a frequency half-band $q{F}_p/2=41.68$ MHz. This allows us to take more information into account, but we see that we take into account some samples corresponding to a lower quality of the filter response.

This result is not surprising, because our acquisition board was designed and optimized for $q=7$, but for a future design of a new board, the possibility to have q may even be interesting because it puts fewer constraints on the choice of the commercial filters.

Indeed, applying the proposed method to our analog board did not need any change because odd values of q are allowed by our method. If a new analog board were to be designed, our method adds an additional degree of freedom because it allows even values of q also. This would allow, for example, the designer of the analog board to use a low-pass filter with a cut-off frequency of about 34 MHz (using $q=6$) instead of 40 MHz (using $q=7$).

6. Application to Fast Calibration

6.1. Proposed Approach

The objective of calibration is to estimate the true matrix P. Indeed, for real-world applications, using the theoretical matrix (Equation (100)) leads to very poor results [13]. In a previous paper [19], we proposed an approach that uses a single wideband signal for calibration, contrary to previous approaches that required successive injections of sinusoids in the system. In that paper, we presented spectrum reconstruction performances and examples of spectrum reconstruction obtained with our calibration method. In the present paper, we will mainly focus on simplifying and speeding up the method thanks to the linear algebra model.

The signal we used for calibration has a spectrum that is totally flat in the

$\left[-F_{nyq}/2,F_{nyq/2}\right]$ frequency band and random phase. Compared to methods based on iterative injections of sinusoidal waves, our new calibration method is time-saving and is more practical in terms of the simplicity of implementation because it requires a single measurement to perform the calibration. The calibration method uses advanced resynchronization preprocessing. Our calibration method offers slightly better spectrum reconstruction performances compared to reference method [13].

If we know matrices Y and Z, using Equation (98), we can estimate matrix P by:

$$\widehat{P}=Y{Z}^{+}$$

(104)

where $Z^{+}$ is the Moore–Penrose pseudo-inverse [27] of Z. Matrix Y depends on the MWC outputs, and matrix Z depends on its input signal. The problem in a real-world context is that we cannot reliably synchronize the input of the MWC with the ADC sampling, which provides the output, and even if a costly synchronization device was implemented, there are delays in the analog board itself, which are intractable. Then, the input signal must be designed such that a synchronization can be performed numerically. Otherwise, in Equation (104), we would multiply matrices Y and $Z^{+}$ corresponding to desynchronized data, which would make no sense.

In order to allow an efficient numerical synchronization, we used, for the MWC input, a periodical signal with a flat spectrum and random phase. More precisely, the period of this signal corresponds to the chosen block size, that is $N/F_{nyq}$, and one period can be represented by a length-N row vector x. This vector is generated as follows:

• A length-N vector $\overline{x}$ is generated such that, for any of its elements $\overline{x}\left(k\right)$, we have $\left|\overline{x}\left(k\right)\right|=1$ and $Arg\left(\overline{x}\left(k\right)\right)$ is random in $[0,2\pi ]$ under the constraint $Arg\left(\overline{x}(-k)\right)=-Arg\left(\overline{x}\left(k\right)\right)$ (this constraint ensures that x is real).
• x is deduced from $\overline{x}$ by an inverse FFT: $x=\overline{x}{F}_N^{-1}$

Reminding that matrix Z contains the elements of $\overline{x}$ (Equation (101)), the constant modulus $\left|\overline{x}\left(k\right)\right|=1$ ensures that no element of Z is privileged or disadvantaged by the input signal. Furthermore, the random phase ensures that the input signal has good localization properties, which is desired for efficient synchronization. Finally, choosing a periodic signal has a strong advantage: a time shift of a block taken on the input signal is equivalent to a cyclic permutation of vector x.

On the programming point of view, building Z from x is very simple:

xb = fft(x, [], 2);

xbs = circshift(xb,[0 -(r+K)]);

Z(L:-1:1,:) = reshape(xbs,K,L).’;

In the following, we will note $x_0=x$ as the signal pattern and $x_d$ a cyclic permutation, d positions to the right, of the pattern. This means that

$$x_d=x{J}_N^d$$

(105)

The procedure that we propose is as follows:

• Feed the acquisition device with a periodic signal, which is a repetition of a known pattern $x_0$
• Record a samples at the output of the acquisition device (this is vector y), then compute matrix Y using steps 1 to 4 of the reconstruction procedure.
• Perform cyclic permutations of the input pattern $x_0$. For each vector $x_d$, compute matrix $Z_d=Z_{{\overline{x}}_d}$ (using Equations (96) and (101)) and then $\widehat{P}$ (using Equation (104)).

Let us note the residue

$$R_d=Y-{\widehat{P}}_d{Z}_d.$$

(106)

The criterion to determine the best cyclic permutation of the input pattern $x_0$ is the inverse Frobenius norm of $R_d$ (the Frobenius norm is the square root of the sum of the square modules of the elements of the matrix).

The computational cost per iteration (i.e., per value of d tested) can be estimated as follows:

• A FFT is required to compute ${\overline{x}}_d$, that is $LK{log}_2\left(LK\right)$ flops (because $N=LK$).
• Computation of the pseudo-inverse $Z_d^{+}$: $7K{L}^2+20L^3$ flops (see Section 2.7).
• Computation of $\widehat{P}=Y{Z}_d^{+}$ (the sizes of the matrices are $qM\times K$ and $K\times L$): $qMKL$ flops.
• Computation of $\widehat{P}{Z}_d$ requires $qMLK$ flops.
• Computation of the Frobenius norm requires $KL$ flops

Globally, since the computation of the Frobenius norm can be neglected compared to the other terms, the algorithm requires about $L(2qMK+7LK+20L^2+K{log}_2\left(LK\right))$ flops per iteration.

6.2. Fast Update of Matrix Z

While evaluating all possible shifts d, computation of matrix $Z_d=Z_{{\overline{x}}_d}$ requires an N-points FFT to obtain ${\overline{x}}_d$, which requires approximately $Nlo{g}_2\left(N\right)$ multiplication at each iteration. However, we can reduce the complexity just by computing the first matrix and then updating it at each iteration as described below. Let us consider a vector $x_d$, which is a cyclic permutation, d positions to the right of pattern x. We have:

$$x_d=x_0{J}_N^d$$

(107)

Then

$$\begin{array}{cc}\hfill {\overline{x}}_d& =x_0{J}_N^d{F}_N\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill & =x_0{F}_N{F}_N^{-1}{J}_N^d{F}_N\hfill \end{array}$$

(109)

$$\begin{array}{cc}\hfill & ={\overline{x}}_0{D}_{{\overline{\alpha }}_d}\hfill \end{array}$$

(110)

where

$${\alpha }_d=\left[\underset{d}{\underbrace{0\cdots 0}}10\cdots 0\right]$$

(111)

Indeed, since $J_N^d=C_{{\alpha }_d}$, using Equation (15) we have:

$$\begin{array}{cc}\hfill F_N^{-1}{J}_N^d{F}_N& =F_N^{-1}{C}_{{\alpha }_d}{F}_N\hfill \end{array}$$

(112)

$$\begin{array}{cc}\hfill & =D_{{\overline{\alpha }}_d}\hfill \end{array}$$

(113)

Since ${\overline{\alpha }}_d={\alpha }_d{F}_N$, it is easy to see that ${\overline{\alpha }}_d$ is the ${(d+1)}^{th}$ row of $F_N$, that is (see Equations (10) and (11)):

$${\overline{\alpha }}_d=\left[\begin{array}{ccccc}1& {\omega }^d& {\omega }^{2d}& \cdots & {\omega }^{(N-1)d}\end{array}\right]$$

(114)

Then, we have:

$$\begin{array}{cc}\hfill Z_d& =Z_{{\overline{\alpha }}_d•{\overline{x}}_0}\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill & =Z_{{\overline{\alpha }}_d}•Z_0\hfill \end{array}$$

(116)

where • stands for element by element multiplication. This equation can also be written

$$Z_d={\omega }^{-rd}{\mathsf{\Theta }}_d^{*}{Z}_0{\mathsf{\Omega }}_d$$

(117)

with

$${\mathsf{\Theta }}_d=diag\left[\begin{array}{ccccc}1& {\omega }^{Kd}& {\omega }^{2Kd}& \cdots & {\omega }^{(L-1)Kd}\end{array}\right]$$

(118)

$${\mathsf{\Omega }}_d=diag\left[\begin{array}{ccccc}1& {\omega }^d& {\omega }^{2d}& \cdots & {\omega }^{(K-1)d}\end{array}\right]$$

(119)

To see where this formula comes from, we must remind that multiplication by a diagonal matrix on the left (right) multiplies the rows (columns) by the elements of the diagonal. Then, denoting ${\omega }_d={\omega }^d$ we can see that:

$$\begin{array}{cc}\hfill {\omega }_d^{-r}vec{\left({\mathsf{\Theta }}_d\right)}^{*}vec\left({\mathsf{\Omega }}_d\right)& ={\omega }_d^{-r}\left[\begin{array}{c}1\\ {\omega }_d^{-K}\\ \cdots \\ {\omega }_d^{-(L-1)K}\end{array}\right]\left[\begin{array}{cccc}1& {\omega }_d& \cdots & {\omega }_d^K\end{array}\right]\hfill \end{array}$$

(120)

$$\begin{array}{cc}\hfill & ={\omega }_d^{-r}\left[\begin{array}{cccc}1& {\omega }_d& \cdots & {\omega }_d^{K-1}\\ {\omega }_d^{-K}& {\omega }_d^{-K+1}& \cdots & {\omega }_d^{-1}\\ ⋮& ⋮& & ⋮\\ {\omega }_d^{-(L-1)K}& {\omega }_d^{-(L-1)K+1}& \cdots & {\omega }_d^{-(L-2)K-1}\end{array}\right]\hfill \end{array}$$

(121)

$$\begin{array}{cc}\hfill & =Z_{{\overline{\alpha }}_d}\hfill \end{array}$$

(122)

If we evaluate by step g, we can use:

$$Z_d=Z_{{\overline{\alpha }}_g}•Z_{d-g}$$

(123)

Matrix $Z_{{\overline{\alpha }}_g}$ can be precomputed. Therefore, at each iteration, we need only N multiplications, which is less complex than computing ${\overline{x}}_d$ each time. This update requires only $N=LK$ multiplications at each iteration, instead of approximately $Nlo{g}_2\left(N\right)$.

If we want to allow sub-sample precision (i.e., $g<1$), we just have to note $\eta =N/2$ and write ${\overline{\alpha }}_g$ as follows:

$${\overline{\alpha }}_g=\left[\begin{array}{ccccccc}1& {\omega }^g& \cdots & {\omega }^{\eta -1}& {\omega }^{-\eta g}& \cdots & {\omega }^{-g}\end{array}\right]$$

(124)

This is very interesting because sub-sample precision is then allowed without any additional cost due to oversampling (with this method, no oversampling is required). Matrix $Z_{{\overline{\alpha }}_g}$ is computed as follows:

alpha = exp(-i*2*pi*g*[0:eta-1 -eta:-1]/N);

alpha = circshift(alpha,[0 -(r+K)]);

Zalpha(L:-1:1,:) = reshape(alpha,K,L).’;

Then, at each iteration, updating matrix Z is achieved by:

Z = Z .* Zalpha;

6.3. Fast Update of Matrix $Z^{+}$

The approach is similar to the pseudo-inverse: we can reduce the complexity just by computing the first pseudo-inverse matrix and then updating it at each iteration, as described below. This update requires only N multiplications at each iteration.

According to discussions above, denoting $Z_0=US{V}^{*}$ as the SVD of $Z_0$, we have

$$\begin{array}{cc}\hfill Z_d& ={\omega }^{-rd}{\mathsf{\Theta }}_d^{*}\left(US{V}^{*}\right){\mathsf{\Omega }}_d\hfill \end{array}$$

(125)

$$\begin{array}{cc}\hfill & =\left({\omega }^{-rd}{\mathsf{\Theta }}_d^{*}U\right)S{\left({\mathsf{\Omega }}_d^{*}V\right)}^{*}\hfill \end{array}$$

(126)

This equation directly provides the SVD of $Z_d$. It follows that $Z_d^{+}$ is

$$\begin{array}{cc}\hfill Z_d^{+}& =\left({\mathsf{\Omega }}_d^{*}V\right)S^{+}{\left({\omega }^{-rd}{\mathsf{\Theta }}_d^{*}U\right)}^{*}\hfill \end{array}$$

(127)

$$\begin{array}{cc}\hfill & ={\omega }^{rd}{\mathsf{\Omega }}_d^{*}{Z}_0^{+}{\mathsf{\Theta }}_d\hfill \end{array}$$

(128)

This update can be realized by element-by-element multiplication:

$$Z_d^{+}=Z_{{\overline{\alpha }}_d}^{*}•Z_0^{+}$$

(129)

If we evaluate by step g we can update the matrix at each iteration by:

$$Z_d^{+}=Z_{{\overline{\alpha }}_g}^{*}•Z_{d-g}^{+}$$

(130)

where $Z_{{\overline{\alpha }}_g}^{*}$ can be precomputed.

6.4. Fast Synchronization

To obtain a synchronization, we must evaluate the Frobenius norm of $R_d$ at each iteration.

We can evaluate the computational complexity of this fast algorithm as follows:

• Updating matrix $Z_d$ and $Z_d^{+}$ requires $2LK$ flops.
• Computing $Y_d{Z}_0^{+}$ requires $qMKL$ flops.
• Multiplying $Y_d{Z}_0^{+}$ by $Z_0$ requires $qMKL$ flops.
• Computing the Frobenius norm requires $KL$ flops.

Globally, the algorithm requires $2qMKL+3KL\simeq 2qMKL$ flops per iteration, which is to be compared to $L(2qMK+7LK+20L^2+K{log}_2\left(LK\right))$ for the previous version. The gain is, approximately:

$$\begin{array}{cc}\hfill G& \simeq 1+\frac{7L}{2qM}+\frac{10L^2}{qMK}+\frac{{log}_2\left(L\right)}{2qM}+\frac{{log}_2\left(K\right)}{2qM}\hfill \end{array}$$

(131)

$$\begin{array}{cc}\hfill & \simeq 21\hfill \end{array}$$

(132)

Computation time on Octave is 54 s for the slow version and 1.6 s for the fast version, which is then 34 times faster.

As a function of K, the gain decreases until $K=20ln\left(2\right)L^2=127761$ (which is a huge value, not expected for real-world acquisition devices), where it reaches a minimum of 13.4, then increases slowly, behaving asymptotically as ${log}_2\left(K\right)/\left(2qM\right)$, as shown in Figure 11.

7. Experimental Results

Using a step $g>1$ for the evaluation of the synchronization criterion allows a decrease in the computation time by a factor g. A good strategy is to obtain a coarse synchronization with a step $g>1$, and then to perform a fine synchronization around the detected peak. The fine synchronization may even be realized at sub-sample precision, if desired. However, a too large initial step must be avoided because it may lead to missing the synchronization peak during the coarse synchronization. In our experiments, we first used a coarse synchronization with step $g=16$, then a fine synchronization with step $g=1$ around the coarse synchronization peak. Figure 12 shows the obtained synchronization data. Figure 13 and Figure 14 present are magnifications of the synchronization peak to show more details.

If the response of the filter is not taken into account in Equation (78) (i.e., assuming an ideal low-pass filter), the synchronization peak is only slightly lower (4.75 instead of 5.22). This is due to the fact that the low-pass filter used in our analog board has good characteristics (almost flat response, and almost linear phase, in the band of interest). The difference would be higher with a lower quality filter. Anyway, it is always better to integrate the filter response in the equations, as we did, because the additional computational cost is negligible.

Once the signal is synchronized, the estimated matrix $\widehat{P}$ is obtained using Equation (104) (at no additional cost because this computation was already part of the synchronization process). Figure 15 shows the modulus of the matrix elements. Here, since we have $M=4$ physical channels and we have taken $q=7$, the matrix has $qM$ = 28 rows (and $L=96$ columns). The first $q=7$ rows correspond to the first physical channel, then the next q rows to the second physical channel, and so on.

The theoretical (ideal) matrix $P_{th}$ can be computed using Equation (100). Since the analog scrambling sequence is the output of a Digital to Analog Converter (DAC), fed with a digital pseudo-random sequence, it is (ideally) piecewise constant in the time domain. In the spectral domain, this is equivalent to a multiplication by the $sinc$ function in which first zero is at $F_{nyq}$. We take that into account when computing our theoretical matrix in order to be as close as possible to the real-world matrix. Figure 16 shows the modules of the elements of this matrix.

The estimated (real-world) mixture matrix $\widehat{P}$ may be compared with the theoretical (ideal) mixture matrix $P_{th}$. On the basis of the elements modules, we can see that their overall aspects are close despite noticeable differences. In fact, the main differences are in the phases of the elements. If we draw the normalized correlation coefficients between the columns of both matrices (Figure 17), we obtain low values, which confirms significant differences. We remind that a normalized correlation coefficient is the absolute value of the cosine of the angle between two vectors (here the columns of both matrices), then values around 0.5 mean that the angle is about 60 degrees, and thus, the columns are significantly different.

In a previous paper [19], we showed that despite the good quality of our real-world acquisition board, calibration of the system is absolutely required: using the theoretical matrix leads to poor reconstruction performances. Without calibration, the system usually incorrectly detects the active sub-bands, and even when the active sub-bands are correctly identified, the spectrum reconstruction provided by the uncalibrated system is extremely poor, as illustrated in Figure 18.

Figure 19 shows the relative error between computed and observed system output, as a function of output frequency. The computation is performed using a formula similar to the criterion used in [32].

$$\epsilon \left(f\right)=20{log}_{10}\left(\frac{∥outpu{t}_{real}\left(f\right)-outpu{t}_{computed}\left(f\right)∥}{∥outpu{t}_{real}\left(f\right)∥}\right)$$

(133)

Here the frequency band of the acquisition system output signal has been divided into 28 subbands, and $outpu{t}_{real}$ is the real output signal comprised in the subband centered on f. Similarly, $outpu{t}_{computed}\left(f\right)$ is the computed output signal comprised in the subband centered on f, “computed” meaning that it is computed from the input signal using the system model and a given matrix P. The three curves correspond to using different matrices P in this computation. It can be seen that our method and the reference method provide a good prediction of the real output (relative errors around −18 dB), while the theoretical uncalibrated matrix provides poor results (large relative errors). We remind that in any case, even for the theoretical matrix, we always use a low-pass filter that is calibrated: the true frequency response of the filter (Figure 5 and Figure 6) is being taken into account in the equations through the diagonal matrix $D_{\tilde{h}}$.

To further illustrate the interest of the approach, Figure 20 shows the spectrum of the reconstructed signal using different matrices P. On top, it is the true spectrum. Here we have two transmitters in the monitored bandwidth. We can see that our method, as well as the reference method [13], correctly identifies the presence of two active transmitters. On the contrary, when using the theoretical matrix, the number of transmitters and their frequency locations are incorrect.

If we zoom in on Figure 20, we can see (Figure 21) that the frequency location of the second transmitter is more precise when using matrix P provided by our method than when using the matrix provided by the reference method.

Let us now consider a much more difficult case, where six transmitters are simultaneously active in the monitored bandwidth (Figure 22). Using the matrix P provided by our approach leads to a correct estimation of the number of transmitters and of their frequency location, despite imperfect shapes of the spectrum reconstruction for some transmitters. The reference method produced quite good results also, but missed one transmitter (around 330 MHz) and detected a non-existing transmitter (around 180 MHz). Finally, as anticipated, using the theoretical matrix P lead to a reconstruction that is almost totally false.

This shows that calibration is unavoidable. An interest of the extremely fast calibration procedure proposed in this paper is the possibility of performing quick recalibration of the system as soon as the performances appear to decrease. Indeed, many factors, such as temperature, external perturbation, components aging, etc., modify the characteristics of the system, making a recalibration necessary.

8. Conclusions

In this paper, we have established an MWC model solely based on linear algebra. It is very convenient as a basis for fast and efficient programming of simulation, reconstruction and calibration algorithms related to MWC. It suppresses a previous restriction on the channels augmentation factor, hence providing more degrees of liberty to the systems designer. It also allowed us to develop an extremely fast implementation of a previously proposed calibration algorithm, leading to a gain of a factor greater than 20 on the computation time. This fast calibration allows quick recalibration of the system as soon as it becomes necessary. Our future work will include more in-depth exploitation of the advantages and interesting properties of this model.