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

Use filter_block() with flipped parameter.

parent aaff79f9
No related branches found
No related tags found
1 merge request!419Resolve RTSD-265
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/pfb_os/pfb_reconstruction.ipynb
+ 29
− 38
View file @ 9e11c085
%% 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]
* Investigate reconstruction using critically sampled analysis PFB and synthesis PFB * Investigate reconstruction using critically sampled analysis PFB and synthesis PFB
   
Result: Result:
* Reconstruction with SNR = pow(input) / pow(input - reconctructed) > 18 dB does not seem feasible, for any firType. Increasing Ntaps also does not help. * Reconstruction with SNR = pow(input) / pow(input - reconctructed) > 18 dB does not seem feasible, for any firType. Increasing Ntaps also does not help.
* Stop band attenuation is important for loss due to aliasing from other than neighbour subband * Stop band attenuation is important for loss due to aliasing from other than neighbour subband
* Half power gain is important for all pass transfer function across neighbouring subbands * Half power gain is important for all pass transfer function across neighbouring subbands
* LOFAR subband filter achieves SNR = 14 dB with hpFactor = 0.96 to have fcutoff ~= sqrt(0.5) * LOFAR subband filter achieves SNR = 14 dB with hpFactor = 0.96 to have fcutoff ~= sqrt(0.5)
   
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
# Auto reload module when it is changed # Auto reload module when it is changed
%load_ext autoreload %load_ext autoreload
%autoreload 2 %autoreload 2
   
# Add rtdsp module path to $PYTHONPATH # Add rtdsp module path to $PYTHONPATH
import os import os
import sys import sys
module_path = os.path.abspath(os.path.join('../')) module_path = os.path.abspath(os.path.join('../'))
if module_path not in sys.path: if module_path not in sys.path:
sys.path.insert(0, module_path) sys.path.insert(0, module_path)
   
# Import rtdsp # Import rtdsp
from rtdsp.utilities import pow_db from rtdsp.utilities import pow_db
from rtdsp.firfilter import prototype_fir_low_pass_filter, design_fir_low_pass_filter, \ from rtdsp.firfilter import prototype_fir_low_pass_filter, design_fir_low_pass_filter, \
design_fir_low_pass_filter_adjust, \ design_fir_low_pass_filter_adjust, \
raised_cosine_response, square_root_raised_cosine_response, \ raised_cosine_response, square_root_raised_cosine_response, \
filterbank_impulse_response, filterbank_frequency_response filterbank_impulse_response, filterbank_frequency_response
from rtdsp.fourier import dtft, estimate_gain_at_frequency from rtdsp.fourier import dtft, estimate_gain_at_frequency
from rtdsp.multirate import PolyPhaseFirFilterStructure from rtdsp.multirate import PolyPhaseFirFilterStructure
from rtdsp.plotting import plot_power_spectrum, plot_magnitude_spectrum, plot_colorbar, \ from rtdsp.plotting import plot_power_spectrum, plot_magnitude_spectrum, plot_colorbar, \
plot_time_response plot_time_response
``` ```
   
%% Output
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
%% Cell type:code id:a39e99a2 tags: %% Cell type:code id:a39e99a2 tags:
   
``` python ``` python
atol = 1e-8 atol = 1e-8
``` ```
   
%% Cell type:code id:3436bc2a tags: %% Cell type:code id:3436bc2a tags:
   
``` python ``` python
# Filterbank # Filterbank
Ntaps = 8 #8 # number of taps per poly phase FIR filter Ntaps = 8 #8 # number of taps per poly phase FIR filter
Ndft = 16 #16 DFT size Ndft = 16 #16 DFT size
Ncoefs = Ndft * Ntaps # number of coefficients in prototype LPF Ncoefs = Ndft * Ntaps # number of coefficients in prototype LPF
   
# Select FFT or IFFT, end result is the same, but subband frequencies differ in sign # Select FFT or IFFT, end result is the same, but subband frequencies differ in sign
analysisFFT = True # analysis FFT --> synthesis IFFT as in LOFAR analysisFFT = True # analysis FFT --> synthesis IFFT as in LOFAR
# analysisFFT = False # analysis IFFT --> synthesis FFT as in [HARRIS 6, 7] # analysisFFT = False # analysis IFFT --> synthesis FFT as in [HARRIS 6, 7]
``` ```
   
