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

Updated model/pfs_coeff_final.m, but preserved old FTF results for reference... 

Updated model/pfs_coeff_final.m, but preserved old FTF results for reference via reproduce_ftf = true.
parent a3ad554c
Branches
No related tags found
1 merge request!188Updated model/pfs_coeff_final.m, but preserved old FTF results for reference...
Pipeline #23419 passed
Stage: simulation
Stage: synthesis
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/pfs_coeff_final.m
+ 389
− 107
View file @ 3a0f7298
%------------------------------------------------------------------------------
%
% Copyright 2021
% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
%------------------------------------------------------------------------------
%
% Authors :
% Created by FilterTaskForce team 2004 for LOFAR1 [1, 2]
% Updated in 2022 for LOFAR20 by E. Kooistra
% Purpose:
%
% References:
% [1] hdl/applications/lofar1/FilterTaskForce/readme_FilterTaskForce.txt
% [2] hdl/applications/lofar1/FilterTaskForce/pfs_coeff_final.m (Original
% file from LOFAR1)
% [3] LOFAR_MEM_162_digfilter_stopband3.pdf
% [4] data/Coeffs16384Kaiser-quant.dat original LOFAR1 FIR coefficients
% [5] data/Coefficient_16KKaiser.gold FTF reference FIR coefficients
% [6] run_pfir_coeff.m, similar to this FIR filter generation script, used
% for Apertif.
%
% Description:
%
% The LOFAR1 coefficients [4] were probably generated with this script, but
% the exact same values cannot be reproduced, probably due to differences
% in parameter settings. Nevertheless the FIR filter design method in this
% scriptis still valid, so if we should need to generate new coefficients,
% then we can do so with this script.
%
% Using reproduce_ftf = true does reproduce the FTF coefficients from [5].
% However, these contain an error in r_pass, which could also be present the
% LOFAR1 coefficients [4], which means that the r_pass ripple could be
% improved.
%
% The half power subband bandwidth for the LOFAR1 [4] FIR filter is HP_BW =
% 0.9954 * f_sub and for the FTF [5] FIR filter HP_BW = 0.9947 * f_sub.
% Both FIR filters have the same maximum coefficient value of 2^15 - 1 =
% 32767.
%
% With a plain FFT one expects a processing gain of N_point for the signal
% to noise ration (SNR) of a bin, because the wanted signal is in that bin,
% while the noise is still in all N_point bins. However, in poly phase
% filterbank (PFB) simulations I (Eric) noticed that the PFB does not achieve
% a processing gain of N_point in SNR. I think this is due to that the poly
% phase FIR filter part in the PFB has a DC response that is not equal for
% all N_point polyphases, that form the N_point input to the FFT part of the
% PFB. There is a small variation in DC reponse that then causes leakage into
% other bins and thus reduces the processing gain of the PFB to less than
% N_points. This script shows that a more flat DC response over the poly
% phases can be achieved by using a stronger window function for the FIR
% coefficients calculation. However, more windowing would also increases the
% half power bandwidth of the subbands. For LOFAR the reduced PFB processing
% gain is no problem, because the PFB still meets the subband stop band
% attenuation requirement of [3].
%
% The LOFAR1 FIR coefficients [4] have nof_bits = 16 bits. This is enough to
% meet the specification of [3]. Therefore it is not necessary to regenerate
% the FIR coefficients to 18 bit, even though the multipliers in the FPGA
% could support 18b nearly for free.
%
% For comparision between LOFAR1 and LOFAR20 stations it is preferred to
% keep using the same subband FIR coefficients. This is fine, because they
% meet the specification of [3].
%
% Conclusion:
% a) We can and should keep the original LOFAR1 subband FIR filter coefficients
% also for LOFAR20.
% b) We cannot exactly reproduce the LOFAR1 subband FIR filter coefficients, but
% with this script we could generate an equivalent or improved set of subband
% FIR coefficients if we would have to.
%
% Modifications:
% a) Corrected missing -1 in r_pass=10^(r_pass_dB/20) - 1
% b) Change r_stop_dB = -89 dB into -91 dB conform [3]
% c) Clarified plots and corrected the f1_axis center
% d) Support round() instead of uencode() for coefficient quantization
% e) Compare created coefficients with LOFAR1 coefficients
% f) Added plot of DC response per poly phase
% g) Determine half power bandwidth of the pass band
% h) Support loading and analysing LOFAR1 coefficients for comparison
% i) Support changing the half power subband band width via relative_bw = 1
close all;
clear all;
fig = 0;
reproduce_ftf = true;
%reproduce_ftf = false;
%% User settings
if not(reproduce_ftf)
% passband ripple (in dB);
r_pass_dB = 0.5;
r_pass = 10^(r_pass_dB/20) - 1;
% stopband ripple (in dB);
r_stop_dB = -89;
r_stop = 10^(r_stop_dB/20);
% relative subband band width
relative_bw = 1;
% word size for coefficients
nof_bits = 16;
% window function for coefficients
kaiser_beta = 1;
use_window = 'kaiser';
% quantization method
use_uencode = false;
% support analysing LOFAR1 coefficients
use_lofar1 = false;
else
disp(sprintf('NOTE: recreate FTF coefficients result for reference'));
r_pass = 10^(0.5/20); % note missing -1
r_stop = 10^(-89/20); % -89 instead of -91 in [3]
relative_bw = 1;
nof_bits = 16;
kaiser_beta = 1;
use_window = 'kaiser';
use_uencode = true;
use_lofar1 = false;
end
fs =200*10^6;
coeff_max = 2^(nof_bits - 1) - 1;
%% Filter specification
% fft_size;
N = 1024;
% taps per frequency channel
L = 16;
% passband ripple (in dB);
r_pass=10^(0.5/20);
% stopband ripple (in dB);
r_stop=-89;
% upsampling factor
Q = 16;
% word size for coefficients
W=16;
% pass bandwidth
subband_bw = relative_bw / N;
comp_bw = Q * subband_bw;
% computed filter length
M1 = N*L/Q;
% interpolated filter length
M2 = N*L;
Nfft = M2;
% compute filter
c=fircls1(M1-1,Q/N,10^(r_pass/20),10^(r_stop/20),'trace');
% fine frequency resolution
P = 10;
M3 = M2 * P;
% frequency scale
fs1 = N/Q;
fs2 = N;
f1_axis = [-M1/2:M1/2-1]/M1 * fs1; % for compute filter
f2_axis = [-M2/2:M2/2-1]/M2 * fs2; % for interpolated filter
f3_axis = [-M3/2:M3/2-1]/M3 * fs2; % for fine resolution spectrum
figure(30)
Nfft =length(c);
hold on
spectrum_Fircls1= abs(fftshift(fft(c/max(c),Nfft)))/max(abs(fftshift(fft(c/max(c),Nfft))));
plot([-Nfft/2+1:Nfft/2]/Nfft*fs,20*log10(spectrum_Fircls1),'r')
hold off
legend('Other candidate','filter given by fircls1')
xlabel('f(Hz)')
ylabel('Gain(dB)')
title('Compirason of the prototype filter window before FFT interpolation');
%% Plot possible window functions
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot([kaiser(M1, 0.5), kaiser(M1, 1), kaiser(M1, 2), kaiser(M1, 4), hanning(M1), blackman(M1)])
legend('kaiser beta = 0.5', 'kaiser beta = 1', 'kaiser beta = 2', 'kaiser beta = 4', 'hanning', 'blackman')
title('Window functions');
xlabel('Sample')
ylabel('Gain')
grid on;
%% Compute filter
h_comp = fircls1(M1-1, comp_bw, r_pass, r_stop, 'trace');
fft_hcomp = fft(h_comp);
fftshift_comp = fftshift(fft_hcomp);
spectrum_comp = abs(fftshift_comp) / max(abs(fftshift_comp));
spectrum_comp_dB = 20*log10(spectrum_comp);
phase_comp = unwrap(angle(fftshift_comp));
%% Compute windowed filter
if strcmp(use_window, 'blackman')
disp(sprintf('NOTE: use blackman window'));
h_window = blackman(M1);
elseif strcmp(use_window, 'hanning')
disp(sprintf('NOTE: use hanning window'));
h_window = hanning(M1);
elseif strcmp(use_window, 'kaiser')
disp(sprintf('NOTE: use kaiser window (beta = %4.2f)', kaiser_beta));
h_window = kaiser(M1, kaiser_beta);
else % default 'none'
disp(sprintf('NOTE: use no window'));
h_window = ones(M1, 1); % no window is rectangular window
end
h_wcomp = h_window' .* h_comp;
%% Plot frequency response of compute filter
fft_wcomp = fft(h_wcomp);
fftshift_wcomp = fftshift(fft_wcomp);
spectrum_wcomp = abs(fftshift_wcomp) / max(abs(fftshift_wcomp));
spectrum_wcomp_dB = 20*log10(spectrum_wcomp);
phase_wcomp = unwrap(angle(fftshift_wcomp));
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(f1_axis, spectrum_comp_dB, 'r', f1_axis, spectrum_wcomp_dB, 'g')
legend('computed','windowed')
title('Power spectrum of computed and windowed filter');
xlabel('f [Hz]')
ylabel('Gain [dB]')
grid on;
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(f1_axis, phase_comp, 'r', f1_axis, phase_wcomp, 'g')
legend('computed','windowed')
title('Phase spectrum of computed and windowed filter'); % no difference
xlabel('f [Hz]')
ylabel('Phase [rad]')
grid on;
% use fourier interpolation to create final filter
f1=fft(kaiser(length(c),1).'.*c);%hanning(length(c)).'.*
%% Use fourier interpolation to create final filter
f2 = zeros(1, M2);
% copy the lower frequency half.
% copy the lower frequency half of the computed filter into f2.
n = 0:M1/2;
f2(1+n) = fft_wcomp(1+n);
figure(31)
subplot(2,1,1)
plot(-length(f1)/2+1:length(f1)/2,20*log10(fftshift(abs(f1))))
xlabel('f(Hz)')
title('Gain of the filter after FFT')
ylabel('Gain(dB)')
subplot(2,1,2)
plot(-length(f1)/2+1:length(f1)/2,unwrap(angle(fftshift(f1))))
title('Phase of the filter after FFT')
xlabel('f(Hz)')
ylabel('Gain(dB)')
f2(1+n)=f1(1+n);
% to make the impulse response exactly symmetric, we need to do a phase
% correction
% to make the impulse response exactly symmetric, we need to do a phase correction
% apply phase correction (= time shift) to make the impulse response exactly symmetric
f2(1+n) = f2(1+n) .* exp(-sqrt(-1)*(Q-1)*pi*n*1/M2);
% create the upper part of the spectrum from the lower part
f2(M2-n) = conj(f2(2+n));
figure(32)
subplot(2,1,1)
plot(-length(f2)/2+1:length(f2)/2,20*log10(fftshift(abs(f2))))
xlabel('f(Hz)')
title('Gain of the filter after FFT')
ylabel('Gain(dB)')
subplot(2,1,2)
hold on
plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))))
plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))),'r')
hold off
title('Phase of the filter after FFT')
xlabel('f(Hz)')
ylabel('Gain(dB)')
% back to time domain
c2=real(ifft(f2));
h_fourier = real(ifft(f2));
% use same scale
h_fourier = h_fourier * coeff_max / max(h_fourier);
% quantize the coefficients
c3=c2/max(c2);
coeffs_quant = double(uencode(c3,16,max(c3),'signed'));
% q=1/(2^(W-1));
% c3=q*round(c3/q);
% c3=c3*max(c2)/max(c3);
if use_uencode
% uencode is equivalent to floor()
disp(sprintf('NOTE: use uencode() quantization'));
h_quant = double(uencode(h_fourier, nof_bits, coeff_max, 'signed'));
else
% round is equivalent to floor(x + 0.5) provided that exactly 0.5 does not
% because round(0.5) = floor(x + 0.5) = 1, but round(-0.5) = -1.
disp(sprintf('NOTE: use round() quantization'));
h_quant = round(h_fourier);
end
%% Load LOFAR1 coefficients
if use_lofar1
h_quant = load('data/Coeffs16384Kaiser-quant.dat');
h_quant = h_quant';
end
%% Check symmetry of FIR coefficients
n=0:M2-1;
h_reverse = h_quant(M2-n);
d = h_quant == h_reverse;
if sum(d) == length(d)
disp(sprintf('NOTE: quantized FIR coefficients are symmetrical.'));
else
disp(sprintf('WARNING: quantized FIR coefficients are NOT symmetrical.'));
end
% save the coefficients
fid = fopen('data/Coefficient_16KKaiser.dat','w');
fprintf(fid,'%i\n', coeffs_quant);
fprintf(fid,'%i\n', h_quant);
fclose(fid);
fprintf('\n\nWARNING !\n')
fprintf('This coefficient recordered in Coefficient_16K.dat are scaled to one. It is necessary then to multiply them with the number of bits needed for the quanization: 2^N_bits \n\n')
%% Compare with reference FTF coefficients
h_ftf = load('data/Coefficient_16KKaiser.gold');
d = h_ftf' == h_quant;
if sum(d) == length(d)
disp(sprintf('NOTE: quantized FIR coefficients are identical to FTF coefficients.'));
else
if reproduce_ftf
disp(sprintf('ERROR: quantized FIR coefficients are NOT identical to FTF coefficients.'));
else
disp(sprintf('NOTE: quantized FIR coefficients are NOT identical to FTF coefficients.'));
end
end
% freqz(c3,1,M2*32,200*10^6);
% q=1/(2^(W-1));
% c3=q*round(c3/q);
% c3=c3*max(c2)/max(c3);
%% Compare with LOFAR1 coefficients
h_lofar1 = load('data/Coeffs16384Kaiser-quant.dat');
d = h_lofar1' - h_quant;
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(d)
title('Difference between recreated coefficients and LOFAR1 coefficients');
xlabel('Sample')
ylabel('Coefficient')
grid on;
figure(2);
n=0:M2-1;
plot([c2',c2(M2-n)',c3']);
title('impulse response');
legend('unquantized','reversed','quantized');
d = h_lofar1' == h_quant;
if sum(d) == length(d)
disp(sprintf('NOTE: quantized FIR coefficients are identical to LOFAR1 coefficients.'));
else
disp(sprintf('NOTE: quantized FIR coefficients are NOT identical to LOFAR1 coefficients.'));
end
%% Plot (un)quantized impulse response
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot([h_fourier', h_quant']);
legend('unquantized','quantized');
title('Impulse response');
xlabel('Sample')
ylabel('Coefficient')
grid on;
%% Plot quantization error histogram to see effect of use_uencode
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
hist(h_fourier - h_quant,100);
title('Quantization error histogram');
xlabel('Quantization error');
ylabel('Relative probability');
grid on;
figure(2);
hist(c2-c3,100);
title('quantization error histogram');
xlabel('quantization error');
ylabel('relative probability');
figure(10);
freqz(c3,1,M2*4);
title('frequency response');
Nfft = 16384*8;
c3 = c3/max(c3);
%spectrum_fft = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
spectrum = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
frequency = [-Nfft/2+1:Nfft/2]*200*10^6/Nfft;
plot(frequency(length(frequency)/2:end),20*log10(spectrum(length(frequency)/2:end)))
grid on
title('Gain of the filter window with quantized coefficients (16bits) using Fircls1 window with the FFT interpolation method: blackman window')
xlabel('f(Hz)')
ylabel('Gain(dB)')
legend('FFT interpolation')
%% Plot DC response per polyphase
% . the DC response over the polyphases gets more flat with more windowing
dc_fourier = sum(reshape(h_fourier, N, L)');
dc_quant = sum(reshape(h_quant, N, L)');
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot([dc_fourier', dc_quant']);
legend('unquantized','quantized');
title('DC response per polyphase');
xlabel('Polyphase index')
ylabel('DC response')
grid on;
%% Plot frequency response of original and quantized filter
%freqz(h_fourier, 1, M3);
%freqz(h_quant, 1, M3);
fft_org = fft(h_fourier, M3);
fftshift_org = fftshift(fft_org);
spectrum_org = abs(fftshift_org) / max(abs(fftshift_org));
spectrum_org_dB = 20*log10(spectrum_org);
phase_org = unwrap(angle(fftshift_org));
fft_quant = fft(h_quant, M3);
fftshift_quant = fftshift(fft_quant);
spectrum_quant = abs(fftshift_quant) / max(abs(fftshift_quant));
spectrum_quant_dB = 20*log10(spectrum_quant);
phase_quant = unwrap(angle(fftshift_quant));
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
plot(f3_axis, spectrum_org_dB, 'r', f3_axis, spectrum_quant_dB, 'g')
legend('original','quantized')
title('Power spectrum of original and quantized filter');
xlabel('f [Hz]')
ylabel('Gain [dB]')
grid on;
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
% There is a large absolute phase offset between phase_org and phase_quant,
% due to the quantization noise. Remove this absolute phase offset to be
% able to compare the phase within the pass band, to see that it is the
% same and is linear from -0.5:0.5. Use the phase offset at DC so the center
% at index M3/2. Manually zoom in to range -1:1 to see the linear phase from
% -0.5:0.5 in the subband.
phase_offset = phase_org(M3/2) - phase_quant(M3/2);
phase_quant = phase_quant + phase_offset;
plot(f3_axis, phase_org, 'r', f3_axis, phase_quant, 'g')
legend('original','quantized')
title('Phase spectrum of original and quantized filter');
xlabel('f [Hz]')
ylabel('Phase [rad]')
grid on;
%% Determine half power bandwidth of the subband
pass_band = spectrum_quant >= 0.5;
pass_indices = find(pass_band);
pass_index_lo = pass_indices(1);
pass_index_hi = pass_indices(end);
% linear interpolate to relative index at 0.5 power
a = spectrum_quant(pass_index_lo);
b = spectrum_quant(pass_index_lo - 1);
pass_index_lo = pass_index_lo + (0.5 - b) / (a - b);
a = spectrum_quant(pass_index_hi);
b = spectrum_quant(pass_index_hi + 1);
pass_index_hi = pass_index_hi + 1 - (0.5 - b) / (a - b);
pass_bw = (pass_index_hi - pass_index_lo + 1) / M3 * fs2;
%disp(sprintf('%f', pass_index_lo));
%disp(sprintf('%f', pass_index_hi));
disp(sprintf('NOTE: quantized FIR filter normalized subband bandwidth is %6.4f * f_sub', pass_bw));
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment