diff --git a/matlab/one_pfb.m b/matlab/one_pfb.m
index f8434022c4a271275bc509248e0006918f2ee3fb..b24198eceb0e5d3b5df5413e7b1453bdadf411b2 100644
--- a/matlab/one_pfb.m
+++ b/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()
diff --git a/matlab/run_pfb.m b/matlab/run_pfb.m
index 55455b2a997a28eeb1abccb1d175dca9f6eee1d0..30025d3cbd5dfbfb531b9f8b921d718b8c2ac506 100644
--- a/matlab/run_pfb.m
+++ b/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);
diff --git a/matlab/run_pfb_complex.m b/matlab/run_pfb_complex.m
index a9baa0a4a3508e13ea1d0d7bf1bf6f946947aec1..7716e4a640a0d884eed44e613b3d1a6a174ebbf0 100644
--- a/matlab/run_pfb_complex.m
+++ b/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