%% 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 = 'firwin' firType = 'firwin'
#firType = 'firls' #firType = 'firls'
#firType = 'srRemezAdjust' #firType = 'srRemezAdjust'
#firType = 'srRemezCutoff' #firType = 'srRemezCutoff'
#firType = 'srCos' #firType = 'srCos'
``` ```
   
%% Cell type:code id:74eb2365 tags: %% Cell type:code id:74eb2365 tags:
   
``` python ``` python
# Select WG signal # Select WG signal
wgType = 'chirp' wgType = 'chirp'
wgType = 'noise' wgType = 'noise'
``` ```
   
%% 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:20ae2f4a tags: %% Cell type:code id:20ae2f4a tags:
   
``` python ``` python
# Subbands # Subbands
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:
   
# 1 Waveform generator (WG) # 1 Waveform generator (WG)
   
%% Cell type:code id:e08b5ba2 tags: %% Cell type:code id:e08b5ba2 tags:
   
``` python ``` python
if wgType == 'chirp': if wgType == 'chirp':
# Chirp signal # Chirp signal
chirpStartSub = 1 # start chirp subband chirpStartSub = 1 # start chirp subband
chirpNofSub = 10 # number of chirp subbands chirpNofSub = 10 # number of chirp subbands
#chirpNofSub = Nsub #chirpNofSub = Nsub
NSimSub = chirpNofSub * 10 # number of subband periods to simulate NSimSub = chirpNofSub * 10 # number of subband periods to simulate
NSimSamples = NSimSub * Ndft NSimSamples = NSimSub * Ndft
fstart = fsub * chirpStartSub fstart = fsub * chirpStartSub
fend = fsub * (chirpStartSub + chirpNofSub) fend = fsub * (chirpStartSub + chirpNofSub)
tsim = NSimSub * tsub tsim = NSimSub * tsub
t = np.arange(0, tsim, ts) # len(t) == NSimSamples t = np.arange(0, tsim, ts) # len(t) == NSimSamples
wgData = signal.chirp(t, f0=fstart, f1=fend, t1=tsim, method='linear') wgData = signal.chirp(t, f0=fstart, f1=fend, t1=tsim, method='linear')
   
print('Chirp:') print('Chirp:')
print('. chirpStartSub = %d' % chirpStartSub) print('. chirpStartSub = %d' % chirpStartSub)
print('. chirpNofSub = %d' % chirpNofSub) print('. chirpNofSub = %d' % chirpNofSub)
print('. NSimSub = %d' % NSimSub) print('. NSimSub = %d' % NSimSub)
print('. NSimSamples = %d' % NSimSamples) print('. NSimSamples = %d' % NSimSamples)
   
# Marker pulses at subband center frequencies # Marker pulses at subband center frequencies
wgMarkers = np.zeros(NSimSamples) wgMarkers = np.zeros(NSimSamples)
nofStep = (fend - fstart) / fsub nofStep = (fend - fstart) / fsub
tStep = tsim / nofStep tStep = tsim / nofStep
tMark = 0 tMark = 0
while tMark < tsim: while tMark < tsim:
tIndex = int(NSimSamples * tMark / tsim) tIndex = int(NSimSamples * tMark / tsim)
wgMarkers[tIndex] = 1 wgMarkers[tIndex] = 1
tMark += tStep tMark += tStep
   
plt.plot(t, wgData, 'b') plt.plot(t, wgData, 'b')
plt.plot(t, wgMarkers, 'r') plt.plot(t, wgMarkers, 'r')
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)
``` ```
   
%% Cell type:code id:bd4054d7 tags: %% Cell type:code id:bd4054d7 tags:
   
``` python ``` python
if wgType == 'noise': if wgType == 'noise':
NSimSub = Ntaps * 5 # number of subband periods to simulate NSimSub = Ntaps * 5 # number of subband periods to simulate
NSimSamples = NSimSub * Ndft NSimSamples = NSimSub * Ndft
tsim = NSimSub * tsub tsim = NSimSub * tsub
t = np.arange(0, tsim, ts) # len(t) == NSimSamples t = np.arange(0, tsim, ts) # len(t) == NSimSamples
   
