## Maximally Decimated Polyphase Channelizer Help

Started by 3 months ago3 replieslatest reply 3 months ago183 views
#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.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.

[ - ]