Skip to content
Snippets Groups Projects
Commit 3dff5517 authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Merge branch 'RTSD-268' into 'master' 

Resolve RTSD-268

Closes RTSD-268

See merge request !406
parents 5f4c30ae 92f2013e
Branches
No related tags found
1 merge request!406Resolve RTSD-268
Pipeline #81337 passed
Stage: linting
Show whitespace changes
Inline Side-by-side
Showing
with 1159 additions and 332 deletions
applications/lofar2/model/pfb_os/dsp.py
+ 221
− 65
View file @ 3dff5517
......@@ -117,6 +117,77 @@ def impulse_at_zero_crossing(x):
return np.concatenate((np.array([0]), diff))
###############################################################################
# Windowed sinc filter design
###############################################################################
def raised_cosine_response(Ntaps, Nsps, beta):
"""Generate a raised cosine (RC) FIR filter impulse response.
Input:
. Ntaps : FIR filter length
. Nsps: symbol period Tsymbol in number of samples per symbol
. beta : Roll off factor in [0, 1.0], BW = (1 + beta) / Tsymbol, so:
- beta = 0.0: rectangular spectrum with BW = 1 / Tsymbol
- beta = 1.0: cosine spectrum with BW = 2 / Tsymbol
Return:
. hRc : impulse response of the raised cosine filter.
"""
# time axis
Tsymbol = Nsps # normalized sample period Ts = 1.0
tIndices = np.arange(Ntaps)
tCenter = (Ntaps - 1) / 2
t = tIndices - tCenter
# sinc term, can use array assignment because sinc(1 / 0) = 1
hRc = 1 / Tsymbol * np.sinc(t / Tsymbol) # np.sinc(x) = sin(pi x) / (pi x)
# apply cos term, use for loop instead of array assignment, to detect divide by 0
for tI in tIndices:
t = tI - tCenter
if np.abs(t) != Tsymbol / (2 * beta):
hRc[tI] *= np.cos(np.pi * beta * t / Tsymbol) / (1 - (2 * beta * t / Tsymbol)**2)
return hRc
def square_root_raised_cosine_response(Ntaps, Nsps, beta):
"""Generate a square root raised cosine (SRRC) FIR filter impulse response.
Reference:
. [HARRIS section 4.3]
. https://en.wikipedia.org/wiki/Root-raised-cosine_filter
Input:
. Ntaps : FIR filter length
. Nsps: symbol period Tsymbol in number of samples per symbol
. beta : Roll off factor in [0, 1.0]
Return:
. hSrRc : impulse response of the square root raised cosine filter.
"""
# time axis
Tsymbol = Nsps # normalized sample period Ts = 1.0
tIndices = np.arange(Ntaps)
tCenter = (Ntaps - 1) / 2
t = tIndices - tCenter
# numerator term, using array assignment
hSrRc = 1 / Tsymbol * (np.cos(np.pi * (1 + beta) * t / Tsymbol) * 4 * beta * t / Tsymbol +
np.sin(np.pi * (1 - beta) * t / Tsymbol))
# apply denumerator term, use for loop instead of array assignment, to detect divide by 0
for tI in tIndices:
t = tI - tCenter
if t == 0.0:
hSrRc[tI] = 1 / Tsymbol * (1 + beta * (4 / np.pi - 1))
elif np.abs(t) == Tsymbol / (4 * beta):
hSrRc[tI] = 1 / Tsymbol * beta / np.sqrt(2) * \
((1 + 2 / np.pi) * np.sin(np.pi / (4 * beta)) +
(1 - 2 / np.pi) * np.cos(np.pi / (4 * beta)))
else:
hSrRc[tI] /= (1 - (4 * beta * t / Tsymbol)**2) * (np.pi * t / Tsymbol)
return hSrRc
###############################################################################
# FIR Filter design
###############################################################################
......@@ -169,9 +240,12 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0):
return h, f, HF
def prototype_fir_low_pass_filter(
Npoints=1024, Ntaps=16, Ncoefs=1024*16, hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, beta=1):
"""Derive FIR coefficients for prototype low pass filter using firls
def prototype_fir_low_pass_filter(method='firls',
Npoints=1024, Ntaps=16, Ncoefs=1024*16,
hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, beta=1, fs=1.0):
"""Derive FIR coefficients for prototype low pass filter
Use method 'firls' or 'remez'.
Default use LPF specification for LOFAR subband filter. For subband filter
BWbin = fs / Npoints and the number of subbands is Nsub = Npoints / 2.
......@@ -191,6 +265,7 @@ def prototype_fir_low_pass_filter(
Normalize DC gain to 1.0.
Input:
- method: FIR filter design method 'firls' or 'remez'
- Npoints: Number of points of the DFT in the filterbank, or multirate
factor Nup in upsampling or Ndown in downsampling
- Ntaps: Number of taps per polyphase.
......@@ -198,80 +273,134 @@ def prototype_fir_low_pass_filter(
coefs are zero
- hpFactor : Half power bandwidth of the filter relative to BWbin
- transitionFactor: transition bandwidth factor relative to fpass
- stopRippleFactor: stopband ripple facotr relative to pass band ripple
- beta: When beta > 0 then apply Kaiser window on FIR coefficients
- stopRippleFactor: stopband ripple factor relative to pass band ripple
- beta: When beta > 0 then additionally apply a Kaiser window on FIR
coefficients
- fs: sample frequency, for logging
Return:
- hFirls: FIR coefficients calculated with firls and optional Kaiser window
- hInterpolated: Interpolated FIR coefficients for requested Ncoefs
- h: FIR coefficients for requested Ncoefs
"""
# LPF specification for subband filter
BWbin = 1 / Npoints # bandwidth of one bin (is subband frequency)
BWbin = fs / Npoints # bandwidth of one bin (is subband frequency)
# Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0
BWpass = hpFactor * BWbin
fpass = BWpass / 2 # bin at DC: -fpass to +fpass
fpass = BWpass / 2 # pass band frequency of bin at DC: -fpass to +fpass
transitionBW = transitionFactor * fpass # transition bandwidth
fcutoff = 0 # not used, use hpFactor instead
cutoffGain = 0.5
fstop = fpass + transitionBW # stop band frequency
rippleWeights = [1, stopRippleFactor]
# Design subband filter
h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, beta, fs)
return h
# Initial FIR filter length
def design_fir_low_pass_filter(method,
Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1], beta=0, fs=1.0):
"""Derive FIR coefficients for prototype low pass filter
Use method 'firls' or 'remez', fs = 1.0
Derive symmetrical FIR coeffients for linear phase. For symmetrical
coefficient FIR filters the group delay is (Ncoefs - 1) / 2, so
choose Ncoefs is odd to have integers number of samples output delay.
Normalize DC gain to 1.0.
Input:
- method: FIR filter design method 'firls' or 'remez'
- Ncoefs: Actual number of prototype FIR filter coefficients, any extra
coefs are zero
- fpass: pass band frequency
- fstop: stop band frequency
- fs: sample frequency, 0 < fpass < fcutoff < fstop < fNyquist = fs / 2
- fcutoff: when fcutoff > 0, then define cutoff frequency point in transition band, fpass < fcutoff < fstop
- cutoffGain: normalized LPF gain at fcutoff
- rippleWeights: relative ripple factors for pass band, optional fcutoff, and stop band
- beta: When beta > 0, then additionally apply a Kaiser window on FIR
coefficients
Return:
- h: FIR coefficients for requested Ncoefs
"""
fNyquist = fs / 2
fsub = fpass * 2 # double sided LPF band width, when used as prototype for a BPF
# Initial FIR filter length N
# . Use interpolation of factor Q shorter filter to ensure the FIR filter design converges
# and to speed up calculation. N >> NFirlsMax = 1024 is not feasible.
# . The passband ripple and stopband attenuation depend on the transition bandwidth w_tb
# and the weight. Choose 0 ~< w_tb ~< 1.0 fpass, to ensure the FIR filter design converges
# and improve the passband ripple and stopband attenuation. A to large transition band
# also gives the design too much freedom and causes artifacts in the transition.
# . The passband ripple and stopband attenuation depend on the transition bandwidth
# and the rippleWeights. Choose transition band width ~< fpass, to ensure the FIR filter
# design converges and improve the passband ripple and stopband attenuation. A too large
# transition bandwidth also gives the design too much freedom and causes artifacts in
# the transition band.
# . scipy firls() defines fpass relative to fs, so can use fpass as cutoff frequency.
if method == 'firls':
NFirlsMax = 1024
Q = ceil_div(Ncoefs, NFirlsMax)
N = Ncoefs // Q
# If necessary -1 to make odd, because firls only supports odd number of FIR coefficients
if is_even(N):
N = N - 1
if method == 'remez':
NRemezMax = 1024
Q = ceil_div(Ncoefs, NRemezMax)
N = Ncoefs // Q
# Low pass transition band
f_pb = fpass * Q # pass band cut off frequency
w_tb = transitionFactor * fpass * Q # transition bandwidth
f_sb = f_pb + w_tb # stop band frequency
weight = [1, stopRippleFactor] # weight pass band ripple versus stop band ripple
print('> prototype_fir_low_pass_filter()')
print('Pass band specification')
print('. f_pb = ', str(f_pb))
print('. w_tb = ', str(w_tb))
print('. f_sb = ', str(f_sb))
hFirls = signal.firls(N, [0, f_pb, f_sb, 0.5], [1, 1, 0, 0], weight, fs=1.0)
# Apply Kaiser window with beta = 1 like in pfs_coeff_final.m, this improves the
# stopband attenuation near the transition band somewhat
f_pb = fpass * Q # pass band frequency
f_co = fcutoff * Q # optional cutoff frequency in transition band
f_sb = fstop * Q # stop band frequency
# Calculate initial FIR filter for length N
if method == 'firls':
if not fcutoff:
hFir = signal.firls(N, [0, f_pb, f_sb, fNyquist],
[1, 1, 0, 0], rippleWeights, fs=fs)
else:
hFir = signal.firls(N, [0, f_pb, f_co, f_co, f_sb, fNyquist],
[1, 1, cutoffGain, cutoffGain, 0, 0], rippleWeights, fs=fs)
if method == 'remez':
if not fcutoff:
hFir = signal.remez(N, [0, f_pb, f_sb, fNyquist],
[1, 0], rippleWeights, fs=fs)
else:
hFir = signal.remez(N, [0, f_pb, f_co, f_co, f_sb, fNyquist],
[1, cutoffGain, 0], rippleWeights, fs=fs)
# Additionally apply a Kaiser window, with beta = 1 like in pfs_coeff_final.m, this improves
# the stopband attenuation near the transition band somewhat
# . beta: 0 rect, 5 hamming, 6 hanning
if beta:
win = signal.windows.kaiser(N, beta=1)
hFirls *= win
win = signal.windows.kaiser(N, beta)
hFir *= win
# Normalize DC gain
dcSum = np.sum(hFirls)
hFirls *= 1.0 / dcSum
# Symmetrical FIR coeffients for linear phase
print('hFirls:')
print('. f_pb = %f' % f_pb)
print('. w_tb = %f' % w_tb)
print('. f_sb = %f' % f_sb)
print('. beta = %f' % beta)
print('. Q = %d' % Q)
print('. N = %d' % len(hFirls))
print('. DC sum = %f' % np.sum(hFirls))
print('. Symmetrical coefs = %s' % is_symmetrical(hFirls))
if len(hFirls) == Ncoefs:
h = hFirls
dcSum = np.sum(hFir)
hFir *= 1.0 / dcSum
if len(hFir) == Ncoefs:
h = hFir
else:
# Use Fourier interpolation to create final FIR filter coefs when
# Q > 1 or when Ncoefs is even
HFfirls = np.fft.fft(hFirls)
hInterpolated = fourier_interpolate(HFfirls, Ncoefs)
print('hInterpolated:')
print('. Ncoefs = %d' % len(hInterpolated))
print('. DC sum = %f' % np.sum(hInterpolated))
print('. Symmetrical coefs = %s' % is_symmetrical(hInterpolated))
h = hInterpolated
# Q > 1 or when N != Ncoefs
HFfir = np.fft.fft(hFir)
h = fourier_interpolate(HFfir, Ncoefs)
print('> design_fir_low_pass_filter()')
print('. Method = %s' % method)
print('. Q = %f' % Q)
print('. fsub = fpass * 2 = %f' % fsub)
print('. fpass = %f' % fpass)
if fcutoff:
print('. fcutoff = %f' % fcutoff)
print('. cutoffGain = %f' % cutoffGain)
print('. fstop = %f' % fstop)
print('. fNyquist = %f' % fNyquist)
print('. fs = %f' % fs)
print('. Ncoefs = %d' % len(h))
print('. DC sum = %f' % np.sum(h))
print('. Symmetrical coefs = %s' % is_symmetrical(h))
print('')
return h
......@@ -283,7 +412,7 @@ def fourier_interpolate(HFfilter, Ncoefs):
interpolation inserts Ncoefs - N zeros and then performs IFFT to get the
interpolated impulse response.
Use phase correction depenent on interpolation factor Q to have fractional
Use phase correction dependent on interpolation factor Q to have fractional
time shift of hInterpolated, to make it symmetrical. Similar as done in
pfs_coeff_final.m and pfir_coeff.m. Use upper = conj(lower), because that
is easier than using upper from HFfilter.
......@@ -914,18 +1043,23 @@ def resample(x, Nup, Ndown, coefs, verify=False): # interpolate and decimate by
# Plotting
###############################################################################
def plot_time_response(h, name='', markers=False):
def plot_time_response(h, title='', color='r', markers=False):
"""Plot time response (= impulse response, window, FIR filter coefficients).
Input:
. h: time response
. title: plot title
. color: curve color format character
. markers: when True plot time sample markers in curve
"""
if markers:
plt.plot(h, '-', h, 'o')
plt.plot(h, color + '-', h, color + 'o')
else:
plt.plot(h, color + '-')
if title:
plt.title(title)
else:
plt.plot(h, '-')
plt.title('Time response %s' % name)
plt.title('Time response')
plt.ylabel('Voltage')
plt.xlabel('Sample')
plt.grid(True)
......@@ -1091,6 +1225,28 @@ def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
plt.grid(True)
def plot_magnitude_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, voltLim=None):
"""Plot magnitude (= voltage) spectrum
Input:
. f: normalized frequency axis for HF (fs = 1)
. HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
. fmt: curve format string
. fs: sample frequency in Hz, scale f by fs, fs >= 1
"""
flabel = 'Frequency [fs = %f]' % fs
plt.plot(f * fs, np.abs(HF), fmt)
plt.title('Magnitude spectrum')
plt.ylabel('Voltage')
plt.xlabel(flabel)
if fLim:
plt.xlim(fLim)
if voltLim:
plt.ylim(voltLim)
plt.grid(True)
def plot_two_power_spectra(f1, HF1, name1, f2, HF2, name2, fs=1.0, fLim=None, dbLim=None, showRoll=False):
"""Plot two power spectra in same plot for comparison
......
applications/lofar2/model/pfb_os/dsp_study_erko.txt
+ 67
− 26
View file @ 3dff5517
......@@ -33,6 +33,7 @@
# * [JOS4] Spectral Audio Signal Processing, 2011
#
# * [WIKI] https://en.wikipedia.org/wiki/Bilinear_transform
# * [WIKI] https://en.wikipedia.org/wiki/Discrete_cosine_transform
# * [CHIPMUNK] https://github.com/chipmuenk : Python Frequency Design Analysis and DSP
# * [WHDLWHIZ] https://vhdlwhiz.com/articles/ : FIR filter design using DSP blocks
# * [BIQUAD]
......@@ -303,21 +304,30 @@
that the output blocks overlap, to keep the output and input sample rate
the same.
- Correlation equation [JOS1 7.2.5] (for comparison with convolution):
- Correlation
. Correlation is a measure of similarity between two function x(k) and y(k),
for different shifts (lag) in time. Autocorrelation can show the periodicity
of a signal, because then it has similarity for some k > 0.
. Difference between correlation and convolution is that convolution flips
one input, so corr(x, y) = conv(x, flip(y)). Hence if ne input is
symmetrical then correlation and convolution are the same.
* The purpose of convolution is to determine the output of a filter with
impulse response h.
* The purpose of correlation is to determine how much signal y is present
in x for different time delays (lags).
. Correlation equation [JOS1 7.2.5]:
N-1
xy(n) = sum conj(x(k)) y(n + k), time shift n is correlation lag
k=0
The cross power spectrum is [JOS1 8.4]:
. The cross power spectrum is [JOS1 8.4]:
HXY(w_k) = DFT_k(xy(n)) = 1 / N conj(X(w_k)) Y(w_k)
Correlation is a measure of similarity between two function x(k) and y(k),
for different shifts (lag) in time. Autocorrelation can show the periodicity
of a signal, because they it has similarity for some k > 0. To prove that
correlation can be expressed as convolution use a helper function
[WOLFSOUND]:
. To prove that correlation can be expressed as convolution use a helper
function [WOLFSOUND]:
xh[n] = sum_k x[n + k] h[k], with sum_k for k = -inf to +inf
= sum_k x[-(-n - k)] h[k]
......@@ -326,15 +336,16 @@
JOS4 7.2.4]: y(n) = sum_k x(n - k) h(k) = x(n) * h(k)
= x[-n] * h[n])[-n], again with x[-p] = x1[p], so correlation can be
calculated by convolving the time flipped x and then time flip
the result. Beware correlation is not commutative, so xh[n] !=
hx[n].
the result.
- FIR system identification from input-output measurements [JOS1 8.4.5,
PROAKIS 12]
. Convolution is commutative so x * y = y * x, but correlation is not
commutative, so xy != yx
. Use correlation for system identification from input-output measurements
[JOS1 8.4.5, PROAKIS 12]
y = h * x <==> H Y
y = h * x <=> H Y
xy = x cross y <==> conj(X) Y = conj(X) H Y = H |X|^2, so H = Rxy / Rxx
xy = x cross y <=> conj(X) Y = conj(X) H Y = H |X|^2, so H = Rxy / Rxx
4) Hilbert transform (HT) and analytic signal [LYONS 9]
......@@ -666,7 +677,7 @@
The DFT uses exp(-j w) and the IDFT uses exp(+j w), so applying IDFT on x(n)
will also result in a frequency domain representation.
- Matrix formulation, DFT as linear transformation [JOS1, PROAKIS 5.1.3]:
- Matrix formulation of DFT, is DFT as linear transformation [JOS1, PROAKIS 5.1.3]:
DFT:
XN = WN xN
......@@ -825,8 +836,38 @@
= [x_k * flip(w)](m), with x_k(n) = x(n) exp(-j w_k n)
9) Discrete Cosine Transform (DCT) [WIKI]
. The discrete transform of N samples implicitely assumes that the N samples
extend periodically, this causes discontinuities at the edges for the DFT.
. DCT type II has even symmetry on both sides, so block of N inputs extends
periodically like with DFT, but:
- for DCT it extends flipped to avoid a zero-th order discontinuity, the
slope (first order) is typically still discontinuous.
- and for type II it extends symmetrically half way between the end points
. The different types come from how the boudaries are defined.
. The DCT makes the transform converge faster than the DFT, because any
discontinuities in a function reduce the rate of convergence of the Fourier
series.
. DCT is used for image compression (like JPG), by keeping only few
coefficients of the transformed signal.
. DCT is equivalent to DFTs of roughly twice the length, operating on real
data with even symmetry.
. Discrete Sin Transform (DST) is equivalent to the imaginary parts of a
DFT of roughly twice the length, operating on real data with odd symmetry
(since the Fourier transform of a real and odd function is imaginary and
odd). The DCT is more common than the DST, because the with the DST the
boundaries typically still have discontinuities.
. The Modified DCT (MDCT) uses DCT type IV with overlap. With factor 2
overlap it is used for audio compression (like MP3). Factor 2 overlap
makes it computationally equivalent to DFT.
. The inverse of DCT-II is DCT-III multiplied by 2/N.
N-1
DCT-II: X(k) = sum x(n) cos(pi/n (n + 1/2) k), for k = 0, 1, ..., N-1
n=0
9) Multirate processing:
10) Multirate processing:
- Linear Time Variant (LTV) process, because it depends on when the
downsampling and upsampling start.
- Polyphase filtering ensures that only the values that remain are calculated,
......
applications/lofar2/model/pfb_os/filter_design_firls.ipynb
+ 37
− 25
View file @ 3dff5517
%% Cell type:markdown id:6e0a005d tags:
# Try firls FIR filter design method
Author: Eric Kooistra, nov 2023
Purpose:
* Practise DSP [1].
* Try firls least squares FIR filter design method for LPF.
* Try to reproduce LOFAR subband filter FIR coefficients using scipy instead of MATLAB.
MATLAB:
* The pfs_coeff_final.m from the Filter Task Force (FTF) in 2005 use fircls1 with r_pass and r_stop to define the ripple. In addition it post applies a Kaiser window with beta = 1 to make the filter attenuation a bit more deep near the transition.
* The pfir_coeff.m from Apertif also uses fircls1.
* Both use fircls1 with N = 1024 FIR coefficients and then Fourier interpolation to achieve Ncoefs = 1024 * 16 FIR coefficients. Both scripts can not exactly reproduce the actual LOFAR1 coefficients, therefore these are loaded from a file Coeffs16384Kaiser-quant.dat
Python (scipy.signal):
* The windowed sync method, firls leased squares method and remez method all yield comparable results, but firls and remez perform slightly better near the transition band. The firls and remez functions from scipy.signal use transition bandwidth and weights between pass and stop band to influence the transition region and ripple. For remez the ripple is constant in the pass band and stop band, for firls the ripple is largest near the band transition.
Conclusion:
* It is possible to design a good FIR filter using Python scipy. Possibly with some extra help of a filter design and analysis (FDA) tool like pyfda [2].
References:
1. dsp_study_erko, summary of DSP books
2. pyfda, dsp, at https://github.com/chipmuenk
%% Cell type:code id:3563bc63 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
```
%% Cell type:code id:f820b0ac tags:
``` python
import dsp
```
%% Cell type:code id:a131b5b6 tags:
``` python
import importlib
importlib.reload(dsp)
```
%% Output
<module 'dsp' from '/dop466_0/kooistra/git/hdl/applications/lofar2/model/pfb_os/dsp.py'>
%% Cell type:markdown id:2a467746 tags:
# 1 Least squares method
%% Cell type:code id:da2a98e9 tags:
``` python
# passband ripple (in dB);
r_pass_dB = 0.5;
r_pass = 10**(r_pass_dB / 20) - 1;
# stopband ripple (in dB);
r_stop_dB = -89;
r_stop = 10**(r_stop_dB / 20);
print('r_pass = %f' % r_pass)
print('r_stop = %f' % r_stop)
print('r_pass / r_stop = %f' % (r_pass / r_stop))
```
%% Output
r_pass = 0.059254
r_stop = 0.000035
r_pass / r_stop = 1669.996877
%% Cell type:code id:4b23f0c1 tags:
``` python
# LPF specification for LOFAR subband filter
Npoints = 1024 # = number of bins in fs, = DFT size
Ntaps = 16
Ncoefs = Npoints * Ntaps
hpFactor = 0.9
transitionFactor = 0.4
stopRippleFactor = 1000000
beta = 1
hPrototype = dsp.prototype_fir_low_pass_filter(
Npoints, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta)
fs = 1.0
#fs = 200e6 / hpFactor # to check fsub = 195312.5 Hz for LOFAR
hPrototype = dsp.prototype_fir_low_pass_filter('firls',
Npoints, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta, fs)
plt.plot(hPrototype, 'r', hPrototype - np.flip(hPrototype), 'r--')
plt.title('Impulse response')
plt.grid(True)
```
%% Output
f_pb = 0.00703125
w_tb = 0.0028125000000000003
f_sb = 0.00984375
hFirls:
. f_pb = 0.007031
. w_tb = 0.002813
. f_sb = 0.009844
. beta = 1.000000
. Q = 16
. N = 1023
. DC sum = 1.000000
. Symmetrical coefs = True
fNyquist = 0.5
fs = 1.0
hInterpolated.imag ~= 0
hInterpolated:
> design_fir_low_pass_filter()
. Method = firls
. Q = 16
. fsub = fpass * 2 = 8.789063e-04
. fpass = 4.394531e-04
. fstop = 6.152344e-04
. Ncoefs = 16384
. DC sum = 1.000000
. Symmetrical coefs = True
%% Cell type:code id:3ed56c18 tags:
``` python
fLim = (-2, 2)
#fLim = None
dbLim = (-150, 5)
#dbLim = None
h, f, HF = dsp.dtft(hPrototype)
dsp.plot_spectra(f, HF, Npoints, fLim, dbLim)
```
%% Output
%% Cell type:markdown id:e8acbe8f tags:
# 2 Compare firls filter and LOFAR subband filter
%% Cell type:code id:732899c1 tags:
``` python
fLim = (-2, 2)
dbLim = (-140, 5)
dbLim = (-140, 5) # view stop band attenuation
#dbLim = (-5, 5) # view pass band ripple
#plt.figure(0)
fs = Npoints
h, f, HF = dsp.dtft(hPrototype)
dsp.plot_power_spectrum(f, HF, 'b', fs, fLim, dbLim)
#plt.figure(1)
lofarCoefs = dsp.read_coefficients_file('../data/Coeffs16384Kaiser-quant.dat')
lofarCoefs /= np.sum(lofarCoefs)
fs = Npoints
h, f, HF = dsp.dtft(lofarCoefs)
dsp.plot_power_spectrum(f, HF, 'g', fs, fLim, dbLim)
plt.legend(['firls', 'lofarCoefs'])
```
%% Output
<matplotlib.legend.Legend at 0x7fe25ad489a0>
%% Cell type:code id:5f307eee tags:
``` python
```
......
applications/lofar2/model/pfb_os/filter_design_remez.ipynb
+ 67
− 120
View file @ 3dff5517
Source diff could not be displayed: it is too large. Options to address this: view the blob.
applications/lofar2/model/pfb_os/narrowband_noise_generator.ipynb
+ 48
− 61
View file @ 3dff5517
Source diff could not be displayed: it is too large. Options to address this: view the blob.
applications/lofar2/model/pfb_os/one_pfb.ipynb
+ 28
− 35
View file @ 3dff5517
%% Cell type:markdown id:c69a2eb8 tags:
 
# Try polyphase filterbank (PFB)
 
Author: Eric Kooistra, apr 2024
 
Purpose:
* Practise DSP [1].
* Try to reproduce LOFAR subband filterbank in Python instead of MATLAB [2]
 
References:
1. dsp_study_erko, summary of DSP books
2. https://git.astron.nl/rtsd/apertif_matlab/-/blob/master/matlab/one_pfb.m
 
%% Cell type:code id:02689e50 tags:
 
``` python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
```
 
%% Cell type:code id:65235f50 tags:
 
``` python
import dsp
```
 
%% Cell type:code id:acea4f0a tags:
 
``` python
import importlib
importlib.reload(dsp)
```
 
%% Output
 
<module 'dsp' from '/dop466_0/kooistra/git/hdl/applications/lofar2/model/pfb_os/dsp.py'>
 
%% Cell type:code id:a39e99a2 tags:
 
``` python
atol = 1e-8
```
 
%% Cell type:code id:55e0fe37 tags:
 
``` python
# Samples
fs = 1 # sample rate
ts = 1 / fs # sample period
```
 
%% Cell type:code id:3436bc2a tags:
 
``` python
# Subbands
Ndft = 16 # DFT size
Nsub = Ndft // 2 # number of subbands in fs / 2
fsub = fs / Ndft # subband frequency
tsub = 1 / fsub # subband period
```
 
%% Cell type:markdown id:a949b9ff tags:
 
# Waveform generator (WG)
 
%% Cell type:code id:e08b5ba2 tags:
 
``` python
# Time
Nsim = 1000 # number of subband periods to simulate
tsim = Nsim * tsub
t = np.arange(0, tsim, ts)
 
