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

Improve prototype_fir_low_pass_filter() using design_fir_low_pass_filter().

parent f7318e7b
No related branches found
No related tags found
1 merge request!406Resolve RTSD-268
Pipeline #81336 passed
Stage: linting
Hide whitespace changes
Inline Side-by-side
applications/lofar2/model/pfb_os/one_pfb.ipynb
+ 28
− 35
View file @ 92f2013e
%% Cell type:markdown id:c69a2eb8 tags: %% Cell type:markdown id:c69a2eb8 tags:
   
# Try polyphase filterbank (PFB) # Try polyphase filterbank (PFB)
   
Author: Eric Kooistra, apr 2024 Author: Eric Kooistra, apr 2024
   
Purpose: Purpose:
* Practise DSP [1]. * Practise DSP [1].
* Try to reproduce LOFAR subband filterbank in Python instead of MATLAB [2] * Try to reproduce LOFAR subband filterbank in Python instead of MATLAB [2]
   
References: References:
1. dsp_study_erko, summary of DSP books 1. dsp_study_erko, summary of DSP books
2. https://git.astron.nl/rtsd/apertif_matlab/-/blob/master/matlab/one_pfb.m 2. https://git.astron.nl/rtsd/apertif_matlab/-/blob/master/matlab/one_pfb.m
   
%% Cell type:code id:02689e50 tags: %% Cell type:code id:02689e50 tags:
   
``` python ``` python
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy import signal from scipy import signal
``` ```
   
%% Cell type:code id:65235f50 tags: %% Cell type:code id:65235f50 tags:
   
``` python ``` python
import dsp import dsp
``` ```
   
%% Cell type:code id:acea4f0a tags: %% Cell type:code id:acea4f0a tags:
   
``` python ``` python
import importlib import importlib
importlib.reload(dsp) importlib.reload(dsp)
``` ```
   
%% Output %% Output
   
<module 'dsp' from '/dop466_0/kooistra/git/hdl/applications/lofar2/model/pfb_os/dsp.py'> <module 'dsp' from '/dop466_0/kooistra/git/hdl/applications/lofar2/model/pfb_os/dsp.py'>
   
%% Cell type:code id:a39e99a2 tags: %% Cell type:code id:a39e99a2 tags:
   
``` python ``` python
atol = 1e-8 atol = 1e-8
``` ```
   
%% Cell type:code id:55e0fe37 tags: %% Cell type:code id:55e0fe37 tags:
   
``` python ``` python
# Samples # Samples
fs = 1 # sample rate fs = 1 # sample rate
ts = 1 / fs # sample period ts = 1 / fs # sample period
``` ```
   
%% Cell type:code id:3436bc2a tags: %% Cell type:code id:3436bc2a tags:
   
``` python ``` python
# Subbands # Subbands
Ndft = 16 # DFT size Ndft = 16 # DFT size
Nsub = Ndft // 2 # number of subbands in fs / 2 Nsub = Ndft // 2 # number of subbands in fs / 2
fsub = fs / Ndft # subband frequency fsub = fs / Ndft # subband frequency
tsub = 1 / fsub # subband period tsub = 1 / fsub # subband period
``` ```
   
%% Cell type:markdown id:a949b9ff tags: %% Cell type:markdown id:a949b9ff tags:
   
# Waveform generator (WG) # Waveform generator (WG)
   
%% Cell type:code id:e08b5ba2 tags: %% Cell type:code id:e08b5ba2 tags:
   
``` python ``` python
# Time # Time
Nsim = 1000 # number of subband periods to simulate Nsim = 1000 # number of subband periods to simulate
tsim = Nsim * tsub tsim = Nsim * tsub
t = np.arange(0, tsim, ts) t = np.arange(0, tsim, ts)
   
# Chirp signal # Chirp signal
fstart = fsub * 0 fstart = fsub * 0
fend = fsub * (Nsub + 1) fend = fsub * (Nsub + 1)
wgData = signal.chirp(t, f0=fstart, f1=fend, t1=tsim, method='linear') wgData = signal.chirp(t, f0=fstart, f1=fend, t1=tsim, method='linear')
   
