diff --git a/applications/lofar2/model/pfb_os/multirate_mixer.ipynb b/applications/lofar2/model/pfb_os/multirate_mixer.ipynb
index ebf514732ccbac8ac08ede5b35eae09a17b12e7d..1b47f011fe91aa3387415041cf64eb72cac04caf 100644
--- a/applications/lofar2/model/pfb_os/multirate_mixer.ipynb
+++ b/applications/lofar2/model/pfb_os/multirate_mixer.ipynb
@@ -20,7 +20,7 @@
},
{
"cell_type": "code",
- "execution_count": 389,
+ "execution_count": 1,
"id": "8043fa7b",
"metadata": {},
"outputs": [],
@@ -33,19 +33,10 @@
},
{
"cell_type": "code",
- "execution_count": 390,
+ "execution_count": 2,
"id": "fc530dbc",
"metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "The autoreload extension is already loaded. To reload it, use:\n",
- " %reload_ext autoreload\n"
- ]
- }
- ],
+ "outputs": [],
"source": [
"# Auto reload module when it is changed\n",
"%load_ext autoreload\n",
@@ -79,7 +70,7 @@
},
{
"cell_type": "code",
- "execution_count": 391,
+ "execution_count": 3,
"id": "5b37a1dc",
"metadata": {},
"outputs": [
@@ -148,7 +139,7 @@
},
{
"cell_type": "code",
- "execution_count": 392,
+ "execution_count": 4,
"id": "e5680c7b",
"metadata": {},
"outputs": [],
@@ -160,7 +151,7 @@
},
{
"cell_type": "code",
- "execution_count": 393,
+ "execution_count": 5,
"id": "74ca764f",
"metadata": {},
"outputs": [
@@ -189,7 +180,7 @@
},
{
"cell_type": "code",
- "execution_count": 394,
+ "execution_count": 6,
"id": "786af296",
"metadata": {},
"outputs": [],
@@ -220,7 +211,7 @@
},
{
"cell_type": "code",
- "execution_count": 395,
+ "execution_count": 7,
"id": "1bb76ada",
"metadata": {},
"outputs": [
@@ -252,7 +243,7 @@
},
{
"cell_type": "code",
- "execution_count": 396,
+ "execution_count": 8,
"id": "6ebc94aa",
"metadata": {},
"outputs": [
@@ -283,7 +274,7 @@
},
{
"cell_type": "code",
- "execution_count": 397,
+ "execution_count": 9,
"id": "3abeee86",
"metadata": {},
"outputs": [
@@ -339,7 +330,7 @@
},
{
"cell_type": "code",
- "execution_count": 398,
+ "execution_count": 10,
"id": "372445f4",
"metadata": {},
"outputs": [
@@ -371,7 +362,7 @@
},
{
"cell_type": "code",
- "execution_count": 399,
+ "execution_count": 11,
"id": "e74340e4",
"metadata": {},
"outputs": [],
@@ -390,7 +381,7 @@
},
{
"cell_type": "code",
- "execution_count": 400,
+ "execution_count": 12,
"id": "bd2add56",
"metadata": {},
"outputs": [],
@@ -415,17 +406,17 @@
},
{
"cell_type": "code",
- "execution_count": 401,
+ "execution_count": 13,
"id": "7106ad3f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
- "[<matplotlib.lines.Line2D at 0x7fd6133cae80>]"
+ "[<matplotlib.lines.Line2D at 0x7f14108ec160>]"
]
},
- "execution_count": 401,
+ "execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
@@ -465,7 +456,7 @@
},
{
"cell_type": "code",
- "execution_count": 402,
+ "execution_count": 14,
"id": "da53b25e",
"metadata": {},
"outputs": [],
@@ -477,17 +468,17 @@
},
{
"cell_type": "code",
- "execution_count": 403,
+ "execution_count": 15,
"id": "9acf0ec2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
- "<matplotlib.legend.Legend at 0x7fd613deacd0>"
+ "<matplotlib.legend.Legend at 0x7f1410acb760>"
]
},
- "execution_count": 403,
+ "execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
@@ -533,7 +524,7 @@
},
{
"cell_type": "code",
- "execution_count": 404,
+ "execution_count": 16,
"id": "b8036250",
"metadata": {},
"outputs": [
@@ -578,7 +569,7 @@
},
{
"cell_type": "code",
- "execution_count": 405,
+ "execution_count": 17,
"id": "0663df66",
"metadata": {},
"outputs": [],
@@ -594,7 +585,7 @@
},
{
"cell_type": "code",
- "execution_count": 406,
+ "execution_count": 18,
"id": "3a039428",
"metadata": {},
"outputs": [
@@ -638,7 +629,7 @@
},
{
"cell_type": "code",
- "execution_count": 407,
+ "execution_count": 19,
"id": "327236c2",
"metadata": {},
"outputs": [
@@ -675,7 +666,7 @@
},
{
"cell_type": "code",
- "execution_count": 408,
+ "execution_count": 20,
"id": "aefa8615",
"metadata": {},
"outputs": [
@@ -688,6 +679,7 @@
" . Nx = 373\n",
" . Nblocks = 32\n",
" . len(yc) = 32\n",
+ " . Ros = 1.3333333333333333\n",
" . Ndown = 12\n",
" . Ndft = 16\n",
" . k = 1\n",
@@ -698,10 +690,10 @@
{
"data": {
"text/plain": [
- "[<matplotlib.lines.Line2D at 0x7fd613329700>]"
+ "[<matplotlib.lines.Line2D at 0x7f140afb4d60>]"
]
},
- "execution_count": 408,
+ "execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
@@ -738,7 +730,7 @@
},
{
"cell_type": "code",
- "execution_count": 409,
+ "execution_count": 21,
"id": "f814c9b9",
"metadata": {},
"outputs": [
@@ -780,7 +772,7 @@
},
{
"cell_type": "code",
- "execution_count": 410,
+ "execution_count": 22,
"id": "dd7a9503",
"metadata": {},
"outputs": [],
@@ -817,7 +809,7 @@
},
{
"cell_type": "code",
- "execution_count": 411,
+ "execution_count": 23,
"id": "3cf6aa74",
"metadata": {},
"outputs": [],
@@ -828,7 +820,7 @@
},
{
"cell_type": "code",
- "execution_count": 412,
+ "execution_count": 24,
"id": "03da0b30",
"metadata": {},
"outputs": [
@@ -849,7 +841,7 @@
},
{
"cell_type": "code",
- "execution_count": 413,
+ "execution_count": 25,
"id": "a5c60be5",
"metadata": {},
"outputs": [
@@ -904,7 +896,7 @@
},
{
"cell_type": "code",
- "execution_count": 414,
+ "execution_count": 26,
"id": "7049249e",
"metadata": {},
"outputs": [
@@ -938,7 +930,7 @@
},
{
"cell_type": "code",
- "execution_count": 415,
+ "execution_count": 27,
"id": "64cc34f3",
"metadata": {},
"outputs": [],
@@ -954,7 +946,7 @@
},
{
"cell_type": "code",
- "execution_count": 416,
+ "execution_count": 28,
"id": "71a91beb",
"metadata": {},
"outputs": [
@@ -1003,7 +995,7 @@
},
{
"cell_type": "code",
- "execution_count": 417,
+ "execution_count": 29,
"id": "a3aae48a",
"metadata": {},
"outputs": [
@@ -1040,7 +1032,7 @@
},
{
"cell_type": "code",
- "execution_count": 418,
+ "execution_count": 30,
"id": "3c2c8ec5",
"metadata": {},
"outputs": [
@@ -1067,7 +1059,7 @@
},
{
"cell_type": "markdown",
- "id": "a9e76fc1",
+ "id": "ee9daf9f",
"metadata": {},
"source": [
"# 5 Compare with DFT filterbank"
@@ -1075,7 +1067,7 @@
},
{
"cell_type": "markdown",
- "id": "a6bce8dc",
+ "id": "d8bb436d",
"metadata": {},
"source": [
"# 5.1 Analysis"
@@ -1083,8 +1075,8 @@
},
{
"cell_type": "code",
- "execution_count": 419,
- "id": "45850a76",
+ "execution_count": 31,
+ "id": "b0f997d3",
"metadata": {},
"outputs": [
{
@@ -1092,17 +1084,19 @@
"output_type": "stream",
"text": [
"> analysis_dft_filterbank():\n",
- " . len(x) = 384\n",
- " . Nx = 373\n",
- " . Nblocks = 32\n",
- " . Ndown = 12\n",
- " . Ndft = 16\n",
+ " . len(x) = 384\n",
+ " . Nx = 373\n",
+ " . Nblocks = 32\n",
+ " . Ros = 1.3333333333333333\n",
+ " . Ndown = 12\n",
+ " . Ndft = 16\n",
+ " . commutator = cw\n",
"\n"
]
}
],
"source": [
- "Yc = analysis_dft_filterbank(xData, Ndown, Ndft, hPrototype, verbosity=1)\n",
+ "Yc = analysis_dft_filterbank(xData, Ndown, Ndft, hPrototype, commutator='cw', verbosity=1)\n",
"yDownBin = Yc[kLo]\n",
"\n",
"result = np.all(np.isclose(yDown, yDownBin))\n",
@@ -1117,7 +1111,7 @@
},
{
"cell_type": "markdown",
- "id": "772350cd",
+ "id": "d3f49b8e",
"metadata": {},
"source": [
"# 5.2 Synthesis"
@@ -1126,7 +1120,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "9a027306",
+ "id": "47d5bf5b",
"metadata": {},
"outputs": [],
"source": []
diff --git a/applications/lofar2/model/rtdsp/multirate.py b/applications/lofar2/model/rtdsp/multirate.py
index 2388ea34388ea18bae5fa7ba092005c438042650..7c3d14862d2a630fe637f76fcaf960d5166a9ad9 100644
--- a/applications/lofar2/model/rtdsp/multirate.py
+++ b/applications/lofar2/model/rtdsp/multirate.py
@@ -28,10 +28,10 @@
#
# References:
# [1] dsp_study_erko.txt
-# [2] http://localhost:8888/notebooks/pfb_os/up_downsampling.ipynb
+# [2] pfb_os/up_downsampling.ipynb, pfb_os/multirate_mixer.ipynb
#
# Books:
-# . HARRIS 6 title Figure
+# . HARRIS 6
# . LYONS
# . CROCHIERE
@@ -54,6 +54,17 @@ def ones(shape, cmplx=False):
return zeros(shape, cmplx) + 1
+def fold(x):
+ """Circular reverse sequence around index 0.
+
+ x[n] --> flip(x) --> x[N - n - 1]
+ x[n] --> fold(x) -> x[-n] = x[N - n] = roll(flip(x), 1)
+
+ The index is modulo N, so x[N - 0] = x[N] = x[0].
+ """
+ return np.roll(np.flip(x), 1)
+
+
def down(x, D, phase=0):
"""Downsample x[n] by factor D, xD[m] = x[m D], m = n // D
@@ -248,7 +259,7 @@ class PolyPhaseFirFilterStructure:
and <= Nphases for downsampling.
. flipped: False then inData order is inData[n] with newest sample last
and still needs to be flipped. If True then the inData is already
- flipped in time, so with newest sample first.
+ flipped in time, so with newest sample at inData[0].
Return:
. pfsData: block of polyphase FIR filtered output data for Nphases, with
pfsData[p] and p = 0:Nphases-1 from top to bottom.
@@ -308,7 +319,7 @@ def polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros):
"""
Lx = len(x)
Nxp = (Nzeros + Lx) // Ndown # Number of samples per polyphase
- Nx = 1 + Ndown * (Nxp - 1) # Used number of samples from x
+ Nx = Ndown * Nxp - Nzeros # Used number of samples from x
# Load x into polyX with Ndown rows = polyphases
# . prepend x with Nzeros zeros
@@ -817,6 +828,8 @@ def non_maximal_downsample_bpf(x, Ndown, k, Ndft, coefs, verbosity=1):
. yc: Downsampled and downconverted output signal yc[m], m = n // D for
bin k. Complex baseband signal.
"""
+ Ros = Ndft / Ndown
+
# Prepend x with Ndown - 1 zeros, and represent x in Nblocks of Ndown
# samples
Nzeros = Ndown - 1
@@ -849,6 +862,7 @@ def non_maximal_downsample_bpf(x, Ndown, k, Ndft, coefs, verbosity=1):
print(' . Nx =', str(Nx))
print(' . Nblocks =', str(Nblocks))
print(' . len(yc) =', str(len(yc))) # = Nblocks
+ print(' . Ros =', str(Ros))
print(' . Ndown =', str(Ndown))
print(' . Ndft =', str(Ndft))
print(' . k =', str(k))
@@ -937,6 +951,7 @@ def non_maximal_upsample_bpf(xBase, Nup, k, Ndft, coefs, verbosity=1):
print(' . len(xBase) =', str(len(xBase)))
print(' . Nblocks =', str(Nblocks))
print(' . len(yc) =', str(len(yc))) # = Nblocks
+ print(' . Ros =', str(Ros))
print(' . Nup =', str(Nup))
print(' . Ndft =', str(Ndft))
print(' . k =', str(k))
@@ -944,7 +959,7 @@ def non_maximal_upsample_bpf(xBase, Nup, k, Ndft, coefs, verbosity=1):
return yc
-def analysis_dft_filterbank(x, Ndown, Ndft, coefs, verbosity=1):
+def analysis_dft_filterbank(x, Ndown, Ndft, coefs, commutator, verbosity=1):
"""DFT filterbank with Ros = Ndft / Ndown.
Implements WOLA structure for DFT filterbank [CROCHIERE 7.2.5]. Key steps:
@@ -967,23 +982,29 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, verbosity=1):
exp(-j 2pi k m M / K) of the DFT output is equivalent to a circular time
shift by l = (m M) % K samples of the DFT input. Circular time shift DFT
theorem:
-
x[n] <= DFT => X(k), then
x[(n - l) % K] <= DFT => X(k) exp(-j 2pi k l / K)
+ - For counter clockwise commutator use Ndft * IDFT(pfsData).
+ For clockwise commutator use DFT(fold(pfsData)), because
+ DFT(x) = Ndft * IDFT(fold(x)).
Input:
. x: Input signal x[n]
. Ndown: downsample factor and input block size
. Ndft: DFT size, number of polyphases in PFS FIR filter
. coefs: prototype LPF FIR filter coefficients for anti aliasing BPF
+ . commutator: 'ccw' to use IDFT or 'cw' to use DFT
- verbosity: when > 0 print() status, else no print()
Return:
. Yc: Downsampled and downconverted output signal Yc[k, m], m = n // M for
all Ndft bins k. Complex baseband signals.
"""
+ Ros = Ndft / Ndown
+
# Prepend x with Ndown - 1 zeros, and represent x in Nblocks of Ndown
# samples
Nzeros = Ndown - 1
+ # Nzeros = 0
xBlocks, Nx, Nblocks = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros)
# Prepare output
@@ -998,20 +1019,31 @@ def analysis_dft_filterbank(x, Ndown, Ndft, coefs, verbosity=1):
inData = xBlocks[:, b]
pfsData = pfs.filter_block(inData, flipped=True)
+ # Default with 'ccw' keep the pfsData to apply IDFT
+ if commutator == 'cw':
+ # For 'cw' fold (= circular flip) the pfsData to apply the DFT
+ pfsData = fold(pfsData)
+
# Apply time (m) and bin (k) dependend phase shift W_K^(kmM) by circular
# time shift of DFT input
tShift = (b * Ndown) % Ndft
- pfsShifted = np.roll(pfsData, -tShift)
- # IDFT
- Yc[:, b] = Ndft * np.fft.ifft(pfsShifted)
+ # For 'ccw' apply IDFT, for 'cw' apply DFT
+ if commutator == 'cw': # DFT
+ pfsShifted = np.roll(pfsData, tShift)
+ Yc[:, b] = np.fft.fft(pfsShifted)
+ else: # 'ccw', IDFT
+ pfsShifted = np.roll(pfsData, -tShift)
+ Yc[:, b] = np.fft.ifft(pfsShifted) * Ndft
if verbosity:
print('> analysis_dft_filterbank():')
- print(' . len(x) =', str(len(x)))
- print(' . Nx =', str(Nx))
- print(' . Nblocks =', str(Nblocks))
- print(' . Ndown =', str(Ndown))
- print(' . Ndft =', str(Ndft))
+ print(' . len(x) =', str(len(x)))
+ print(' . Nx =', str(Nx))
+ print(' . Nblocks =', str(Nblocks))
+ print(' . Ros =', str(Ros))
+ print(' . Ndown =', str(Ndown))
+ print(' . Ndft =', str(Ndft))
+ print(' . commutator =', commutator)
print('')
return Yc