# Broadband noise # Broadband noise
rng = np.random.default_rng(seed=7) rng = np.random.default_rng(seed=7)
mu = 0.0 mu = 0.0
sigma = 1.0 sigma = 1.0
wgData = rng.normal(mu, sigma, NSimSamples) wgData = rng.normal(mu, sigma, NSimSamples)
   
# Marker pulses at subband periods # Marker pulses at subband periods
wgMarkers = np.zeros(NSimSamples) wgMarkers = np.zeros(NSimSamples)
tStep = tsub * 3 tStep = tsub * 3
nofStep = int(tsim / tStep) nofStep = int(tsim / tStep)
tIndex = 0 tIndex = 0
while tIndex < NSimSamples: while tIndex < NSimSamples:
wgMarkers[tIndex] = 1 wgMarkers[tIndex] = 1
tIndex += Ndft tIndex += Ndft
   
plt.plot(t, wgData, 'b') plt.plot(t, wgData, 'b')
plt.plot(t, wgMarkers, 'r') plt.plot(t, wgMarkers, 'r')
plt.title('Broadband noise') plt.title('Broadband noise')
plt.xlabel('t [ts = ' + str(ts) + ']') plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage') plt.ylabel('voltage')
#plt.xlim(0, 100) #plt.xlim(0, 100)
plt.grid(True) plt.grid(True)
``` ```
   
%% Output %% Output
   
   
%% Cell type:markdown id:acca4f19 tags: %% Cell type:markdown id:acca4f19 tags:
   
# 2 Prototype FIR low pass filter # 2 Prototype FIR low pass filter
   
%% 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
``` ```
   
%% Cell type:code id:6108ff59 tags: %% Cell type:code id:6108ff59 tags:
   
``` python ``` python
# Use LOFAR subband prototype FIR filter, to minimize aliasing # Use LOFAR subband prototype FIR filter, to minimize aliasing
if firType == 'firls': if firType == 'firls':
hpFactor = 0.90 hpFactor = 0.90
transitionFactor = 0.4 transitionFactor = 0.4
hpFactor = 0.97 hpFactor = 0.97
transitionFactor = 0.2 transitionFactor = 0.2
stopRippleFactor = 1 # 1000000 stopRippleFactor = 1 # 1000000
kaiserBeta = 0 # 1 kaiserBeta = 0 # 1
hPrototype = prototype_fir_low_pass_filter('firls', hPrototype = prototype_fir_low_pass_filter('firls',
Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, kaiserBeta, fs) Ndft, Ntaps, Ncoefs, hpFactor, transitionFactor, stopRippleFactor, kaiserBeta, fs)
``` ```
   
%% 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
hpFactor = 1.11 hpFactor = 1.11
hpFactor = 1.0 hpFactor = 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
fcutoff = fpass fcutoff = fpass
hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann', fs=fs) hPrototype = signal.firwin(Ncoefs, fcutoff, window='hann', fs=fs)
``` ```
   
%% Cell type:code id:da37c5b9 tags: %% Cell type:code id:da37c5b9 tags:
   
``` python ``` python
# Use square root gain at fcutoff prototype FIR filter, to minimize aliasing # Use square root gain at fcutoff prototype FIR filter, to minimize aliasing
if firType == 'srRemezAdjust': if firType == 'srRemezAdjust':
cutoffBeta = 0.2 cutoffBeta = 0.2
fcutoff = fsub / 2 # center frequency of transition band fcutoff = fsub / 2 # center frequency of transition band
fpass = (1 - cutoffBeta) * fcutoff # pass band frequency fpass = (1 - cutoffBeta) * fcutoff # pass band frequency
fstop = (1 + cutoffBeta) * fcutoff # stop band frequency fstop = (1 + cutoffBeta) * fcutoff # stop band frequency
# Gain at fcutoff center frequency of the transition band # Gain at fcutoff center frequency of the transition band
cutoffGain = np.sqrt(0.5) cutoffGain = np.sqrt(0.5)
#cutoffGain = 0.5 #cutoffGain = 0.5
# weight pass band ripple versus stop band ripple # weight pass band ripple versus stop band ripple
rippleWeights = [1, 1] rippleWeights = [1, 1]
   
