Matched Filtering with Polyphase Filters

These notes use the thesis Design and Implementation of a Broadband 16QAM/16APSK Modem with a Flexibly Variable Rate as study material.

Implementing variable-rate pulse shaping

The pulse-shaping expression

When the overall response of a baseband system is equivalent to an ideal low-pass filter, intersymbol interference can be eliminated. An ideal low-pass filter cannot be implemented in practice, however. A root-raised-cosine filter is therefore used at the transmitter for pulse shaping, and another root-raised-cosine filter is used at the receiver for matched filtering. This makes the overall channel equivalent to a filter with a response close to the desired low-pass characteristic: the response is greatest at the sampling instant for the current symbol ($t=0$) and zero at the sampling instants of the other symbols ($t=\pm k T_{s}$).

The time-domain response of a root-raised-cosine filter is

Deriving rate conversion through digital resampling

What does variable rate mean, and why is rate conversion necessary? Why do we need high-speed parallel modulation with a variable symbol rate?

1. The sampling rate of the transmitter’s DAC or receiver’s ADC is usually fixed at the start of the design so that the implementation operates correctly. The symbol rate, however, may be adjusted for the application. At a fixed sampling rate, this changes the oversampling factor, which can sometimes be fractional. The pulse-shaping filter therefore needs a dedicated design with variable-rate support to accommodate input signals at different rates.
2. The processing clock inside an FPGA generally does not exceed 200 MHz. Achieving modulation at gigabits per second requires a fully parallel architecture, so the modulator must be parallelized.

According to the Nyquist sampling theorem, provided the sampling frequency exceeds twice the signal bandwidth, the original signal $y(t)$ can be reconstructed from its samples $x(nT_s)$:

Here, $x(t)$ is the input signal, $h(t)$ is the filter’s time-domain response, and $T_s$ is the symbol period.

The physical meaning is that obtaining the output signal $y(t)$ requires convolving the input signal with the filter. Notice that the input is discrete, with sampling instants $nT_s$, while the filter’s time-domain response is centered at $t$. This is a familiar expression, but understanding it thoroughly is essential to following the derivation below.

If we resample this reconstructed signal—that is, obtain its value at a desired point—the first idea is to replace $t$ in the expression above with the time of that point, such as $kT_{sp}$. This gives

Here, $T_{sp}$ can be understood as the DAC sampling period. Physically, we want to modulate the symbols at the specified sampling period $T_{sp}$, and the expression above gives $y(kT_{sp})$. The information-symbol period is currently $T_s$, corresponding to the lower rate, while the DAC samples with period $kT_{sp}$, corresponding to the higher rate. We therefore use the expression to pulse-shape the slower information symbols into a faster signal for the DAC to output.

We can see that $x(nT_s)$ is independent of $x(T_{sp})$ and depends only on $h(t)$.

To simplify the expression, let $m_{k}=\left\lfloor\frac{k T_{sp}}{T_{s}}\right\rfloor$ and $\mu_{k}=\frac{k T_{sp}}{T_{s}}-m_{k}$. Substituting gives

Then let $m=m_{k}-n$, giving

After this simplification, $y(kT_{sp})$ depends only on $x\left(\left(m_{k}-m\right) T_{s}\right)$ and $h\left(\left(\mu_{k}+m\right) T_{s}\right)$. The integer parameter involved is $m_k, m$, and the fractional parameter is $\mu_k$. For $x(t)$, the index must be an integer, so $m_k-m$ satisfies the requirement and selects an ordinary element from the $x$ sequence. For $h$, the fractional part means that we first oversample $h$ and then perform ordinary selection. The position represented by $m$ is the correct sampling point before the offset is applied; adding $\mu_k$ introduces the offset.

A question

To me, the two expressions above look equivalent. But once their ranges are specified, are they still equivalent?

An infinite summation cannot be implemented, so the filter can be truncated:

For both $x$ and $h$, $T_s$ indicates selection within a sequence and acts as a scale factor. We can therefore set that scaling relationship aside and focus on mapping $x((m_k-m)-T_s)$ to $x(m_k-m)$.

Designing the variable-rate pulse-shaping filter

On an FPGA, $h_{T}\left(\mu_{k}+m\right)$ cannot be calculated in real time, so a lookup table is used. The oversampling factor determines the precision of $\mu_k$. If the precision of $\mu_k$ is set to $\frac{1}{2^{-12}}$, then at a sampling rate of 5 Gsps the variable-rate step size is $1.22Msps$.

Following the previous section, suppose that the effects of symbols more than three symbol intervals before or after the current symbol are negligible. Truncating the filter gives

Implementation:

Use MATLAB’s rcosdesign function to generate coefficients with α = 0.35, an oversampling factor of 4096, and a span of six symbols across both sides. Quantize the coefficients to 12 bits, then split them by symbol into six groups of 4096 coefficients. Store the groups in six ROM cores, each 12 bits wide and 4096 entries deep. Scale both $\mu_k$ and $m$ by $2^{12}$ to obtain address indices, and retrieve the six coefficients required for the calculation.

At a sampling rate of 5 Gsps and a clock frequency of 200 MHz or less, 32 parallel paths are required. The expression can be rewritten as