# Chirp signal
fstart = fsub * 0
fend = fsub * (Nsub + 1)
wgData = signal.chirp(t, f0=fstart, f1=fend, t1=tsim, method='linear')
 
plt.plot(t, wgData)
plt.title('Linear chirp, f = ' + str(fstart) + ' to f = ' + str(fend))
plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage')
#plt.xlim(700, 900)
plt.grid(True)
```
 
%% Output
 
 
%% Cell type:markdown id:acca4f19 tags:
 
# Prototype FIR low pass filter
 
%% Cell type:code id:ca6b8c9d tags:
 
``` python
# Select LPF type
firType = 'rect'
firType = 'firls'
#firType = 'firwin'
 
# Specifications
Ntaps = 8
Ncoefs = Ndft * Ntaps
if firType != 'rect':
hpFactor = 0.93
#hpFactor = 1.0
if firType == 'firls':
transitionFactor = 0.4
stopRippleFactor = 1000
beta = 1
```
 
%% Cell type:code id:181d6c07 tags:
 
``` python
# Use rectangular prototype FIR filter, for transparant response of DFT filterbank
if firType == 'rect':
hPrototype = np.zeros(Ncoefs)
hPrototype[0 : Ndft] = 1 / Ndft
 
print('hInterpolated:')
print('. Ncoefs = %d' % Ndft)
print('. DC sum = %f' % np.sum(hPrototype))
```
 
%% Cell type:code id:6108ff59 tags:
 
``` python
# Use subband prototype FIR filter, to minimize aliasing
if firType == 'firls':
hPrototype = dsp.prototype_fir_low_pass_filter(
Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta)
hPrototype = dsp.prototype_fir_low_pass_filter(firType,
Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta, fs)
```
 
%% Output
 
> prototype_fir_low_pass_filter()
Pass band specification
. f_pb = 0.0290625
. w_tb = 0.011625000000000002
. f_sb = 0.0406875
hFirls:
. f_pb = 0.029063
. w_tb = 0.011625
. f_sb = 0.040688
. beta = 1.000000
. Q = 1
. N = 127
. DC sum = 1.000000
. Symmetrical coefs = True
hInterpolated.imag ~= 0
hInterpolated:
> design_fir_low_pass_filter()
. Method = firls
. Q = 1.000000
. fsub = fpass * 2 = 0.058125
. fpass = 0.029063
. fstop = 0.040688
. fNyquist = 0.500000
. fs = 1.000000
. Ncoefs = 128
. DC sum = 1.000000
. Symmetrical coefs = True
 
%% Cell type:code id:3ec78793 tags:
 
``` python
# Use windowed sync prototype FIR filter, is Portnoff window
if firType == 'firwin':
BWbin = 1 / Ndft # bandwidth of one bin
# . Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0
BWpass = hpFactor * BWbin
fpass = BWpass / 2 # bin at DC: -fpass to +fpass
# scipy defines cutoff frequency relative to fnyquist = fs / 2 instead of fs.
fcutoff = 2 * fpass
hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann')
```
 
%% Cell type:code id:8e867248 tags:
 
``` python
# Plot impulse response
plt.figure(1)
plt.plot(hPrototype)
plt.title('Impulse response')
#plt.xlim((50, 60))
plt.grid(True)
 
# Plot transfer function
plt.figure(2)
fLim = (-2, 2)
#fLim = None
dbLim = (-140, 5)
h, f, HF = dsp.dtft(hPrototype)
dsp.plot_power_spectrum(f, HF, 'g', fs * Ndft, fLim, dbLim)
```
 
%% Output
 
 
 
%% Cell type:code id:eee415a1 tags:
 
``` python
# Polyphase FIR filter
pfsDown = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='down')
pfsUp = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='up')
```
 
%% Cell type:markdown id:d3628f9f tags:
 
# Analysis filterbank
 
%% Cell type:code id:9468907a tags:
 
``` python
pfsDownData = np.zeros((Nsim, Ndft))
dftData = np.zeros((Nsim, Ndft), dtype=np.complex_)
 
# pfir + dft
for b in range(Nsim):
# Input signal
inData = wgData[b * Ndft : (b + 1) * Ndft]
# Filtered signal
pfsDownData[b] = pfsDown.filter_block(inData)
# Frequency domain data
dftData[b] = np.fft.fft(pfsDownData[b])
#dftData[b] = Ndft * np.fft.ifft(pfsDownData[b])
```
 
%% Cell type:code id:4d046175 tags:
 
``` python
# Plot subband spectra
spectra = np.zeros((Nsim, Nsub))
for b in range(Nsim):
subbands = dftData[b][0:Nsub]
spectra[b] = dsp.pow_db(np.abs(subbands))
 