hPrototype = design_fir_low_pass_filter_adjust('remez', hPrototype = design_fir_low_pass_filter_adjust('remez',
Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, fs=fs, resolution=1e-4) Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, fs=fs, resolution=1e-4)
``` ```
   
%% Cell type:code id:4f11bd88 tags: %% Cell type:code id:4f11bd88 tags:
   
``` python ``` python
# Use square root gain at fcutoff prototype FIR filter, to minimize aliasing # Use square root gain at fcutoff prototype FIR filter, to minimize aliasing
if firType == 'srRemezCutoff': if firType == 'srRemezCutoff':
cutoffBeta = 0.2 cutoffBeta = 0.2
fcutoff = fsub / 2 # center frequency of transition band fcutoff = fsub / 2 # center frequency of transition band
fpass = (1 - cutoffBeta) * fcutoff # pass band frequency fpass = (1 - cutoffBeta) * fcutoff # pass band frequency
fstop = (1 + cutoffBeta) * fcutoff # stop band frequency fstop = (1 + cutoffBeta) * fcutoff # stop band frequency
# Gain at fcutoff center frequency of the transition band # Gain at fcutoff center frequency of the transition band
cutoffGain = np.sqrt(0.5) cutoffGain = np.sqrt(0.5)
# weight pass band ripple versus stop band ripple, and of zero width interval in the transition band # weight pass band ripple versus stop band ripple, and of zero width interval in the transition band
rippleWeights = [1, 1, 1] rippleWeights = [1, 1, 1]
   
hPrototype = design_fir_low_pass_filter('remez', hPrototype = design_fir_low_pass_filter('remez',
Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, fs=fs) Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, fs=fs)
``` ```
   
%% Cell type:code id:a3b429a8 tags: %% Cell type:code id:a3b429a8 tags:
   
``` python ``` python
if firType == 'srCos': if firType == 'srCos':
rollOffBeta = 0.2 rollOffBeta = 0.2
Nsps = Ndft Nsps = Ndft
hPrototype = square_root_raised_cosine_response(Ncoefs, Nsps, rollOffBeta) hPrototype = square_root_raised_cosine_response(Ncoefs, Nsps, rollOffBeta)
``` ```
   
%% 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 = dtft(hPrototype) h, f, HF = dtft(hPrototype)
plot_power_spectrum(f, HF, 'g', fs * Ndft, fLim, dbLim) plot_power_spectrum(f, HF, 'g', fs * Ndft, fLim, dbLim)
   
plt.figure(3) plt.figure(3)
voltLim = None voltLim = None
plot_magnitude_spectrum(f, HF, 'g', fs * Ndft, fLim, voltLim) plot_magnitude_spectrum(f, HF, 'g', fs * Ndft, fLim, voltLim)
   
# Log gain at fsub / 2 # Log gain at fsub / 2
fIndex, fValue, fGain = estimate_gain_at_frequency(f, HF, fsub / 2) fIndex, fValue, fGain = estimate_gain_at_frequency(f, HF, fsub / 2)
print('estimate_gain_at_frequency f = 0.5 [fsub]') print('estimate_gain_at_frequency f = 0.5 [fsub]')
print('. fIndex = %d' % fIndex) print('. fIndex = %d' % fIndex)
print('. fValue = %e [fsub]' % (fValue / fsub)) print('. fValue = %e [fsub]' % (fValue / fsub))
print('. fGain**2 = %e' % fGain**2) print('. fGain**2 = %e' % fGain**2)
``` ```
   
%% Output %% Output
   
estimate_gain_at_frequency f = 0.5 [fsub] estimate_gain_at_frequency f = 0.5 [fsub]
. fIndex = 1088 . fIndex = 1088
. fValue = 5.000000e-01 [fsub] . fValue = 5.000000e-01 [fsub]
. fGain**2 = 2.493174e-01 . fGain**2 = 2.493174e-01
   
   
   
   
%% Cell type:code id:55808f1c tags: %% Cell type:code id:55808f1c tags:
   
``` python ``` python
# Filterbank aggregate frequency response # Filterbank aggregate frequency response
h, f, HFprototype = dtft(hPrototype) h, f, HFprototype = dtft(hPrototype)
HFbank = filterbank_frequency_response(HFprototype, Ndft) HFbank = filterbank_frequency_response(HFprototype, Ndft)
   
fLim = None fLim = None
fLim = (0, 2) fLim = (0, 2)
dbLim = None dbLim = None
#dbLim = (-20, 5) #dbLim = (-20, 5)
#dbLim = (-0.2, 0.2) #dbLim = (-0.2, 0.2)
voltLim = None voltLim = None
voltLim = (0, 2) voltLim = (0, 2)
#voltLim = (0.9, 1.1) #voltLim = (0.9, 1.1)
plt.figure(1) plt.figure(1)
plot_power_spectrum(f, HFbank, 'g', fs * Ndft, fLim, dbLim) plot_power_spectrum(f, HFbank, 'g', fs * Ndft, fLim, dbLim)
plt.figure(2) plt.figure(2)
plot_magnitude_spectrum(f, HFbank, 'g', fs * Ndft, fLim, voltLim) plot_magnitude_spectrum(f, HFbank, 'g', fs * Ndft, fLim, voltLim)
``` ```
   
%% Output %% Output
   
   
   
%% Cell type:code id:6eb3b494 tags: %% Cell type:code id:6eb3b494 tags:
   
``` python ``` python
# Filterbank aggregate time response # Filterbank aggregate time response
hbank = filterbank_impulse_response(hPrototype, Ndft) hbank = filterbank_impulse_response(hPrototype, Ndft)
   
plot_time_response(hbank, title='Filterbank impulse response', color='r', markers=False) plot_time_response(hbank, title='Filterbank impulse response', color='r', markers=False)
``` ```
   
%% Output %% Output
   
   
%% Cell type:code id:eee415a1 tags: %% Cell type:code id:eee415a1 tags:
   
``` python ``` python
# Polyphase FIR filter, can use same for analysis PFB (down) and synthesis PFB (up) # Polyphase FIR filter, can use same for analysis PFB (down) and synthesis PFB (up)
pfs = PolyPhaseFirFilterStructure(Ndft, hPrototype) pfs = PolyPhaseFirFilterStructure(Ndft, hPrototype)
``` ```
   
%% Cell type:markdown id:d3628f9f tags: %% Cell type:markdown id:d3628f9f tags:
   
# 3 Analysis filterbank # 3 Analysis filterbank
   
%% Cell type:code id:9468907a tags: %% Cell type:code id:9468907a tags:
   
``` python ``` python
pfs.reset_poly_delays() pfs.reset_poly_delays()
pfsDownData = np.zeros((NSimSub, Ndft)) pfsDownData = np.zeros((NSimSub, Ndft))
dftData = np.zeros((NSimSub, Ndft), dtype=np.complex_) dftData = np.zeros((NSimSub, Ndft), dtype=np.complex_)
   
# pfir + dft # pfir + dft
for b in range(NSimSub): for b in range(NSimSub):
# Input signal # Input signal
inData = wgData[b * Ndft : (b + 1) * Ndft] inData = wgData[b * Ndft : (b + 1) * Ndft]
# Filter to downsample time signal # Filter to downsample time signal
pfsDownData[b] = pfs.filter_block(inData) pfsDownData[b] = pfs.filter_block(inData, flipped=False)
# Frequency domain data # Frequency domain data
# . The order of the pfs polyphases 0 : Ndft-1 fits the input order for FFT (and IFFT) # . The order of the pfs polyphases 0 : Ndft-1 fits the input order for FFT (and IFFT)
if analysisFFT: if analysisFFT:
dftData[b] = np.fft.fft(pfsDownData[b]) # LOFAR dftData[b] = np.fft.fft(pfsDownData[b]) # LOFAR
else: else:
dftData[b] = Ndft * np.fft.ifft(pfsDownData[b]) # [HARRIS Fig 6.14, 6.21] dftData[b] = Ndft * np.fft.ifft(pfsDownData[b]) # [HARRIS Fig 6.14, 6.21]
``` ```
   
%% Cell type:code id:4d046175 tags: %% Cell type:code id:4d046175 tags:
   