Here, $n$ is the time index and $i$ is the parallel-path index. $m_{n}(i)=\left\lfloor\frac{(32 n+i) T_{s p}}{T_{s}}\right\rfloor$ and $\mu_{n}(i)=\frac{(32 n+i) T_{s p}}{T_{s}}-m_{n}(i)$.

The implementation block diagram is shown below. The input symbol data first passes through a FIFO and then enters DataCache. COEGen calculates $m_n(i)$ and $\mu_n(i)$, which select six pairs of symbols and filter coefficients from DataCache and FilterRom. Multiplying and summing those pairs produces one pulse-shaped output path.

Some parameter choices for the variable-rate pulse-shaping module

The highest sampling rate used in this project is $5Gsps$, and the highest symbol rate is $1200Msps$. Under these conditions, $\frac{32T_{sp}}{T_s}=7.68$. Thus, with 32 parallel paths, $m_n(1)$ and $m_n(32)$ differ by at most 8, so $m_n(i)$ ranges over $[0, 7]$. Accounting for $m_{n}(1)-2$ and $m_{n}(32)+3$, a total of 13 symbols is needed for 32-path parallel pulse shaping. DataCache therefore buffers one clock cycle of the eight parallel symbols arriving from the FIFO. Together with the data from the next cycle, this forms 16 parallel symbols and meets the filter’s rate requirement.

For register addressing, we represent $\frac{(32 n+i) T_{s p}}{T_{s}}$ using 53 bits to achieve the required precision: three integer bits and fifty fractional bits. In other words, $m_n(i)$ uses the three integer bits, while $\mu_n(i)$ uses the highest twelve fractional bits.

For the FIFO readout strategy, refer directly to the source material.

Implementing a high-speed parallel matched filter

The $Z$-transform

Applying $\delta(n)$ to a system gives its impulse response $h(n)$; for example, $h(n)=\left[\begin{array}{lll}1 & 2 & 1\end{array}\right]$.

Now suppose another signal is applied to the system, for example $x(n)=[1\quad 2\quad 1]=\delta(n)+2 \delta(n-1)+\delta(n-2)$.

The output should be $y(n)=[1\quad 4\quad 6\quad 4\quad 1]=h(n)+2 h(n-1)+h(n-2)=x(n)* h(n)$, obtained by convolution.

On inspection, this is simply polynomial multiplication of $[1\quad 2\quad 1]$ and $[1\quad 2\quad 1]$!

We therefore rewrite the signal $x(n)$ as $X(Z)=1+2 Z^{-1}+Z^{-2}$ and likewise write the impulse response as $H(Z)=1+2 Z^{-1}+Z^{-2}$. Then $Y(Z) = 1+4 Z^{-1}+6 Z^{-2}+4 Z^{-3}+Z^{-4} = X(Z)H(Z)$. This is called the $Z$-transform.

⭐️ $Z^{-1}$ corresponds to data from the next clock cycle—the future.

Derivation

An Nth-order FIR filter can be expressed in the time domain as

Its Z-transform is

Following the thesis, the $z$-transforms of sequences $x(n)$ and $h(n)$ can be separated into the $Z$-transforms of their even- and odd-indexed subsequences:

Similarly,

The corresponding $y(n)$ sequence can also be written as the sum of even- and odd-indexed subsequences:

In matrix form,

For L parallel paths, this extends to

Implementing the high-speed parallel matched filter

Here, there are 32 parallel paths, while the I and Q inputs each have only 16 paths. The filter coefficient length is 25, corresponding to rcosdesign(0.35, 4, 6, 'sqrt').

For simplicity, first assume a coefficient length of 32, derive the implementation, and then simplify it to the length-25 case. For the I branch, the input $x(n)$ contains only the even-indexed sequence: it is nonzero only for $x(2k)$ and identically zero for $x(2k+1)$. Nevertheless, it still corresponds to two output groups of length 16:

Equivalently,

The same reasoning applies to the Q branch:

Expanding the elementwise products shows that the first column of the $H$ matrix is multiplied by $X_0$, the third column of the $H$ matrix by $X_2$, and the 31st column of the $H$ matrix by $X_{30}$. The products to be used subsequently are

Because the filter coefficients are symmetric, $H_0X_0$ and $H_{31}X_0$ are actually the same, so only half of the filter coefficients are needed. In addition, half of the matrix operations involve a delay of 32 clocks; with 32-way parallel processing, this amounts to a delay of only one processing clock cycle. The results can therefore be split into two parts. One part is delayed before entering the pipelined adder, while the other enters it directly. Pipelined addition completes the low-pass filtering operation. The block diagram is shown below:

For the matched-filter implementation, the input is likewise 32-way parallel, but its order is only 25. Thus, $H_{25}-H_{31}$ are all zero, while the remaining coefficients $H_0-H_{24}$ are nonzero. The algorithm can therefore be written as

Its implementation is almost identical to that of the low-pass filter, except that operations associated with the zero coefficients can be omitted.

For convenient reference when implementing this in MATLAB, the matrix expressions and summation formulas for the 32-way parallel I and Q branches are listed below:

The 64-way parallel case is also listed here: For the I branch,

The corresponding matrix expression is

Similarly, for the Q branch,