diff --git a/applications/lofar2/model/pfb_os/multirate_mixer.ipynb b/applications/lofar2/model/pfb_os/multirate_mixer.ipynb index 01eede15225a2600b1e14bbb2603155d8d328438..8fe5820fd4df1663e83c85e08a5f599bcdf767dd 100644 --- a/applications/lofar2/model/pfb_os/multirate_mixer.ipynb +++ b/applications/lofar2/model/pfb_os/multirate_mixer.ipynb @@ -7,7 +7,7 @@ "source": [ "# Multirate mixer\n", "\n", - "Author: Eric Kooistra, May 2024\n", + "Author: Eric Kooistra, May - Dec 2024\n", "\n", "Purpose:\n", "* Practise DSP [1].\n", @@ -18,20 +18,31 @@ "\n", "For a real input the mixer can either downconvert the positive or the negative frequeny band. Convention is to downconvert the positive frequency band [2, 3]. This leads to using IDFT in the analysis DFT filterbank and IDFT in the synthesis filterbank [2, 3]. By using IDFT(x) = DFT(fold(x)) it is possible to change to a DFT in the analysis filterbank, this corresponds to using a clockwise commutator [3] and type 3 polyphase filter structure [4].\n", "\n", + "Sampling and anlogue signal yields the low pass spectrum and its replicas around multiples of fs. This endles replication is why the discrete samples can represent and anlogue signal. Downsampling of the digitized signal again causes aliasing of all replicas to baseband. Normally the interest is in the baseband replica. However using a BPF instead of an LPF any of the replicas can be selected for baseband. Similar for upsampling (synthesis) the baseband signal can be upconverted to the wanted replica by using a BPF instead of an LPF. The first step in upsampling is inserting zeros to increase the sample rate, without changing the digital signal. However this increased sample rate does make the replicas available in the digital domain.\n", + "\n", + "With a single channel mixer the downsampling by M implies that at the downsampled rate there is processing power to process M - 1 more channels. This fits a critically sampled DFT filterbank with K = M. Whereby the DFT should typically be an FFT.\n", + "\n", + "The single channel mixer and DFT filterbank can share a prototype LPF for all channels, because the downsampling and upsampling performs the demodulation and modulation, to and from baseband at 0 Hz. When Ros = K / M > 1 then there remains and offset frequency that is compensated by a phasor term W_K^(knM). With the DFT filterbank this phasor that can be implemented as a circular shift (modulo K) at the input of the DFT. \n", + "\n", + "For perfect reconstruction by an analysis-synthesis filterbank the aliasing must cancel and distortion must be zero. Considering only adjancent channels is called pseudo QMF (two-channel), for the other channels the stop band attenuation suppresses the crosstalk. For zero distortion the transfer function of the channels must be power complementary. With an DFT filterbank (= complex modulated filterbank) the aliasing does not cancel. For critically sampled filterbanks this leads to the Modified DFT (MDFT) or Cosine modulated fiterbank [4]. By using oversampling the aliasing can be avoided. Hence using some oversampling and power complementary channel filters the oversampled DFT filterbank can achieve almost perfect reconstruction.\n", + "\n", + "Remarks:\n", + "* With oversampling the subband bandwidth can be increased to have a flat reponse up to fsub / 2. The transition band is then from fsub / 2 to Ros * fsub / 2 The subband filter is then not power complementary (around fsub / 2), so the subbands are not suitable for synthesis. It may be feasible achieve almost perfect reconstruction by first removing the transition band from the subbands, by separating the subbands into fine channels using a DFT, then zero the fine channels of the transistion band, and then synthesize the subband using IDFT.\n", + "\n", "Features:\n", "* Use full rate model of single channel down converter and up converter as expected exact golden reference result for the efficient polyphase implementation and WOLA implementation and refBunton.\n", "* Use SNR of input and difference with time aligned reconstructed output, to show that perfect reconstruction depends on the pass band gain being one and the stopband gain approaching zero. Hence the SNR for center bin sine wave inputs improves, towards perfect reconstruction, when Ntaps is increased.\n", "* Use analysis DFT filterbank and synthesis IDFT filterbank and verify that it yields the same result as with the single channel pipeline.\n", "* Compare model with PFB reconstruction by reconstruct.m of John Bunton for Ntaps = 12 and Ndft = 192, using refBunton.\n", - "* Support Ncoefs <= Ndelays = Ntaps * Ndft, to verify correct implementation of coefs and delay line\n", - "* Support WG with offset bin center frequency or with noise, to verify that result is correct for any frequency, not only the center frequency\n", + "* Support Ncoefs <= Ndelays = Ntaps * Ndft, to verify correct implementation of coefs and delay line.\n", + "* Support WG with offset bin center frequency (e.g. wgSub = 1.5) or with noise (via SNR_WG_dB < 100), to verify that result is correct for any frequency, not only the center frequency.\n", "* Support refBunton with asymmetrical hPrototype, to verify correct order implementation of coefs in analysis and synthesis. The refBunton has slightly asymmetrical hPrototype due to that it uses interpolation to increase the number of coefficient after designing the FIR filter with fircls1().\n", "\n", "References:\n", "1. dsp_study_erko, summary of DSP books\n", "2. chapter 7 in [CROCHIERE]\n", "3. chapter 6 downconverter, 7 upconverter, 9 filterbank in [HARRIS]\n", - "4. chapter 1.1.2, 8.6 in [FLIEGE]" + "4. chapter 1.1.2 (commutator), 7.3 (pseudo QMF), 7.3.3 (MDFT), 7.4 (Cos filterbank) 8.6 in [FLIEGE]" ] }, { @@ -152,8 +163,9 @@ " Ndft = 16 # DFT size\n", " #Ntaps = 4\n", " #Ndft = 8\n", - " Ntaps = 12\n", - " Ndft = 192\n", + " #Ntaps = 12\n", + " #Ndft = 192\n", + " Ntaps = 32\n", " # Generic hPrototype is symmetrical, so then coefs flip makes no difference\n", " flipAPrototype = False\n", "Ndelays = Ndft * Ntaps\n", @@ -189,6 +201,7 @@ "else:\n", " Ndown = Ndft * 3 // 4 # oversampled PFB\n", " #Ndown = Ndft * 7 // 8\n", + " #Ndown = Ndft * 15 // 16\n", " #Ndown = Ndft // 4\n", " #Ndown = Ndft # Critically sampled PFB\n", "Ros = Fraction(Ndft, Ndown)\n", @@ -331,9 +344,11 @@ " # HFbank firwin hanning filter, but for perfect reconstruction HFpowerbank needs to\n", " # be flat.\n", " # Optimum flat HFpowerbank\n", - " if Ntaps == 8 and Ndft == 16:\n", + " if Ndft == 16 and Ntaps == 8:\n", " hpFactor = 1.108\n", - " elif Ntaps == 12 and Ndft == 192:\n", + " elif Ndft == 16 and Ntaps == 32:\n", + " hpFactor = 1.027\n", + " elif Ndft == 192 and Ntaps == 12:\n", " hpFactor = 1.065\n", " else:\n", " # Default for optimum flat HFbank\n", @@ -529,7 +544,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/tmp/ipykernel_45544/4263550011.py:2: UserWarning: The filter's denominator is extremely small at frequencies [0.270, 0.295, 0.319, 0.344, 0.368, 0.393, 0.417, 0.442, 0.466, 0.491, 0.515, 0.540, 0.565, 0.589, 0.614, 0.638, 0.663, 0.687, 0.712, 0.736, 0.761, 0.785, 0.810, 0.834, 0.859, 0.884, 0.908, 0.933, 0.957, 0.982, 1.006, 1.031, 1.055, 1.080, 1.104, 1.129, 1.154, 1.178, 1.203, 1.227, 1.252, 1.276, 1.301, 1.325, 1.350, 1.374, 1.399, 1.424, 1.448, 1.473, 1.497, 1.522, 1.546, 1.571, 1.595, 1.620, 1.644, 1.669, 1.694, 1.718, 1.743, 1.767, 1.792, 1.816, 1.841, 1.865, 1.890, 1.914, 1.939, 1.963, 1.988, 2.013, 2.037, 2.062, 2.086, 2.111, 2.135, 2.160, 2.184, 2.209, 2.233, 2.258, 2.283, 2.307, 2.332, 2.356, 2.381, 2.405, 2.430, 2.454, 2.479, 2.503, 2.528, 2.553, 2.577, 2.602, 2.626, 2.651, 2.675, 2.700, 2.724, 2.749, 2.773, 2.798, 2.823, 2.847, 2.872, 2.896, 2.921, 2.945, 2.970, 2.994, 3.019, 3.043, 3.068, 3.093, 3.117], around which a singularity may be present\n", + "/tmp/ipykernel_27035/4263550011.py:2: UserWarning: The filter's denominator is extremely small at frequencies [0.270, 0.295, 0.319, 0.344, 0.368, 0.393, 0.417, 0.442, 0.466, 0.491, 0.515, 0.540, 0.565, 0.589, 0.614, 0.638, 0.663, 0.687, 0.712, 0.736, 0.761, 0.785, 0.810, 0.834, 0.859, 0.884, 0.908, 0.933, 0.957, 0.982, 1.006, 1.031, 1.055, 1.080, 1.104, 1.129, 1.154, 1.178, 1.203, 1.227, 1.252, 1.276, 1.301, 1.325, 1.350, 1.374, 1.399, 1.424, 1.448, 1.473, 1.497, 1.522, 1.546, 1.571, 1.595, 1.620, 1.644, 1.669, 1.694, 1.718, 1.743, 1.767, 1.792, 1.816, 1.841, 1.865, 1.890, 1.914, 1.939, 1.963, 1.988, 2.013, 2.037, 2.062, 2.086, 2.111, 2.135, 2.160, 2.184, 2.209, 2.233, 2.258, 2.283, 2.307, 2.332, 2.356, 2.381, 2.405, 2.430, 2.454, 2.479, 2.503, 2.528, 2.553, 2.577, 2.602, 2.626, 2.651, 2.675, 2.700, 2.724, 2.749, 2.773, 2.798, 2.823, 2.847, 2.872, 2.896, 2.921, 2.945, 2.970, 2.994, 3.019, 3.043, 3.068, 3.093, 3.117], around which a singularity may be present\n", " w, gd = signal.group_delay((hPrototype, [1.0]), w=1024, fs=1)\n" ] }, @@ -702,7 +717,7 @@ { "data": { "text/plain": [ - "<matplotlib.legend.Legend at 0x7f1042f35fd0>" + "<matplotlib.legend.Legend at 0x7fb4e035dee0>" ] }, "execution_count": 21,