``` python ``` python
# Plot subband spectra # Plot subband spectra
spectra = np.zeros((NSimSub, Nsub)) spectra = np.zeros((NSimSub, Nsub))
for b in range(NSimSub): for b in range(NSimSub):
subbands = dftData[b][0:Nsub] subbands = dftData[b][0:Nsub]
spectra[b] = pow_db(np.abs(subbands)) spectra[b] = 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) + ']')
xstart = 0 xstart = 0
xend = Nsub xend = Nsub
if wgType == 'chirp': if wgType == 'chirp':
xend = min(chirpStartSub + chirpNofSub + 100, Nsub - 1) xend = min(chirpStartSub + chirpNofSub + 100, Nsub - 1)
plt.xlim((xstart, xend)) plt.xlim((xstart, xend))
plot_colorbar(plt) plot_colorbar(plt)
   
plt.figure(2) plt.figure(2)
fn = np.arange(0, Nsub) / Nsub / 2 fn = np.arange(0, Nsub) / Nsub / 2
Nb = 2 # number of spectra to plot Nb = 2 # number of spectra to plot
if wgType == 'chirp': if wgType == 'chirp':
if Nb > chirpNofSub: if Nb > chirpNofSub:
Nb = chirpNofSub Nb = chirpNofSub
for b in range(0, NSimSub, int(NSimSub / Nb)): for b in range(0, NSimSub, int(NSimSub / Nb)):
subbands = dftData[b][0:Nsub] subbands = dftData[b][0:Nsub]
plot_power_spectrum(fn, subbands, fs=Nsub * 2) plot_power_spectrum(fn, subbands, fs=Nsub * 2)
xstart = 0 xstart = 0
if wgType == 'chirp': if wgType == 'chirp':
xend = min(chirpStartSub + chirpNofSub + 10, Nsub -1) xend = min(chirpStartSub + chirpNofSub + 10, Nsub -1)
plt.xlim((xstart, xend)) plt.xlim((xstart, xend))
``` ```
   
%% Output %% Output
   
   
   
%% Cell type:markdown id:f39d83a6 tags: %% Cell type:markdown id:f39d83a6 tags:
   
# 4 Synthesis filterbank # 4 Synthesis filterbank
   
%% Cell type:code id:c8d0560c tags: %% Cell type:code id:c8d0560c tags:
   
``` python ``` python
pfs.reset_poly_delays() pfs.reset_poly_delays()
pfsUpData = np.zeros((NSimSub, Ndft)) pfsUpData = np.zeros((NSimSub, Ndft))
reconstructedData = np.zeros(NSimSub * Ndft) reconstructedData = np.zeros(NSimSub * Ndft)
   
# idft + pfir # idft + pfir
for b in range(NSimSub): for b in range(NSimSub):
# Frequency domain data # Frequency domain data
freqData = dftData[b] freqData = dftData[b]
# Time domain data: timeData = pfsDownData # Time domain data: timeData = pfsDownData
if analysisFFT: if analysisFFT:
timeData = Ndft * np.fft.ifft(freqData) # LOFAR timeData = Ndft * np.fft.ifft(freqData) # LOFAR
else: else:
timeData = np.fft.fft(freqData) # [HARRIS Fig 7.16] timeData = np.fft.fft(freqData) # [HARRIS Fig 7.16 has ifft]
# Check that domain data is real # Check that domain data 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;
# Filter to upsample time domain data # Filter to upsample time domain data
pfsUpData[b] = Ndft * pfs.filter_block(timeData.real) pfsUpData[b] = Ndft * pfs.filter_block(timeData.real, flipped=False)
# Reconstructed time signal is the pfs output from the Ndft polyphases. # Reconstructed time signal is the pfs output from the Ndft polyphases.
# . For upsampling the commutator selects the pfs polyphases 0 : Ndft-1, # . For upsampling the commutator selects the pfs polyphases 0 : Ndft-1,
# which fits the data output order in time # which fits the data output order in time
reconstructedData[b * Ndft : (b + 1) * Ndft] = pfsUpData[b] reconstructedData[b * Ndft : (b + 1) * Ndft] = pfsUpData[b]
``` ```
   
%% Cell type:code id:e01f671d tags: %% Cell type:code id:e01f671d tags:
   
``` python ``` python
# Group delay of output signal # Group delay of 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')
intGroupDelay = int(groupDelay) intGroupDelay = int(groupDelay)
``` ```
   
%% Output %% Output
   
groupDelay = 112.0 samples groupDelay = 112.0 samples
   
%% Cell type:code id:c0573913 tags: %% Cell type:code id:c0573913 tags:
   
``` python ``` python
# Compare reconstructed output with input # Compare reconstructed output with input
tt = t[0 : NSimSamples - intGroupDelay] tt = t[0 : NSimSamples - intGroupDelay]
inpData = wgData[0 : NSimSamples - intGroupDelay] inpData = wgData[0 : NSimSamples - intGroupDelay]
inpMarkers = wgMarkers[0 : NSimSamples - intGroupDelay] inpMarkers = wgMarkers[0 : NSimSamples - intGroupDelay]
outData = reconstructedData[intGroupDelay :] outData = reconstructedData[intGroupDelay :]
diffData = inpData - outData diffData = inpData - outData
diffDataMax = np.max(np.abs(diffData)) diffDataMax = np.max(np.abs(diffData))
SNRdb = pow_db(np.std(wgData) / np.std(diffData)) SNRdb = pow_db(np.std(wgData) / np.std(diffData))
# Expected SNR = 11.96 dB for Ntap = 8, Ndft = 16, 'firwin', hpFactor = 1.0, 'noise' # Expected SNR = 11.96 dB for Ntap = 8, Ndft = 16, 'firwin', hpFactor = 1.0, 'noise'
print('SNRdb = %.2f [dB]' % SNRdb) print('SNRdb = %.2f [dB]' % SNRdb)
if SNRdb < 100: if SNRdb < 100:
plt.figure(1, figsize=(11, 8)) plt.figure(1, figsize=(11, 8))
plt.title('Reconstructed signal for %s' % wgType) plt.title('Reconstructed signal for %s' % wgType)
plt.plot(tt, inpData, 'r') plt.plot(tt, inpData, 'r')
plt.plot(tt, outData, 'g') plt.plot(tt, outData, 'g')
plt.plot(tt, inpMarkers, 'b--') plt.plot(tt, inpMarkers, 'b--')
# use loc='lower center' to avoid time UserWarning on default loc='best' in case of large len(t) # use loc='lower center' to avoid time UserWarning on default loc='best' in case of large len(t)
plt.legend(['wgData', 'reconstructedData', 'wgMarkers'], loc='lower center') plt.legend(['wgData', 'reconstructedData', 'wgMarkers'], loc='lower center')
x = tStep * 1 x = tStep * 1
dx = tStep * 0.1 dx = tStep * 0.1
dStep = tStep * 1 dStep = tStep * 1
plt.xlim((x - dx, x + dx + dStep)) plt.xlim((x - dx, x + dx + dStep))
plt.xlabel('t [ts = ' + str(ts) + ']') plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage') plt.ylabel('voltage')
plt.grid(True) plt.grid(True)
   
plt.figure(2, figsize=(11, 8)) plt.figure(2, figsize=(11, 8))
plt.title('Difference signal for %s' % wgType) plt.title('Difference signal for %s' % wgType)
plt.plot(tt, diffData, 'm') plt.plot(tt, diffData, 'm')
plt.plot(tt, inpMarkers, 'b--') plt.plot(tt, inpMarkers, 'b--')
# use loc='lower center' to avoid time UserWarning on default loc='best' in case of large len(t) # use loc='lower center' to avoid time UserWarning on default loc='best' in case of large len(t)
plt.legend(['diffData', 'wgMarkers'], loc='lower center') plt.legend(['diffData', 'wgMarkers'], loc='lower center')
x = tStep * 1 x = tStep * 1
dx = tStep * 0.1 dx = tStep * 0.1
dStep = tStep * 1 dStep = tStep * 1
#plt.xlim((x - dx, x + dx + dStep)) #plt.xlim((x - dx, x + dx + dStep))
plt.xlabel('t [ts = ' + str(ts) + ']') plt.xlabel('t [ts = ' + str(ts) + ']')
plt.ylabel('voltage') plt.ylabel('voltage')
plt.grid(True) plt.grid(True)
``` ```
   
%% Output %% Output
   
SNRdb = 11.96 [dB] SNRdb = 11.96 [dB]
   
   
   
%% Cell type:code id:8bd92167 tags: %% Cell type:code id:8bd92167 tags:
   
``` python ``` python
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment