Automated search for half power bandwidth adjustment hp_factor.

applications/apertif/matlab/run_pfir_coeff.m

@@ -41,8 +41,11 @@
 %     to choose fir1.
 %   - The half power channel bandwidth is adjusted via hp_factor to have
 %     flat power gain between channels. The hp_factor depends on L and
-%     r_pass. The hp_factor was determined empirically by rerunning this
-%     script for different values and .
+%     r_pass. The hp_factor is determined automatically when hp_search =
+%     true. For floating point coeff_w=0 the hp_factor can achieve an
+%     arbitrary accurate solution. For fixed point coeff_w>0 the solution
+%     can only be as good as possible, so then decreasing hp_resolution
+%     does not change the coeff anymore.
 %   - For fircls1 the passband ripple and stopband ripple can be specified.
 %     For fir1 the passband is flat and the stopband attenuation increases
 %     more and more for more distant frequencies. The fir1 can achieve
@@ -66,7 +69,9 @@ application = 'apertif_subband';    % calculate Apertif subband filterbank coeff
 
 %% - Application specific settings
 % Default settings
-hp_factor = 1;                                          % no channel half power bandwidth adjustment
+hp_factor = 1;
+hp_search = false;                                      
+hp_search = true;                                       % enable automatic channel half power bandwidth adjustment
 if strcmp(application, 'test_bypass')
     N = 64;
     L = 8;
@@ -105,38 +110,51 @@ else
     if strcmp(config.design, 'fir1')
         r_pass = 0;  % not used for fir1
         r_stop = 0;  % not used for fir1
-        % Set channel half power bandwidth adjustment hp_factor (determined empirically)
-        if L==16, hp_factor = 1.049922; end
-        if L== 8, hp_factor = 1.10073; end
-        if L== 4, hp_factor = 1.2; end
     else  % fircls1
         % Pass and stop band deviation for fircls1
         r_pass = 1e-4; % choose 0.0001 to have more flat reponse between channels
                        % for 0.001 it is better to choose fir1
         r_stop = 1e-4; % choose 0.0001 to avoid artefact at ends of impulse response
                        % that occur for 0.001
-        % Set channel half power bandwidth adjustment hp_factor (determined empirically)
-        if r_pass == 1e-3
-            if L==16, hp_factor = 1.054992; end
-            if L== 8, hp_factor = 1.1138; end
-            if L== 4, hp_factor = 1.24;  end
-        elseif r_pass == 1e-4
-            if L==16, hp_factor = 1.055837; end
-            if L== 8, hp_factor = 1.1154; end
-            if L== 4, hp_factor = 1.205; end
-        end
-        % config parameters that are specific for fircls1
+                       
+        % Config parameters that are specific for fircls1
         %config.design_flag = 'trace';
         config.design_flag = '';
         %config.interpolate = 'resample';
         %config.interpolate = 'fourier';
         config.interpolate = 'interpft';
     end
-    % Adjust channel half power bandwidth, so need db(sqrt(2)) = -3 dB gain instead of half gain
+    
+    % Optional FIR filter bandwidth adjust
+    if hp_search==false
+        % No channel half power bandwidth adjustment
         BWchan = hp_factor/N;
 
         % Calculate the FIR coefficients
         coeff = pfir_coeff(N, L, BWchan, r_pass, r_stop, coeff_w, config);
+    else
+        % Automatically adjust channel half power bandwidth, so need db(sqrt(2)) = -3 dB gain instead of half gain
+        hp_radix = 2;  % choose 2 for binary search, or 10 for 'decimal' search
+        hp_step = 1/hp_radix;
+        hp_resolution = 0.000001;
+        while true
+            BWchan = hp_factor/N;
+            coeff = pfir_coeff(N, L, BWchan, r_pass, r_stop, coeff_w, config);
+            hf_abs = abs(fftshift(fft(coeff / sum(coeff), N*L))); 
+            fi_0 = N*L/2+1;              % frequency index of center of channel, so 0 Hz
+            fi_p = fi_0 + L/2;           % frequency index of edge of channel at +f_chan/2
+            hf_abs_p = hf_abs(fi_p);     % gain edge of channel at +f_chan/2
+            disp(sprintf('hp_factor = %10.8f, hf_abs_p = %10.8f, hp_step = %10.8f\n', hp_factor, hf_abs_p, hp_step));
+            if hf_abs_p < sqrt(0.5)
+                hp_factor = hp_factor + hp_step;
+            else
+                if hp_step <= hp_resolution, break; end
+                hp_factor = hp_factor - hp_step;
+                hp_step = hp_step/hp_radix;
+                hp_factor = hp_factor + hp_step;
+            end
+        end
+    end
 end
 NL = N*L;                % Total number of FIR filter coefficients (taps)