Resolve RTSD-265

Closes RTSD-265

applications/lofar2/model/pfb_bunton_annotated/Gen_filter12.m

+% 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