plt.plot(t, wgData) plt.plot(t, wgData)
plt.title('Linear chirp, f = ' + str(fstart) + ' to f = ' + str(fend)) plt.title('Linear chirp, f = ' + str(fstart) + ' to f = ' + str(fend))
plt.xlabel('t [ts = ' + str(ts) + ']') plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage') plt.ylabel('voltage')
#plt.xlim(700, 900) #plt.xlim(700, 900)
plt.grid(True) plt.grid(True)
``` ```
   
%% Output %% Output
   
   
%% Cell type:markdown id:acca4f19 tags: %% Cell type:markdown id:acca4f19 tags:
   
# Prototype FIR low pass filter # Prototype FIR low pass filter
   
%% Cell type:code id:ca6b8c9d tags: %% Cell type:code id:ca6b8c9d tags:
   
``` python ``` python
# Select LPF type # Select LPF type
firType = 'rect' firType = 'rect'
firType = 'firls' firType = 'firls'
#firType = 'firwin' #firType = 'firwin'
   
# Specifications # Specifications
Ntaps = 8 Ntaps = 8
Ncoefs = Ndft * Ntaps Ncoefs = Ndft * Ntaps
if firType != 'rect': if firType != 'rect':
hpFactor = 0.93 hpFactor = 0.93
#hpFactor = 1.0 #hpFactor = 1.0
if firType == 'firls': if firType == 'firls':
transitionFactor = 0.4 transitionFactor = 0.4
stopRippleFactor = 1000 stopRippleFactor = 1000
beta = 1 beta = 1
``` ```
   
%% Cell type:code id:181d6c07 tags: %% Cell type:code id:181d6c07 tags:
   
``` python ``` python
# Use rectangular prototype FIR filter, for transparant response of DFT filterbank # Use rectangular prototype FIR filter, for transparant response of DFT filterbank
if firType == 'rect': if firType == 'rect':
hPrototype = np.zeros(Ncoefs) hPrototype = np.zeros(Ncoefs)
hPrototype[0 : Ndft] = 1 / Ndft hPrototype[0 : Ndft] = 1 / Ndft
   
print('hInterpolated:') print('hInterpolated:')
print('. Ncoefs = %d' % Ndft) print('. Ncoefs = %d' % Ndft)
print('. DC sum = %f' % np.sum(hPrototype)) print('. DC sum = %f' % np.sum(hPrototype))
``` ```
   
%% Cell type:code id:6108ff59 tags: %% Cell type:code id:6108ff59 tags:
   
``` python ``` python
# Use subband prototype FIR filter, to minimize aliasing # Use subband prototype FIR filter, to minimize aliasing
if firType == 'firls': if firType == 'firls':
hPrototype = dsp.prototype_fir_low_pass_filter( hPrototype = dsp.prototype_fir_low_pass_filter(firType,
Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta) Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, beta, fs)
``` ```
   
%% Output %% 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.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 . Ncoefs = 128
. DC sum = 1.000000 . DC sum = 1.000000
. Symmetrical coefs = True . Symmetrical coefs = True
   
%% Cell type:code id:3ec78793 tags: %% Cell type:code id:3ec78793 tags:
   
``` python ``` python
# Use windowed sync prototype FIR filter, is Portnoff window # Use windowed sync prototype FIR filter, is Portnoff window
if firType == 'firwin': if firType == 'firwin':
BWbin = 1 / Ndft # bandwidth of one bin BWbin = 1 / Ndft # bandwidth of one bin
# . Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0 # . Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0
BWpass = hpFactor * BWbin BWpass = hpFactor * BWbin
fpass = BWpass / 2 # bin at DC: -fpass to +fpass fpass = BWpass / 2 # bin at DC: -fpass to +fpass
# scipy defines cutoff frequency relative to fnyquist = fs / 2 instead of fs. # scipy defines cutoff frequency relative to fnyquist = fs / 2 instead of fs.
fcutoff = 2 * fpass fcutoff = 2 * fpass
hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann') hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann')
``` ```
   
%% Cell type:code id:8e867248 tags: %% Cell type:code id:8e867248 tags:
   
