Skip to content
Snippets Groups Projects
Commit a5e88f2e authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Add John Buntons code, but with annotations by Eric Kooistra. Same results.

parent 233d29af
No related branches found
No related tags found
1 merge request!419Resolve RTSD-265
Hide whitespace changes
Inline Side-by-side
Showing
with 355 additions and 0 deletions
applications/lofar2/model/pfb_bunton_annotated/Gen_filter12.m 0 → 100644
+ 100
− 0
View file @ a5e88f2e
% 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
applications/lofar2/model/pfb_bunton_annotated/README.txt 0 → 100644
+ 54
− 0
View file @ a5e88f2e
>>> 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
applications/lofar2/model/pfb_bunton_annotated/polyphase_analysis.m 0 → 100644
+ 38
− 0
View file @ a5e88f2e
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))
applications/lofar2/model/pfb_bunton_annotated/polyphase_synthesis_b.m 0 → 100644
+ 41
− 0
View file @ a5e88f2e
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
applications/lofar2/model/pfb_bunton_annotated/reconstruct.m 0 → 100644
+ 122
− 0
View file @ a5e88f2e
%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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment