diff --git a/applications/apertif/matlab/run_pfir_coeff.m b/applications/apertif/matlab/run_pfir_coeff.m
index eaa885061b09faa395a8916ac985e4f702612fb9..e0783f1d4970938b68537a270fcc188cf42972d1 100644
--- a/applications/apertif/matlab/run_pfir_coeff.m
+++ b/applications/apertif/matlab/run_pfir_coeff.m
@@ -60,6 +60,19 @@
% depends on coeff_w r_pass ~= 1 + 0.5/2^coeff_w.
%
% See pfir_coeff_adjust.m for more information.
+%
+% FIR file management:
+% Use svn status -q data/ to see whether rerun of this script recreates
+% the same FIR file as stored in SVN data/
+%
+% Run diff compare scripts to verify saved FIR files with stored FIR
+% files:
+% >. $RADIOHDL/libraries/dsp/filter/src/python/diff_lofar_coefs
+% >. $RADIOHDL/libraries/dsp/filter/src/python/diff_pfir_coefs
+%
+% Run recreate scripts to (re)create FIR stored mif files from dat files:
+% >. $RADIOHDL/libraries/dsp/filter/src/python/recreate_4wb_mifs
+% >. $RADIOHDL/libraries/dsp/filter/src/python/recreate_pfir_mifs
clear all;
close all;
@@ -212,18 +225,21 @@ fs = 800e6/1024; % Apertif subband sample rate = 800.0 MHz / 1024
f_chan = fs/N; % channel rate = subband rate or channel rate dependent on fs setting
% FIR filter analysis
-hf_abs = abs(fftshift(fft(coeff / sum(coeff), NL)));
+ZP = 1;
+ZP = 16; % use zero padding factor ZP > 1 to increase frequency resolution of hf_abs
+coeff_zero_padded = [coeff / sum(coeff), zeros(1, (ZP-1)*NL)];
+hf_abs = abs(fftshift(fft(coeff_zero_padded, NL*ZP)));
pp = 1:N; % polyphases index
hi = 1:NL; % coefficients index
hx = hi / N; % coefficients axis in tap units
-fi = hi; % frequency index
-fx = (fi - NL/2-1) / L; % frequency axis in channel units after fft_shift()
+fi = 1:NL*ZP; % frequency index
+fx = (fi - NL*ZP/2-1) / L; % frequency axis in channel units after fft_shift()
-fi_0 = NL/2+1; % frequency index of center of channel so fx = 0
-fi_n = fi_0 - round(L/2); % frequency index of -f_chan/2
-fi_p = fi_0 + round(L/2); % frequency index of +f_chan/2
+fi_0 = NL*ZP/2+1; % frequency index of center of channel so fx = 0
+fi_n = fi_0 - round(L*ZP/2); % frequency index of -f_chan/2
+fi_p = fi_0 + round(L*ZP/2); % frequency index of +f_chan/2
hf_abs_0 = db(hf_abs(fi_0)); % power gain at center of channel
hf_abs_n = db(hf_abs(fi_n)); % power gain at -f_chan/2
@@ -253,6 +269,17 @@ if strcmp(config.design, 'fircls1')
end
+%% Create common file_name_prefix
+file_name_prefix = 'run_pfir_coeff_m_';
+if config.hp_adjust
+ file_name_prefix = [file_name_prefix, 'flat_hp_'];
+end
+if config.dc_adjust
+ file_name_prefix = [file_name_prefix, 'no_dc_'];
+end
+file_name_prefix = [file_name_prefix, config.design]
+file_name_prefix = sprintf([file_name_prefix, '_%dtaps_%dpoints_%db'], L, N, coeff_w)
+
%% Save quantized integer FIR filter coefficients
% N.B. - read '.txt'-file with 'wordpad' or use \r\n
if coeff_w>0
@@ -271,15 +298,7 @@ if coeff_w>0
'DC gain per polyphase = %f\n'] , ...
max(coeff), max(int_coeff), q_max, q_max-max(int_coeff), coeff_dc_gain));
end
- file_name_prefix = 'data/run_pfir_coeff_m_';
- if config.hp_adjust
- file_name_prefix = [file_name_prefix, 'flat_hp_'];
- end
- if config.dc_adjust
- file_name_prefix = [file_name_prefix, 'no_dc_'];
- end
- file_name_prefix = [file_name_prefix, config.design];
- file_name = sprintf([file_name_prefix, '_%dtaps_%dpoints_%db.dat'], L, N, coeff_w);
+ file_name = ['data/', file_name_prefix, '.dat'];
fid = fopen(file_name , 'w');
fprintf(fid ,'%d\n', int_coeff);
fclose(fid);
@@ -366,19 +385,41 @@ if strcmp(application, 'lofar_subband') && config.dc_adjust==true
end
% Plot FIR-filter amplitude characteristic using fft()
+% . full plot in [dB]
fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
-plot(fx, db(hf_abs)); % db() = 20*log10 for voltage
-xlo = fx(1);
-xhi = fx(NL);
-xlo = -3;
-xhi = 3;
-%xlim([xlo xhi]);
+plot(fx/ZP, db(hf_abs)); % db() = 20*log10 for voltage
+xlo = fx(1)/ZP;
+xhi = fx(NL*ZP)/ZP;
+%xlo = -1*ZP;
+%xhi = 1*ZP;
+xlim([xlo xhi]);
+ylo = -120;
+yhi = 10;
+ylim([ylo yhi]);
grid on;
-title(['FIR filter transfer function for ', config.design, strNL]);
+title(['Full FIR filter transfer function for ', config.design, strNL]);
xlabel('Frequency [channels]');
ylabel('Magnitude [dB]');
+file_name = ['plots/', file_name_prefix, '_fir_transfer_function.jpg'];
+
+% . zoomed plot in power
+fig=fig+1;
+figure('position', [300+fig*20 200-fig*20 1000 800]);
+figure(fig);
+plot(fx/ZP, hf_abs.^2);
+xlo = fx(1);
+xhi = fx(NL*ZP);
+xlo = -0.7;
+xhi = 0.7;
+xlim([xlo xhi]);
+grid on;
+title(['Zoomed FIR filter transfer function for ', config.design, strNL]);
+xlabel('Frequency [channels]');
+ylabel('Magnitude power');
+file_name = ['plots/', file_name_prefix, '_fir_transfer_function.jpg'];
+print(file_name, '-djpeg')
% Plot FIR-filter amplitude and phase characteristic using freqz()
fig=fig+1;
@@ -398,10 +439,10 @@ fig=fig+1;
figure('position', [300+fig*20 200-fig*20 1000 800]);
figure(fig);
subplot(2,1,1)
-hf_abs1 = circshift(hf_abs, L, 2); % shift to next channel (in dimension 2 because hf_abs is a row vector)
+hf_abs1 = circshift(hf_abs, L*ZP, 2); % shift to next channel (in dimension 2 because hf_abs is a row vector)
plot(fx, db(hf_abs), 'k', fx, db(hf_abs1), 'r'); % db() = 20 * log10 for voltage
-xlo = -1;
-xhi = 2;
+xlo = -1*ZP;
+xhi = 2*ZP;
xlim([xlo xhi]);
grid on;
title(['FIR filter transfer function between channels for ', config.design, strNL]);
@@ -412,8 +453,8 @@ subplot(2,1,2)
plot(fx, hf_abs.^2, 'r.', ...
fx, hf_abs1.^2, 'b.', ...
fx, hf_abs.^2 + hf_abs1.^2, 'k');
-xlo = -1;
-xhi = 2;
+xlo = -1*ZP;
+xhi = 2*ZP;
xlim([xlo xhi]);
%ylim([0.8 1.2]);
grid on;