2
$\begingroup$

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 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 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.

filter Response

channelizer

In the above plots I have added 3 signals with frequencies 0.21, 0,51, and 0.81 and the 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?

changing the above code to

freq = [
    0.2,
    0.5,
    0.8,
]

produces:

channelizer

0.19, 0.21, etc works as expected.

$\endgroup$
1
  • 2
    $\begingroup$ Why don't you try with a simpler channelizer? maybe with less bands? would be easier to debug... $\endgroup$
    – Jdip
    Commented May 16 at 3:45

1 Answer 1

0
$\begingroup$

I also posted this question on dsprelated.com and got a great answer.

The signals are down-converted by the center frequency of each channel. So a tone at the center frequency will be a DC value with 0 frequency.

A related realization is that in the first plot I have 0.21, 0.51, and 0.81, frequencies channelized so that in each of the channels the resulting signal is a 0.01 Hz signal.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.