plt.figure(1)
plt.imshow(spectra, aspect='auto')
plt.title('Power spectra in dB')
plt.xlabel('f [fsub = ' + str(fsub) + ']')
plt.ylabel('t [tsub = ' + str(tsub) + ']')
dsp.plot_colorbar(plt)
 
plt.figure(2)
f = np.arange(0, Nsub)
Nb = 8 # number of spectra to plot
if Nb > Nsub:
Nb = Nsub
for b in range(0, Nsim, int(Nsim / Nb)):
subbands = dftData[b][0:Nsub]
dsp.plot_power_spectrum(f, subbands)
```
 
%% Output
 
 
 
%% Cell type:markdown id:f39d83a6 tags:
 
# Synthesis filterbank
 
%% Cell type:code id:c8d0560c tags:
 
``` python
# Using only IDFT
pfsUpData = np.zeros((Nsim, Ndft))
reconstructedData = np.zeros(Nsim * Ndft)
 
# idft
for b in range(Nsim):
# Frequency domain data
freqData = dftData[b]
timeData = Ndft * np.fft.ifft(freqData) # = pfsDownData
#timeData = np.fft.fft(freqData)
# Check that time signal is real
if np.max(np.abs(timeData.imag)) > atol:
plt.plot(timeData.imag)
break;
# Filtered signal
pfsUpData[b] = Ndft * pfsUp.filter_block(timeData.real)
#pfsUpData[b] = timeData.real
# Output signal
reconstructedData[b * Ndft : (b + 1) * Ndft] = pfsUpData[b]
```
 
%% Cell type:code id:c0573913 tags:
 
``` python
# Reconstructed output signal
if firType == 'rect':
groupDelay = 0 # because only first tap of Ntaps is none zero
else:
groupDelay = (Ncoefs - Ndft) / 2
groupDelay *= 2 # for hDown and hUp
print('groupDelay = ' + str(groupDelay) + ' samples')
outData = np.roll(reconstructedData, -int(groupDelay))
tDelta = wgData - outData
tDeltaMax = np.max(np.abs(tDelta))
print('max(abs(tDelta) = %e' % tDeltaMax)
if tDeltaMax > atol:
plt.plot(t, wgData, 'r')
plt.plot(t, outData, 'g')
#plt.plot(t, tDelta, 'b')
plt.title('Reconstructed delta signal')
plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage')
x = 500
plt.xlim((x, x + 800))
#plt.ylim([-0.05, 0.05])
plt.grid(True)
```
 
%% Output
 
groupDelay = 112.0 samples
max(abs(tDelta) = 1.090654e+00
 
 
%% Cell type:code id:aab1751c tags:
 
``` python
```
......
applications/lofar2/model/pfb_os/square_root_nyquist_filter.ipynb 0 → 100644
+ 691
− 0
View file @ 3dff5517
Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment