Polyphase Decimation Filter parallel inputs

I have an FPGA that's streaming in samples from a high speed ADC. The FPGA clock rate cannot run at the sample rate, so the samples are streamed in parallel (4 samples per clock cycle).

I'd like to design a polyphase decimation filter to split the input spectrum into 4 downsampled subbands, with an effective rate of fs/4 per "channel". In that way, for each subband channel, there's only 1 sample per FPGA clock cycle.

this post seemed to be doing something similar, but when I tried the accepted solution, it's almost working, but I am seeing something strange:

So each color represents one out of the 4 output channels. I input a LFM chirp across the entire input spectrum and plotted which output channel has energy. As you can see, each channel is picking two 1/8th subbands instead of a single 1/4th subband.

Any ideas how to fix this?

code for LFM generation:

  f1 = 1e6;         %Start Freq
f2 = 1250e6;      %Stop Freq
duration = 200e-6;  %Pulse Duration
fc = (f2+f1)/2;
fs_dac = fs;
B  = f2 - f1;
tau = duration
sampCount = tau*fs_dac;
kVect = [0:sampCount];
A = 2*pi*(fc-B/2);
C = pi*B/tau;
Ht = kVect;
J = (2*pi*f1)/fs;
L = C/fs^2;
totalWaveform = exp(1i*(Ht.*(J+L.*Ht)));


Code for filter design:

  h0 = firls(numTaps, [0 0.25 0.25+0.01 1], [1 1 0 0]);
h1 = h0 .* exp(-1j*2*pi*1.*[0:numTaps+1]'/(4));
h2 = h0 .* exp(-1j*2*pi*2.*[0:numTaps+1]'/(4));
h3 = h0 .* exp(-1j*2*pi*3.*[0:numTaps+1]'/(4));


Here's a diagram of what I desire:

and here is how I'm attempting to do the polyphase filtering:

x0 = x(1:4:end);
x1 = x(2:4:end);
x2 = x(3:4:end);
x3 = x(4:4:end);

y0 = filter(h0(1:4:end),1,x0);
y1 = filter(h1(2:4:end),1,x1);
y2 = filter(h2(3:4:end),1,x2);
y3 = filter(h3(4:4:end),1,x3);

• What's the sample rate of the DAC and the FPGA ? Jul 22, 2022 at 0:01

Problem Statement and Bottom Line of Solution

As described, the OP wishes to channelize the input bandwidth that extends from DC to fs/2 (f = 0 to +1.25 GHz) into four channels each having one quarter of the bandwidth, and each of these is to be frequency translated to the new first Nyquist zone covering DC to fs/8 (f = 0 to 312.5 MHz), consistent with the sampling rate for each of the four data paths of 625 MHz, as depicted in the graphic below.

It was assumed that real inputs and outputs are needed. An optimized solution to achieve this result is given in the diagram below with implementation further detailed as "Approach 2" in the section on "Polyphase Channelization Approach". This solution can provide either real outputs as shown, or complex outputs where the resulting waveforms will be the analytic signals (all negative frequency axis values will be zero). When real outputs are desired, the Band1 and Band4 processing paths are completely real (real input, real filter bank), real output.

The filter is a bank of four polyphase filters all running at rate $$R/4$$. The exponent $$m$$ is the sample index at the R/4 rate, and the signal $$-(-1)^m$$ serves to correct a spectrum reversal on Band2 and Band4 while maintaining proper synchronization with the other channels. Given the rotations of +/-j, rotating and then selecting the real part results in using only the real filter or imaginary filter in each path. See the code and further details of this implementation toward the end of this post in the section "Demonstration of Approach 2".

Recreating OP's Result

I was able to recreate the OP's result with the following assumptions of some missing details:

• fs = 2.5 GHz (sampling rate)

• The FPGA sampling rate is 625 MHz (fs/4).

• The input signal is real and the desired 1/4 rate output signals are to be real.

The OP tested the channelization by sweeping a chirp signal which I confirmed was constructed correctly to sweep from DC to 1.25 GHz, properly accounting for the effect of the ramp itself on generating the instantaneous frequency that starts at DC and ends at 1.25 GHz. The OP did not detail the creation of the filter beyond the code provided, so I suspect that the filters were incorrectly implemented to each run directly at the output rate of 625 MHz, instead of decimating the filter coefficients from higher rate filters designed for the the 2.5 GHz rate (at least with this assumption I get the same result). Below is the frequency response in the first Nyquist zone for the four filters created, matching the result for the swept chirp signal (as these responses would repeat periodically in each higher Nyquist zone, which is what we see in the OP's result with these same frequency translations):

The above is the direct result showing the frequency response for the four filters created using the code given in the OP if we assume those filters are running directly at the decimated sampling rate of 625 MHz.

Polyphase Channelization Approach

To accomplish what the OP intends, we could instead proceed as follows:

Approach 1: Use Four Quarter-Band Filters

• Design quarter-band filters for each of the four bands and implement each channel as a polyphase decimator as described in this post for each of the four output bands. The polyphase decimators would each run at rate $$f_s/4$$ and appropriately frequency translate the respective bands to the new output rate for each channel running at $$f_s/4$$. This implementation is depicted in the graphic below.

For each band the grid shown represents the bank of 4 real polyphase FIR filters where the coefficients of each quarter-band filter are loaded in row to column form. The commutators all start at the bottom filter in each bank, where the first input sample is loaded into all bottom filters in parallel. The commutator then moves up (at the input rate) and the next sample is loaded. Once all filters are loaded (after 4 input samples) the filter outputs are computed, the data is shifted in and then the cycle repeats. Thus the filters all run at one quarter of the input rate.

The result of each polyphase decimator is mathematically equivalent to filtering the input waveform with its respective quarter-band filter at the input rate, and then selecting every fourth sample at the output (down-sampling), with the advantage of having all operations running at the lower rate. Note that this process will create spectral inversion whenever the waveform prior to decimation occupies an even Nyquist Zone, as depicted in the graphic below. For real signals, this is corrected by changing the sign for every other sample, as done at the Band2 and Band4 outputs in the implementation diagram above. (m represents the sample index at the output rate starting with index 0).

Approach 2: Freq Translate the Signal through One Filter

• The bandpass and highpass filter can be constructed by rotating the coefficients (multiply by a complex rotating phasor in time domain, which results in a frequency shift in the frequency domain), which frequency translates the filter to the passband of interest, as done in the OP. Alternatively, we can frequency translate the input signal to the filter by rotating the input for each channel. The advantage of this is that we can ultimately move the complex rotation to the filter outputs prior to each summation of the polyphase filters and thus reuse the same filter outputs across all the bands.

The first step to see this approach in comparison to Approach 1 is the intermediate approach depicted in the graphic below. Band 1 proceeds the same way as before, but for Band 2, 3 and 4 we first frequency translate the input signal so that each band of interest passes through the Band 1 filter (we move the signal to the filter instead of moving the filter to the signal). Since we desire a solution for a real input and ultimately a real output, we must use complex single-sideband filters to filter out the negative sideband after frequency translation for Band2 and Band3 (otherwise when we finally take the real part of the resulting output we would get a interference alias back into our final output.). Since the negative sideband is adjacent for Band1 and Band4, we can process those bands with a completely real filter-bank and data paths.

The details on how the rotations were derived is depicted in the graphic below.

The particular translations were chosen for the resulting product term of +1, -1, +j and -j due to the simplicity of implementing those operations in hardware and the resulting reduction in the filter processing required. For example, the resulting rotations by $$j$$ and $$-j$$ involve swapping the complex $$I$$ and $$Q$$ outputs for that path and multiplying by $$(-1)^n$$:

rotation by $$-j$$: $$Q[n]=-I[n](-1)^n$$ $$I[n]=Q[n](-1)^n$$

rotation by $$+j$$: $$Q[n]=I[n](-1)^n$$ $$I[n]=Q[n](-1)^n$$

Further, given a real input and real output; for paths rotated by $$\pm j$$, the resulting implementation will alternate between the real and imaginary coefficients for the four filters in the polyphase filter bank, using only one of those filters in each of the four polyphase paths.

I show fatter traces representing the complex datapaths, along with an operation to select the real components at the filter output. This shows the generic approach if the rotations weren't in multiples of $$j$$. The Band2 and Band3 filters consist of the real Band1 filter for the real components of the datapath, and its Hilbert Transform for the quadrature components of the datapath; resulting in 8 total filters in the filter bank. However, since only the real outputs will be selected at the output, half of the filters will be used (two of the filters are the real filters and the other two are the quadrature filters). Which two depends on the phasing of the input: the diagram I drew is associated with $$n$$ in the diagram starting at index 1, resulting in the second and fourth filter in the bank to use the quadrature filter and the first and third filter to use the real filter.

The significance of using complex filtering for this functional block diagram is depicted below showing the result of filtering after the signal translations. Shown is the selection of Band2 after multiplying the input signal by $$j^n$$ with $$n$$ as the index at the input sample rate, starting with $$n=1$$, which translates the spectrum to the right by the width of two bands, which is then passed through the Band1 filter. If the Band1 filter was real, the negative frequency image of that filter would select the translated Band3 as well, resulting in aliasing. By using a complex filter, we only select Band2, and once filtered, we can then take the real part of that complex result to get a real output representing the real input spectrum for band2.

At this point the approach shown is merely educational as a step in understanding the ultimate "Approach 2" solution. As detailed by fred harris (see reference at bottom of this post), we can achieve significant simplifications by moving the phase rotators from the filter inputs to the outputs of the individual filter banks in the polyphase implementation. Thus the filter inputs will now all be real, and the same to each filter, which means we only need to use one filter bank and rotate the outputs with a complex multiplier before each summer. With different rotations we can achieve the different outputs for each Band, which results in the final solution, repeated below:

This implementation consists of only one filter bank containing four complex polyphase filters with a real input; thus each polyphase filter consists of one "I" and one "Q" filter, resulting in 8 filters total to support the complex filter bank. The fatter traces in the diagram represent complex (I and Q) data-paths. As mentioned earlier, only the real filters are needed for the processing of Band1 and Band4, and alternating real and quadrature filters are used for processing Band2 and Band3. Since the imaginary filter for the top and bottom filter in the polyphase structure is never used, 6 individual filter elements each running at R/4 are actually needed.

For the processing of Band2 and Band3; once the complex signal is properly rotated, the summation is only done on the real component which selects the real output (which simplifies to only using the real or quadrature filter for the successive outputs), and doubled to match the scaling of Band1 and Band4. The resulting Band2 and Band4 output is multiplied by $$-(-1)^n$$ to correct for spectral reversal in the processing.

The complex coefficients for the polyphase filters are required in order to process real inputs and provide ultimately real outputs without any aliasing.

Demonstration of Approaches

To test and demonstrate the recommended approaches, I created a waveform that contains four sinusoids of the same magnitude and frequencies 200 MHz, 410 MHz, 800 MHz and 1.2 GHz. The graphic below demonstrates the expected result:

The code for the test waveform is as follows (each is scaled differently to identify in both frequency and amplitude in final result):

# Test with a real tone in each sub-band
nsamps = 2**16
time = np.arange(nsamps)/fs
tone1 = 0.01 * np.cos(2*np.pi*200e6*time)
tone2 = 0.05 * np.cos(2*np.pi*410e6*time)
tone3 = 0.1 * np.cos(2*np.pi*800e6*time)
tone4 = 0.5 * np.cos(2*np.pi*1200e6*time)
signal = tone1 + tone2 + tone3 + tone4


Demonstration of Approach 1

The quarter-band filters could be directly created with the least squares algorithm for each targeted band, but I used the approach similar to the OP of creating a base low-pass filter and rotating the coefficients to frequency translate the filters to each band. I took the real of the resulting coefficients to demonstrate the implementation with a real input signal and real filters:

# create base filter
numTaps = 199
hb = sig.firls(numTaps, [0, .1, .15, 1], [1, 1, 0, 0])
hb = np.append(hb,0) # make length divisible by 4 for polyphase decomposition

# frequency translate (rotate coefficients) to each band center
# and use real component of resulting coefficients
wn = np.exp(1j*np.pi/8 * np.arange(numTaps+1))
h0 = 2 * np.real(hb * wn)       # Band1 Filter Coefficients
h1 = 2 * np.real(hb * wn**3)    # Band2 Filter Coefficients
h2 = 2 * np.real(hb * wn**5)    # Band3 Filter Coefficients
h3 = 2 * np.real(hb * wn**7)    # Band4 Filter Coefficients


The resulting frequency response for the Channelizer Filters is shown below:

# Demonstration of polyphase structure;
# This is identical to sig.lfilter(h, 1, signal)
def poly(h, signal):
filt0 = sig.lfilter(h[0::4], 1, signal[3::4])
filt1 = sig.lfilter(h[1::4], 1, signal[2::4])
filt2 = sig.lfilter(h[2::4], 1, signal[1::4])
filt3 = sig.lfilter(h[3::4], 1, signal[0::4])
return filt0 + filt1 + filt2 + filt3

nd = np.arange(len(signal)//4)

band1 = poly(h0, signal)
band2 = poly(h1, signal) * (-1)**nd    # flip spectrum for band2 and band4
band3 = poly(h2, signal)
band4 = poly(h3, signal) * (-1)**nd    # flip spectrum for band2 and band4


The resulting spectrum with the test waveform confirming proper operation is shown below, as noted above the amplitudes were scaled to increase in order of the band number to also confirm we aren't seeing an alias from an incorrect band:

Demonstration of Approach 2

The polyphase filter is constructed identical to the Band1 filter in Approach 1, without taking the real part to get the coefficients needed for the complex filter (this is essentially the analytic filter for Band1, as the real Band1 filter from above and its Hilbert Transform as the imaginary components).

# create complex band1 filter
numTaps = 199
hb = sig.firls(numTaps, [0, .1, .15, 1], [1, 1, 0, 0])
hb = np.append(hb,0) # make length divisible by 4 for polyphase implementation
wn = np.exp(1j*np.pi/8 * np.arange(numTaps+1))
h0c = hb * wn


The frequency response for this filter was shown in both real and complex form at the introduction of Approach 2. Below is the code demonstrating the polyphase implementation and output combining.

# Rotations are at filter outputs prior to combining, allowing then for real input but still need complex filter, take real of the result to get real output:

p0_out = sig.lfilter(h0c[0::4], 1, signal[3::4])
p1_out = sig.lfilter(h0c[1::4], 1, signal[2::4])
p2_out = sig.lfilter(h0c[2::4], 1, signal[1::4])
p3_out = sig.lfilter(h0c[3::4], 1, signal[0::4])

band1 = 2*np.real((p0_out + p1_out + p2_out + p3_out))

nd = np.arange(len(signal)//4)
# band 2 and 4 are spectrally inverted so multiplied by (-(-1)**nd)
band2 = 2*np.real((p0_out -1j*p1_out - p2_out + 1j*p3_out)*(-(-1)**nd) )
band3 = 2*np.real((p0_out +1j*p1_out - p2_out - 1j*p3_out))
band4 = 2*np.real((p0_out - p_out + p2_out - p3_out) * (-(-1)**nd) )


An additional test waveform was created to confirm proper phase relationships of all outputs and consistency with all implementations presented here. This waveform was just off of the low end of each Nyquist boundary by 31 MHz, 32 MHz, 33 MHz and 34 MHz respectively (so a sum of input tones at 31 MHz, 344.5 MHz, 658 MHz and 971.5 MHz). The result in frequency and time is plotted below and matched all implementations exactly.

Additional Notes:

When the intention is to implement structures in parallel that are identical to those operating at the full input rate, this would not be the recommended approach to take as more processing than necessary is involved and it introduces cross-over distortion as band edges are traversed. Instead, down-sample the input into four separate channels where each channel is every fourth sample, with each channel offset by one sample from the other channels. In this case, each channel contains information from all four bands, and subsequent processing would include samples from all channels to emulate the identical full rate structure. Please see this other similar post that explains this in more detail as well as provides the approach to create parallel filter implementations for this purpose:

Parallel decimating FIR

See fred harris' book "Multirate Signal Processing for Communication Systems" for further details on polyphase channelizers. The details in this post cover implementations for real inputs and outputs that can be real or complex.

• Comments are not for extended discussion; this conversation has been moved to chat.
– Peter K.
Jul 23, 2022 at 19:06