#python #Polyphase #Channelizer

I'm hoping someone can help me with this maximally decimated polyphase channelizer I've coded up. This is based on the fitler_ten_a code from the book "Multirate Signal Processing for Communication Systems", https://www.routledge.com/Multirate-Signal-Processing-for-Communication-Systems/Harris/p/book/9788770222105.

Here is my code. Hopefully someone can offer some critiques and point out anything I may have missed

# %%
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal


# %%

def dbscale(data):
    v0 = np.abs(data) / np.max(np.abs(data))
    v0 = 20.0 * np.log10(v0)
    return v0


def polyphase_filter_response(polyphase_filter):
    N = polyphase_filter.size * 100
    w = np.linspace(0, 1, N)
    p = 2 * np.pi * np.cumsum(w)
    data = np.exp(1j * p)

    channelizer_out = compute_channelizer(data, polyphase_filter)

    a = 20 * np.log10(np.abs(channelizer_out))
    w = np.linspace(0, 1, a.shape[1])

    plt.figure()
    for p in range(a.shape[0]):
        plt.plot(w, a[p], label=f'{p}')
    plt.legend()
    plt.title('Digital channelizer frequency response')
    plt.xlabel('Frequency [rad/sample]')
    plt.ylabel('Amplitude [dB]')
    plt.title('Polyphase Filter Response')


def compute_channelizer(data, polyphase_filter):
    num_data_samples = data.size
    num_channels, coef_per_stage = polyphase_filter.shape
    reg = np.zeros((num_channels, coef_per_stage), dtype=np.complex128)
    channelizer = np.zeros((num_channels, num_data_samples // num_channels), dtype=np.complex128)
    num_loops = 0
    for nn in range(0, num_data_samples, num_channels):
        reg[:, 1:] = reg[:, 0:-1]
        reg[:, 0] = data[nn:nn + num_channels]
        filter_output = np.zeros(num_channels, dtype=np.complex128)
        for mm in range(num_channels):
            # convolve the input with the filter and put it in reverse order in
            # filter output so the final stage is an ifft
            filter_output[num_channels - mm - 1] = reg[mm, ::-1].dot(polyphase_filter[mm])
        channelizer[:, num_loops] = np.fft.ifft(filter_output, norm='forward')
        num_loops += 1
    return channelizer

# %% setup


# for some reason frequencies like 0.2, 0.3, ... produce a DC value in the
# channelizer output. I assume this is some artifact from the sampling. This
# result is seems to be also present in the original matlab code.
freq = [
    0.21,
    0.51,
    0.81,
]
amp = [
    0.2,
    0.5,
    0.8,
]# sampling is 1 Hz in this example
fs = 1
ts = 1 / fs
num_channels = 10
coef_per_stage = 17

num_data_samples = 1200
t = np.arange(num_data_samples) * ts
t_channelizer = np.arange(num_data_samples // num_channels) * ts * num_channels
in_signal = np.zeros(num_data_samples, dtype=np.complex128)
for i, f in enumerate(freq):
    in_signal += amp[i] * np.exp(1j * 2 * np.pi * f * t)

# %% Filter Setup
order = num_channels * coef_per_stage
bands = np.array([0, 40, 60, 99, 100, 149, 150, 199, 200, 249, 250, 299, 300, 349, 350, 399, 400, 449, 450, 500]) / 500
gain = np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
weight = np.array([1, 100, 140, 180, 220, 260, 300, 340, 380, 420], dtype=np.float32)
filter_coef = signal.remez(order, bands, desired=gain, weight=weight, fs=2, grid_density=61)

# %% filter response
w, h = signal.freqz(filter_coef, 1, 2 ** 14)
h_abs = np.abs(h)
h_abs /= np.max(h_abs)
plt.figure()
plt.plot(fs / 2 * w / np.pi, 20 * np.log10(h_abs))
plt.grid(True)
plt.ylabel("Magnitude (dB)")
plt.xlabel("Frequency (Hz)")
plt.title("Filter Respone")

# %% Polyphase Filter - num_channels stages and coef_per_stage samples per stage

polyphase_filter = filter_coef.reshape((num_channels, coef_per_stage), order='F').copy()

# %% Polyphase Filer Response

polyphase_filter_response(polyphase_filter)

# %% apply the polyphase channelizer to the input

channelizer = compute_channelizer(in_signal, polyphase_filter)

# %% channelizer Plots plots

plt.figure(figsize=(24, 8))

plt.subplot(2, 6, 1)
plt.plot(t, in_signal.real, label='I')
plt.plot(t, in_signal.imag, label='Q')
plt.legend()
plt.title('Time Series Input')
in_sig_fft = dbscale(np.fft.fft(in_signal))
fft_freq = np.arange(in_signal.size) / (in_signal.size * ts)
plt.subplot(2, 6, 2)
plt.plot(fft_freq, in_sig_fft)
plt.ylim(-100, 10)
plt.title('Spectrum Input')

for i in range(num_channels):
    plt.subplot(2, 6, 3 + i)
    plt.plot(t_channelizer, channelizer[i].real, label='I')
    plt.plot(t_channelizer, channelizer[i].imag, label='Q')
    plt.ylim(-1.1, 1.1)
    plt.title(f'Channel {i}')

plt.show()

The Code seems to mostly work. Here are some of the plots from my code.

In the above plots I have added 3 signals with frequencies 0.21, 0,51, and 0.81 and they seem to be separated into their respective channels. However if I change the frequencies to 0.2, 0.5, 0.8 then the channel outputs are near constant. I don't understand. It may just be an artifact of the sampling rates?

after changed the following code I get this result.

freq = [
    0.2,
    0.5,
    0.8,
]

values like 0.19 and 0.21 work, but 0.2 gives the response above.

Isn't that exactly what you expect? If the signal is a pure sinewave, with frequency exactly equal to the channel center, then after down-conversion (channelization) that's at DC (0 Hz).

A 0 Hz sinewave is just a constant value, which is what you observe.

I think you are right! I just needed someone to spell it out for me. Thank you! I'm new to channelizers.

I'm quite sure I faced exactly the same confusion the first time I saw this phenomenon. Glad to help.