Merge branch 'L2SDP-793' into 'master'

Clarified flipped order of FIR coefficients in PFB.

Closes L2SDP-793

See merge request desp/apertif_matlab!1

matlab/one_pfb.m

@@ -29,9 +29,13 @@
 %    . PFIR coeffcients:                      one_pfb_m_pfir_coeff_fircls1_16taps_128points_16b.dat
 %    . PFIR coeffcients, WG, PFIR, PFFT data: one_pfb_m_chirp_8b_16taps_128points_16b_16b.dat
 %
+% References:
+%
 % Description :
+% [1] Multirate Digital Signal Processing (Crochiere)
 %
 % * General
+%   Section 7 in [1] describes a poly phase filterbank (PFB).
 %   The data path (DP) is modelled per block of data. The block of data
 %   are counted by the block sequence number that thus acts as a
 %   timestamp. The data path consist of DSP functions. Each DSP
@@ -83,6 +87,68 @@
 %   filterbank is defined by the FFT size (tb.subband_fft_size). The block
 %   size after the subband filterbank is defined by the number of subbands
 %   (tb.nof_subbands).
+%
+% * Flipped order of FIR coefficients:
+%   The FIR filter in the PFB is defined by a prototype filter h[] with N_fft
+%   * N_taps = 1024 *16 = 16384 coefficients. The impulse response of this
+%   prototype filter will show the coefficients in h[0:16383] order. The
+%   FIR coefficients of h[] in the PFB are flipped per column of N_fft = 1024
+%   coefficients, for all N_taps = 16 columns. In VHDL simulation this shows
+%   with tb_verify_pfb_response.vhd and in the Matlab model this flip is also
+%   done by ctrl_pfir_subband.coeff = flipud(ctrl_pfir_subband.coeff).
+%   With the flip the one_pfb.m yields SNR subband / spurious max = 64.305310
+%   dB in subband 64 for an input WG at 64.015625. Without the flip the SNR
+%   subband / spurious max = 29.857938 dB, so much worse. Hence flipping the
+%   FIR coefficients is needed, and not due to an implementation detail in
+%   the VHDL. The run_pfb.m is equivalent to one_pfb.m and is used to create
+%   reference input and expected output data for the wpfb implementation in
+%   VHDL by wpfb_unit_wide.vhd. The tb_tb_wpfb_unit_wide.vhd verifies multiple
+%   sets of reference data including sinus, chirp, noise and real input or
+%   complex input and wideband (f_sample > f_clk) or same rate (f_sample =
+%   f_clk). This tb guarantees that the VHDL agrees with the Matlab model.
+%   In the Matlab model the data blocks and data time are defined as:
+%
+%    WG            FIR with four taps                   FFT
+%    index     t                                          t
+%        0 -1023   1023   2047  3071  4095  --> + --> -1023
+%              .      .      .     .     .
+%              .      .      .     .     .
+%             -2      2      .     .     .  --> + -->    -2
+%             -1      1      .     .     .  --> + -->    -1
+%     1023     0      0   1024  2048  3072  --> + -->     0
+%              d      h      h     h     h                d
+%
+%   The waveform generator (WG) generates a block of 1024 samples with index
+%   0:1023, where sample at index 1023 is the newest and sample at index 0 is
+%   the oldest. In a filter the newest sample needs to be multiplied with h[0]
+%   and the older samples are multiplied by the subsequent coefficients,
+%   because it is a convolution. Therefore the FIR coefficients need to be
+%   flipped up/down per column, to allow doing the filter as a d .* h vector
+%   multiply in pfir.m. The FIR filter sums the rows for the N_taps = 16 to
+%   yield the input for the FFT. The FFT operates on blocks of data with same
+%   index and time range as the WG. The data output of the FIR filter fits
+%   this input range of the FFT. Therefore no data flipping is needed.
+%
+%   In summary:
+%   * The WG data and FFT input data is processed in blocks with the newest
+%     sample at the bottom and the oldest sample at the top. In the FIR filter
+%     the newest block is at the left and the oldest at the rigth. Therefor
+%     the FIR filter coefficients also have to be ordered from bottom left to
+%     top right.
+%   * The column size is determined by N_fft of the FFT. The FFT is calculated
+%     each time a block is shifted in. For a critically sampled PFB the block
+%     size is N_blk = N_fft, so then it also looks like the blocks only shift
+%     from left to right. Using N_blk < N_fft would yield an oversampled PFB
+%     with oversampling factor R_os = N_fft/N_blk. Then the blocks shift from
+%     top to bottom in each column and from left to right.
+%
+%   In the tb_verify_pfb_response.vhd the input stimuli is a block of N_fft =
+%   1024 ones followed by N_taps-1 blocks with zeros. In time the oldest data
+%   will appear first in the simulator Wave Window, so therefore the
+%   fil_re_scope signal in the tb_verify_pfb_response.vhd will show the FIR
+%   coefficients h[0:16383] in order h[1023:0], h[2047:1024], ...,
+%   h[16383:15360], so flipped per block.
+%
 clear all;
 close all;
 fig=0;
@@ -205,7 +271,7 @@ if strcmp(tb.model_filterbank, 'LOFAR')
     ctrl_pfir_subband.nof_taps           = 16;                                                            % Number of taps
     ctrl_pfir_subband.nof_coefficients   = ctrl_pfir_subband.nof_polyphases*ctrl_pfir_subband.nof_taps;   % Number of filter coefficients (taps)
     ctrl_pfir_subband.data_w             = 16;
-    ctrl_pfir_subband.config.design      = 'lofar file';
+    ctrl_pfir_subband.config.design      = 'lofar_file';
     hfir_subband_coeff = load('data/Coeffs16384Kaiser-quant.dat');
     hfir_subband_coeff = hfir_subband_coeff/max(hfir_subband_coeff);
     hfir_subband_coeff = hfir_subband_coeff';           % Use column vector, same format as by pfir_coeff()

matlab/run_pfb.m

@@ -236,6 +236,7 @@ x = fliplr(x);
 
 % Flip ctrl_pfir_subband.coeff per poly phase, because that is the order
 % in which the pfir() model and HDL expect the FIR coefficients
+% See one_pfb.m for more detailed clarification.
 ctrl_pfir_subband.coeff = reshape(hfir_subband_coeff, ctrl_pfir_subband.nof_polyphases, ctrl_pfir_subband.nof_taps);
 ctrl_pfir_subband.coeff = flipud(ctrl_pfir_subband.coeff);
 ctrl_pfir_subband.Zdelays = zeros(ctrl_pfir_subband.nof_polyphases, ctrl_pfir_subband.nof_taps-1);

matlab/run_pfb_complex.m

@@ -203,7 +203,7 @@ hfir_channel_coeff = pfir_coeff(ctrl_pfir_channel.nof_polyphases, ...
                                 ctrl_pfir_channel.coeff_w, ...
                                 ctrl_pfir_channel.config);
 ctrl_pfir_channel.coeff = reshape(hfir_channel_coeff, ctrl_pfir_channel.nof_polyphases, ctrl_pfir_channel.nof_taps);
-ctrl_pfir_channel.coeff = flipud(ctrl_pfir_channel.coeff);
+ctrl_pfir_channel.coeff = flipud(ctrl_pfir_channel.coeff);  % See one_pfb.m for more detailed clarification.
 ctrl_pfir_channel.Zdelays = zeros(ctrl_pfir_channel.nof_polyphases, ctrl_pfir_channel.nof_taps-1);
 ctrl_pfir_channel.gain = 1;   % no gain adjustment in PFIR, just apply the coefficients to the input data