Added saved plots. Added MATLAB tool requirements. Updated Figure labels. After review by Gijs.

applications/lofar2/model/pfs_coeff_final.m

@@ -21,7 +21,10 @@
 % Authors :
 %   Created by FilterTaskForce team 2004 for LOFAR1 [1, 2]
 %   Updated in 2022 for LOFAR20 by E. Kooistra
-% Purpose:
+% Purpose: Generate and analyse subband FIR filter coefficients and compare
+%          with LOFAR1 coefficients [4, 5]
+%
+% MATLAB: Tested with MATLAB R2018a + Signal Processing Toolbox
 %
 % References:
 %   [1] hdl/applications/lofar1/FilterTaskForce/readme_FilterTaskForce.txt
@@ -93,6 +96,7 @@
 % 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
+% j) Save plots in plots/
 
 close all;
 clear all;
@@ -101,6 +105,7 @@ fig = 0;
 reproduce_ftf = true;
 %reproduce_ftf = false;
 
+
 %% User settings
 if not(reproduce_ftf)
     % passband ripple (in dB);
@@ -152,13 +157,14 @@ M2 = N*L;
 % fine frequency resolution
 P = 10;
 M3 = M2 * P;
-% frequency scale
+% frequency scales, normalized to f_sub = 1
 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
 
+
 %% Plot possible window functions
 fig=fig+1;
 figure('position', [300+fig*20 200-fig*20 1000 800]);
@@ -169,6 +175,8 @@ title('Window functions');
 xlabel('Sample')
 ylabel('Gain')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 
 %% Compute filter
@@ -180,6 +188,7 @@ 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'));
@@ -196,6 +205,7 @@ else % default 'none'
 end
 h_wcomp = h_window' .* h_comp;
 
+
 %% Plot frequency response of compute filter
 fft_wcomp = fft(h_wcomp);
 fftshift_wcomp = fftshift(fft_wcomp);
@@ -209,9 +219,11 @@ 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]')
+xlabel('f [Hz] (normalized)')
 ylabel('Gain [dB]')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 fig=fig+1;
 figure('position', [300+fig*20 200-fig*20 1000 800]);
@@ -219,12 +231,15 @@ 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]')
+xlabel('f [Hz] (normalized)')
 ylabel('Phase [rad]')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 
 %% Use fourier interpolation to create final filter
+% . interpolation is used because fircls1 algorithm cannot handle >> 1000 taps
 f2 = zeros(1, M2);
 % copy the lower frequency half of the computed filter into f2.
 n = 0:M1/2;
@@ -253,6 +268,7 @@ else
     h_quant = round(h_fourier);
 end
 
+
 %% Load LOFAR1 coefficients
 if use_lofar1
     h_quant = load('data/Coeffs16384Kaiser-quant.dat');
@@ -270,11 +286,13 @@ else
   disp(sprintf('WARNING: quantized FIR coefficients are NOT symmetrical.'));
 end
 
-% save the coefficients
+
+%% Save the coefficients
 fid = fopen('data/Coefficient_16KKaiser.dat','w');
 fprintf(fid,'%i\n', h_quant);
 fclose(fid);
 
+
 %% Compare with reference FTF coefficients
 h_ftf = load('data/Coefficient_16KKaiser.gold');
 d = h_ftf' == h_quant;
@@ -288,6 +306,7 @@ else
     end
 end
 
+
 %% Compare with LOFAR1 coefficients
 h_lofar1 = load('data/Coeffs16384Kaiser-quant.dat');
 d = h_lofar1' - h_quant;
@@ -297,8 +316,10 @@ figure(fig);
 plot(d)
 title('Difference between recreated coefficients and LOFAR1 coefficients');
 xlabel('Sample')
-ylabel('Coefficient')
+ylabel('Coefficient difference')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 d = h_lofar1' == h_quant;
 if sum(d) == length(d)
@@ -318,6 +339,9 @@ title('Impulse response');
 xlabel('Sample')
 ylabel('Coefficient')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
+
 
 %% Plot quantization error histogram to see effect of use_uencode
 fig=fig+1;
@@ -328,6 +352,8 @@ title('Quantization error histogram');
 xlabel('Quantization error');
 ylabel('Relative probability');
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 
 %% Plot DC response per polyphase
@@ -344,6 +370,8 @@ title('DC response per polyphase');
 xlabel('Polyphase index')
 ylabel('DC response')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 
 %% Plot frequency response of original and quantized filter
@@ -369,9 +397,11 @@ 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]')
+xlabel('f [Hz] (normalized)')
 ylabel('Gain [dB]')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 fig=fig+1;
 figure('position', [300+fig*20 200-fig*20 1000 800]);
@@ -387,9 +417,11 @@ 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]')
+xlabel('f [Hz] (normalized)')
 ylabel('Phase [rad]')
 grid on;
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')
 
 
 %% Determine half power bandwidth of the subband
@@ -408,3 +440,5 @@ 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));
+file_name = sprintf('plots/pfs_coeff_final_%d.jpg', fig);
+print(file_name, '-djpeg')