``` python ``` python
# Plot impulse response # Plot impulse response
plt.figure(1) plt.figure(1)
plt.plot(hPrototype) plt.plot(hPrototype)
plt.title('Impulse response') plt.title('Impulse response')
#plt.xlim((50, 60)) #plt.xlim((50, 60))
plt.grid(True) plt.grid(True)
   
# Plot transfer function # Plot transfer function
plt.figure(2) plt.figure(2)
fLim = (-2, 2) fLim = (-2, 2)
#fLim = None #fLim = None
dbLim = (-140, 5) dbLim = (-140, 5)
h, f, HF = dsp.dtft(hPrototype) h, f, HF = dsp.dtft(hPrototype)
dsp.plot_power_spectrum(f, HF, 'g', fs * Ndft, fLim, dbLim) dsp.plot_power_spectrum(f, HF, 'g', fs * Ndft, fLim, dbLim)
``` ```
   
%% Output %% Output
   
   
   
%% Cell type:code id:eee415a1 tags: %% Cell type:code id:eee415a1 tags:
   
``` python ``` python
# Polyphase FIR filter # Polyphase FIR filter
pfsDown = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='down') pfsDown = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='down')
pfsUp = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='up') pfsUp = dsp.PolyPhaseFirFilterStructure(Ndft, hPrototype, commutator='up')
``` ```
   
%% Cell type:markdown id:d3628f9f tags: %% Cell type:markdown id:d3628f9f tags:
   
# Analysis filterbank # Analysis filterbank
   
%% Cell type:code id:9468907a tags: %% Cell type:code id:9468907a tags:
   
``` python ``` python
pfsDownData = np.zeros((Nsim, Ndft)) pfsDownData = np.zeros((Nsim, Ndft))
dftData = np.zeros((Nsim, Ndft), dtype=np.complex_) dftData = np.zeros((Nsim, Ndft), dtype=np.complex_)
   
# pfir + dft # pfir + dft
for b in range(Nsim): for b in range(Nsim):
# Input signal # Input signal
inData = wgData[b * Ndft : (b + 1) * Ndft] inData = wgData[b * Ndft : (b + 1) * Ndft]
# Filtered signal # Filtered signal
pfsDownData[b] = pfsDown.filter_block(inData) pfsDownData[b] = pfsDown.filter_block(inData)
# Frequency domain data # Frequency domain data
dftData[b] = np.fft.fft(pfsDownData[b]) dftData[b] = np.fft.fft(pfsDownData[b])
#dftData[b] = Ndft * np.fft.ifft(pfsDownData[b]) #dftData[b] = Ndft * np.fft.ifft(pfsDownData[b])
``` ```
   
%% Cell type:code id:4d046175 tags: %% Cell type:code id:4d046175 tags:
   
``` python ``` python
# Plot subband spectra # Plot subband spectra
spectra = np.zeros((Nsim, Nsub)) spectra = np.zeros((Nsim, Nsub))
for b in range(Nsim): for b in range(Nsim):
subbands = dftData[b][0:Nsub] subbands = dftData[b][0:Nsub]
spectra[b] = dsp.pow_db(np.abs(subbands)) spectra[b] = dsp.pow_db(np.abs(subbands))
   
plt.figure(1) plt.figure(1)
plt.imshow(spectra, aspect='auto') plt.imshow(spectra, aspect='auto')
plt.title('Power spectra in dB') plt.title('Power spectra in dB')
plt.xlabel('f [fsub = ' + str(fsub) + ']') plt.xlabel('f [fsub = ' + str(fsub) + ']')
plt.ylabel('t [tsub = ' + str(tsub) + ']') plt.ylabel('t [tsub = ' + str(tsub) + ']')
dsp.plot_colorbar(plt) dsp.plot_colorbar(plt)
   
plt.figure(2) plt.figure(2)
f = np.arange(0, Nsub) f = np.arange(0, Nsub)
Nb = 8 # number of spectra to plot Nb = 8 # number of spectra to plot
if Nb > Nsub: if Nb > Nsub:
Nb = Nsub Nb = Nsub
for b in range(0, Nsim, int(Nsim / Nb)): for b in range(0, Nsim, int(Nsim / Nb)):
subbands = dftData[b][0:Nsub] subbands = dftData[b][0:Nsub]
dsp.plot_power_spectrum(f, subbands) dsp.plot_power_spectrum(f, subbands)
``` ```
   
%% Output %% Output
   
   
   
%% Cell type:markdown id:f39d83a6 tags: %% Cell type:markdown id:f39d83a6 tags:
   
# Synthesis filterbank # Synthesis filterbank
   
%% Cell type:code id:c8d0560c tags: %% Cell type:code id:c8d0560c tags:
   
``` python ``` python
# Using only IDFT # Using only IDFT
pfsUpData = np.zeros((Nsim, Ndft)) pfsUpData = np.zeros((Nsim, Ndft))
reconstructedData = np.zeros(Nsim * Ndft) reconstructedData = np.zeros(Nsim * Ndft)
   
# idft # idft
for b in range(Nsim): for b in range(Nsim):
# Frequency domain data # Frequency domain data
freqData = dftData[b] freqData = dftData[b]
timeData = Ndft * np.fft.ifft(freqData) # = pfsDownData timeData = Ndft * np.fft.ifft(freqData) # = pfsDownData
#timeData = np.fft.fft(freqData) #timeData = np.fft.fft(freqData)
# Check that time signal is real # Check that time signal is real
if np.max(np.abs(timeData.imag)) > atol: if np.max(np.abs(timeData.imag)) > atol:
plt.plot(timeData.imag) plt.plot(timeData.imag)
break; break;
# Filtered signal # Filtered signal
pfsUpData[b] = Ndft * pfsUp.filter_block(timeData.real) pfsUpData[b] = Ndft * pfsUp.filter_block(timeData.real)
#pfsUpData[b] = timeData.real #pfsUpData[b] = timeData.real
# Output signal # Output signal
reconstructedData[b * Ndft : (b + 1) * Ndft] = pfsUpData[b] reconstructedData[b * Ndft : (b + 1) * Ndft] = pfsUpData[b]
``` ```
   
%% Cell type:code id:c0573913 tags: %% Cell type:code id:c0573913 tags:
   
``` python ``` python
# Reconstructed output signal # Reconstructed output signal
if firType == 'rect': if firType == 'rect':
groupDelay = 0 # because only first tap of Ntaps is none zero groupDelay = 0 # because only first tap of Ntaps is none zero
else: else:
groupDelay = (Ncoefs - Ndft) / 2 groupDelay = (Ncoefs - Ndft) / 2
groupDelay *= 2 # for hDown and hUp groupDelay *= 2 # for hDown and hUp
print('groupDelay = ' + str(groupDelay) + ' samples') print('groupDelay = ' + str(groupDelay) + ' samples')
outData = np.roll(reconstructedData, -int(groupDelay)) outData = np.roll(reconstructedData, -int(groupDelay))
tDelta = wgData - outData tDelta = wgData - outData
tDeltaMax = np.max(np.abs(tDelta)) tDeltaMax = np.max(np.abs(tDelta))
print('max(abs(tDelta) = %e' % tDeltaMax) print('max(abs(tDelta) = %e' % tDeltaMax)
if tDeltaMax > atol: if tDeltaMax > atol:
plt.plot(t, wgData, 'r') plt.plot(t, wgData, 'r')
plt.plot(t, outData, 'g') plt.plot(t, outData, 'g')
#plt.plot(t, tDelta, 'b') #plt.plot(t, tDelta, 'b')
plt.title('Reconstructed delta signal') plt.title('Reconstructed delta signal')
plt.xlabel('t [ts = ' + str(ts) + ']') plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage') plt.ylabel('voltage')
x = 500 x = 500
plt.xlim((x, x + 800)) plt.xlim((x, x + 800))
#plt.ylim([-0.05, 0.05]) #plt.ylim([-0.05, 0.05])
plt.grid(True) plt.grid(True)
``` ```
   
%% Output %% Output
   
groupDelay = 112.0 samples groupDelay = 112.0 samples
max(abs(tDelta) = 1.090654e+00 max(abs(tDelta) = 1.090654e+00
   
   
%% Cell type:code id:aab1751c tags: %% Cell type:code id:aab1751c tags:
   
``` python ``` python
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment