diff --git a/applications/lofar2/model/pfb_bunton_annotated/Gen_filter12.m b/applications/lofar2/model/pfb_bunton_annotated/Gen_filter12.m
new file mode 100644
index 0000000000000000000000000000000000000000..763902cc8e49faf79fe9cbca1637a228e09b9f42
--- /dev/null
+++ b/applications/lofar2/model/pfb_bunton_annotated/Gen_filter12.m
@@ -0,0 +1,100 @@
+% Code : John Bunton, CSIRO
+% Comment : Eric Kooistra, ASTRON
+
+close all
+
+% B = fircls1(N, WO, DEVP, DEVS)
+% LPF FIR filter design by constrained least-squares
+% . N is filter order, so length(B) = N + 1 is number of coefficients
+% . W0 is LPF cutoff frequency between 0 and 1 = fNyquist = fs/2
+% . DEVP maximum linear deviation (ripple) in pass band
+% . DEVS maximum linear deviation (ripple) in stop band
+% Parameter values:
+% . length(c) = 192 coefs
+% . W0 ~= fNyquist / 16. For 1/Kfft = 1/16 power at 0.5 bin is -5.941 dB, so
+%   gain is 10^(-5.941 / 20) = 0.5046. For 1/15.3 gain is 10^(-3.362 / 20) =
+%   0.6790. Probably adjusted W0 to have gain = sqrt(0.5) =  0.7071, so that
+%   analysis LPF and synthesis LPF together have gain 0.5 at 0.5 bin. Use
+%   1/15.19 to have gain is 10^(-3.012 / 20) = 0.7070. However from the
+%   magnitude response of bin 0 and 1 it shows that 1/15.19 does yield gain
+%   1 at 0.5 bin, but with a raise in average ripple level at the bin edges.
+%   With 1/15.3 the average ripple level is almost constant 1 from one bin
+%   to the next.
+% . The DC gain is slightly > 1.0. Do not normalize the DC gain by sum(B),
+%   because fircls1() yields such a DC gain that the LPF pass band power
+%   level is equivalent to that of an ideal (brick wall) LPF with W0.
+Kfft = 16;  % FFT block size
+Ntaps = 12;  % number of FFT blocks
+Kcoefs = Kfft * Ntaps;
+Norder = Kcoefs - 1;
+W0 = 1 / 15.3;
+%W0 = 1 / 15.19;
+c = fircls1(Norder, W0, .02, .0003);
+
+figure(1)
+subplot(3, 1, 1)
+k = Kcoefs - 1;
+plot((0:k) / k * Ntaps, c)
+xlabel('Time in units of FFT length','FontSize',10)
+title('Filter impulse response')
+grid
+
+% Use DTFT as fft with zero padding for factor R = 4 finer frequency
+% resolution. Number of frequency point per bin is Kcoefs / Kfft =
+% Ntaps, so with R = 4 this yields 12 * 4 = 48 frequency points per bin.
+R = 4;
+dtft_c = fft([c c*0 c*0 c*0]);
+
+subplot(3, 1, 2)
+% db(volt) = 10 * log10(volt^2) = 20 * log10(volt)
+%plot((0:48)/48, db(dtft_c(1:49)))
+plot(((0:96)-48)/48, [db(dtft_c(49:-1:2)), db(dtft_c(1:49))])
+hold on
+%plot((48:-1:0)/48, db(dtft_c(1:49)), '.')
+plot(((0:96)-48)/48, db(dtft_c(97:-1:1)), '.')
+xlabel('Frequency (FFT bins)', 'FontSize', 10)
+ylabel('Magnitude (dB)', 'FontSize', 10)
+title('Response of Bin 0 and 1 to monochromatic source')
+hold off
+axis([-1 1 -80 10])
+grid
+
+subplot(3, 1, 3)
+% Plot magnitude squared response, to have correlator (power) response.
+% Alternatively this also accounts for combined magntiude response of analysis
+% LPF and synthesis LPF in case of reconstruction.
+plot((0:48)/48, abs(dtft_c(1:49).^2) + abs(dtft_c(49:-1:1).^2))
+grid
+%axis([0 1 0.98 1.04])
+xlabel('Frequency (FFT bins)','FontSize',10)
+ylabel('Magnitude','FontSize',10)
+title ('Sum of correlator responses for Bin 0 and 1')
+
+% The FIR filter design does not converge for large Kcoefs, therefore use
+% interpolation to create longer FIR filter cl, with Ncoefs = Kcoefs * Q =
+% 2304, using Q = 12. The number of FFT blocks remains Ntaps = 12, because
+% the shape of cl is the same as for c. The FFT block size now becomes Nfft =
+% Kfft * Q = 16 * 12 = 192.
+Q = 12;
+% Y = interpft(X, N) returns a vector Y of length N obtained by interpolation
+% in the Fourier transform of X.
+cl = interpft(c, length(c) * Q) / Q;
+dc_gain_c = sum(c);
+dc_gain_cl = sum(cl);
+Ncoefs = Kcoefs * Q;
+Nfft = Kfft * Q;
+
+figure(2)
+fft_cl = fft(cl);
+% Plot magnitude frequency response from DC = 0 to m bins
+m = 2;
+Nppb = Ntaps * R;  % number of frequency points per bin
+plot((0:m*Nppb)/Nppb, db(dtft_c(1:m*Nppb + 1)), 'b')
+hold on
+Nppb = Ntaps;  % number of frequency points per bin
+plot((0:m*Nppb)/Nppb, db(fft_cl(1:m*Nppb + 1)), 'r')
+legend('c', 'cl')
+xlabel('Frequency (FFT bins)','FontSize',10)
+ylabel('Magnitude','FontSize',10)
+title ('Response of fircls1 filter c and interpolated cl')
+grid
diff --git a/applications/lofar2/model/pfb_bunton_annotated/README.txt b/applications/lofar2/model/pfb_bunton_annotated/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f7029fb2ff965e811ea43ddad5b59370be4c3201
--- /dev/null
+++ b/applications/lofar2/model/pfb_bunton_annotated/README.txt
@@ -0,0 +1,54 @@
+>>> On 17/11/2015 at 02:47, in message
+<40E3C1B05457DB43B21A4F4C5F191A392F8E71EA@exmbx03-cdc.nexus.csiro.au>,
+<John.Bunton@csiro.au> wrote:
+> Hi Marco and Wallace,
+>
+> I know you see the reconstruction of the polyphase filterbank data to a
+> higher sample rate at somewhat unknown and risky.
+> Attached are Matlab files to demonstrate a reconstruction.
+>
+> To use them
+> First generate a prototype filter by executing Gen_filter12.m
+>                 This will also plot up the impulse response of the filter,
+> the response of two adjacent channel  and the response to a monochromatic
+> signal.
+>
+> Next run reconstruct.m  with the following settings
+> ImpulseResponse = 1;
+> Correct = 1;
+>  this will calculate an impulse response and a correction filter for the
+> transfer function error
+>
+> Finally run reconstruct.m with
+> ImpulseResponse = 0;
+> Correct = 1;
+>                 With the settings in the attached files
+> Figure 1 shows the input (blue) and reconstructed signal (red), first 10,000
+> points
+> Figure 2 shows the difference between the reconstructed and input
+> Figure 2 show the spectrum on the input and the difference.
+>
+> Input is a sinewave, line 30 if you want to play with the frequency
+> Comment out lines 29 and 30 to generate a random signal.
+>
+> Things to change are the prototype filter line 1 of Gen_filter12 and
+> Oversampling ratio line 40 of reconstruct.m , note variable k at line 62
+> must be set correctly.
+>
+> There is a trade-off to be made.  With a sufficiently long FIR on the
+> analysis and synthesis filterbank or sufficiently high oversampling ration we
+> do not need the transfer function correction filter.   However there are 1024
+> filterbanks compared to the required 32 correction filters.  So I expect to
+> have a correction filter and for the filterbanks to be similar to what is in
+> the attached files.
+>
+> Regards,
+> John
+>
+> John Bunton PhD BE BSc
+> Project Engineer, SKA LOW Correlator and Beamformer
+> Astronomy and Space Science
+> CSIRO
+> E John.bunton@csiro.au T +61 1 9372 4420 M 04 2907 8320
+> PO Box 76 Epping, 1710, NSW, Australia
+> www.csiro.au | www.atnf.csiro.au/projects/askap
diff --git a/applications/lofar2/model/pfb_bunton_annotated/polyphase_analysis.m b/applications/lofar2/model/pfb_bunton_annotated/polyphase_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..cb6fc0152d9abdb000c1c603111d5493b661d118
--- /dev/null
+++ b/applications/lofar2/model/pfb_bunton_annotated/polyphase_analysis.m
@@ -0,0 +1,38 @@
+function outBins = polyphase_analysis(inData, filt, Nfft, Nstep)
+
+% John Bunton CSIRO 2003
+
+% inData = input time data
+% filt = coefficients of prototype lowpass FIR filter
+% Nfft = length of FFT, prefilter length = length(filt) / Nfft, if not then
+%        'filt' is padded with zeros to a multiple of Nfft
+% Nstep = increment along inData between FFT output blocks, so downsample
+%         inData by Nstep, use Nstep = Nfft for critical sampling of outBins
+% outBins = output data 2D one dimension time the other frequency
+
+% Support padding Ncoefs with zeros, if filt is too short
+Ntaps = ceil(length(filt) / Nfft);
+Ncoefs = Ntaps * Nfft;
+bfir = (1:Ncoefs) * 0;  % = zeros(1, Ncoefs)
+bfir(1:length(filt)) = filt;
+
+nof_blocks = floor((length(inData) - Ncoefs) / Nstep);
+outBins = zeros(nof_blocks, Nfft);
+
+% Remove comment to try interleaved filterbank, produces critically sampled
+% outputs as well as intermediate frequency outputs
+%Nfft = Nfft * 2;
+%Ntaps = Ntaps / 2;
+
+for k = 0:nof_blocks-1
+    % Filter input using bfir as window, with downsampling rate Nstep
+    temp = bfir .* inData((1:Ncoefs) + Nstep*k);
+    % Sum the FIR filter taps
+    temp2 = (1:Nfft) * 0;
+    for m = 0:Ntaps-1
+        temp2 = temp2 + temp((1:Nfft) + Nfft*m);
+    end
+    % FFT
+    outBins(k+1, 1:Nfft) = fft(temp2);
+end
+%plot(real(temp2))
diff --git a/applications/lofar2/model/pfb_bunton_annotated/polyphase_synthesis_b.m b/applications/lofar2/model/pfb_bunton_annotated/polyphase_synthesis_b.m
new file mode 100644
index 0000000000000000000000000000000000000000..3e1ae4f3660abc79e726da267af901ef472228a0
--- /dev/null
+++ b/applications/lofar2/model/pfb_bunton_annotated/polyphase_synthesis_b.m
@@ -0,0 +1,41 @@
+function outData = polyphase_synthesis_b(inBins, filt, Nstep)
+
+% polyphase synthesis on blocks of filterbank data
+%       John Bunton original version Oct 2003
+
+% inBins = input bin data, 2D array of multiple blocks of FFT data.
+%          Nfft = number of frequency bins
+% filt = prototype lowpass filter
+% Nstep = increment along data at output of IFFT, so upsample inBins by
+%         Nstep, use Nstep = Nftt for critically sampled inBins
+% outData = reconstructed output time data
+%
+s = size(inBins);
+nof_blocks = s(1);  % number of FFT blocks
+Nfft = s(2);  % length of FFT
+
+% Support padding Ncoefs with zeros, if filt is too short
+Ntaps = ceil(length(filt) / Nfft);
+Ncoefs = Ntaps * Nfft;
+bfir = (1:Ncoefs) * 0;
+bfir(1:length(filt)) = filt;
+
+% IFFT
+for n = 1:nof_blocks
+    temp(n, 1:Nfft) = real(ifft(inBins(n,:)));
+end
+
+outData = (1:nof_blocks*Nstep + Ncoefs) * 0;
+for n = 1:nof_blocks
+    % Copy the data at each tap
+    temp2 = (1:Ncoefs) * 0;
+    for m = 0:Ntaps-1
+        temp2((1:Nfft) + m*Nfft) = temp(n, 1:Nfft);
+    end
+    % Filter input using bfir as window
+    temp2 = temp2 .* bfir;
+    % Overlap add to sum the taps, with upsampling rate Nstep
+    range = (1:Ncoefs) + (n-1)*Nstep;
+    outData(range) = outData(range) + temp2;  % for continuous time output
+    %outData(n,:) = temp2;  % for output in blocks
+end
diff --git a/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m b/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m
new file mode 100644
index 0000000000000000000000000000000000000000..22c31fcdca64fc4d7a37127c3c4f9bcc05663c52
--- /dev/null
+++ b/applications/lofar2/model/pfb_bunton_annotated/reconstruct.m
@@ -0,0 +1,122 @@
+%reconstruction of signal after polyphase filterbank
+% John Bunton created Dec 2005
+%
+% modified to add variables ImpulseResponse and Correct
+% Impulse Response = 1 calculates the impulse response of the system
+% the result is in variable outData.  It also causes the transfer function
+% correction filter to be calculated, Result in variable correct
+% Correct = 1 applies the correction filter to the reconstructed signal
+% outData
+% John Bunton 2015
+
+close all
+
+ImpulseResponse = 0;  % set to 1 to calculate correction filter
+Correct = 1;          % set to 1 to allow correction filter to be applied
+
+len = 5000;
+z = rand(len,1);
+z = z > 0.5;
+z = 2 * (z - 0.5);
+
+% cl is the prototype lowpass filter
+Ncoefs = length(cl);  % = 2304
+Ntaps = 12;
+Nfft = Ncoefs / Ntaps;  % = 192
+
+% Create random block wave signal with block size Nhold = 16 samples, first
+% insert 16-1 zeros between every sample in z. This yield length
+% (len - 1) * 16 + 1. Then convolve with ones(16, 1) to replicate each value
+% of z instead of the zeros. The length of conv(A, B) is length(A) +
+% length(B) - 1, so length of inData becomes ((len - 1) * 16 + 1) + (16 - 1) =
+% len * 16 = 80000.
+clear inData
+Nhold = 16;
+inData(1:Nhold:len*Nhold) = z;
+inData = conv(inData, ones(Nhold, 1));
+k = len * Nhold;  % = length(inData)
+
+%inData = cos(((1:k).^2 ) * len / (k)^2);  % comment out these two line to have random signal input
+%inData = cos(2*pi*(1:k) * 0.031 / 16);
+
+if ImpulseResponse == 1
+    inData = inData .* 0;   % for testing impulse response of system
+    inData(k / 2) = 1;  % and allows calculation of correction filter
+end
+
+% Set up parameters for the analysis and synthesis filterbank
+Ros = 32 / 24;
+% fix() rounds to nearest integer towards zero.
+Nstep = fix(Nfft / Ros);  % = 144
+clf
+
+% these is the actual analysis and synthesis operations
+outBins = polyphase_analysis(inData, cl, Nfft, Nstep);
+outData = polyphase_synthesis_b(outBins, cl, Nstep);
+
+% adjust gain
+% . unit gain in prototype LPF Gen_filter12.m
+adjust = Nfft;  % normalize gain per bin for analysis LPF
+adjust = adjust / Ros;  % compensate for oversampling factor
+adjust = adjust * Nfft;  % normalize gain per bin for synthesis LPF
+outData = adjust * outData;
+
+if ImpulseResponse == 1
+     % correction filter to fix frequency response error, also
+     % corrects gain because input inData has unit level
+    correct = real(ifft(1./fft(outData)));
+end
+
+% Choose parameter k to suit oversampling ratio
+% oversampling ratio 32/23 then k=0
+% oversampling ratio 32/24 then k=39
+% oversampling ratio 32/25 then k=72
+% oversampling ratio 32/26 then k=3
+% oversampling ratio 32/27 then k=48
+k=39;
+corr = correct;
+% Correct for time shift for transfer function correction filter
+corr(1+k:length(corr)) = corr(1:length(corr)-k);
+
+% Correct for errror in transfer function
+if Correct == 1
+    outData = ifft(fft(outData) .* fft(fftshift(corr)));
+end
+
+% Take the difference to see the error
+diff = outData - inData(1:length(outData));
+
+% SNR
+% First do rng('default') to reset the random number generator, then
+% run reconstruct. The original code by John Bunton then yields:
+% . rng('default'); reconstruct, with Correct = 0: SNR = 30.56 dB
+% . rng('default'); reconstruct, with Correct = 1: SNR = 62.24 dB
+SNR = 10 * log10(sum(abs(inData(10000:70000) .^ 2)) / sum(abs(diff(10000:70000) .^ 2)));
+sprintf('SNR = %6.2f [dB]', SNR)
+
+%plot(diff)
+%plot(db(fft(outData))) %to see frequency error
+%plot(db(fft(diff(length(inData)/8:7*length(inData)/8))))
+
+figure(1)
+hold off
+plot(db(fft(inData(10000:70000)/30000)),'r')
+hold on
+plot(db(fft(diff(10000:70000)/30000)),'b')
+grid
+title('Spectrum of input and error')
+
+figure(2)
+plot(diff)
+title('Difference between input and reconstructed')
+
+figure(3)
+hold off
+plot(inData(1:10000))
+hold on
+title('Input and Reconstructed signal')
+if (ImpulseResponse == 1)
+    plot(outData)
+else
+    plot(outData(1:10000),'r')
+end