I have FIR filters cascade, shown at scheme below. All filters have same coefficients. This cascade only approximation for the problem description.
I use "rule of thumb" for word size increase in decimating filters is to increase word size 1 bit for each decimate by four (or log2(DF)/2 where DF is Decimation Factor).
Full FIR filter output bit width calculated as f_0_full = ceil(log2(sum(abs(coe)))) + f_0_inp, where coe - coefficient vector, and f_0_inp = 16 - filter #0 input bit width.
Сoefficients for filters:
Fs = 450e+6/4; % Sampling Frequency
Fpass = 20e+6; % Passband
N = 42; % Filter order
coe = firhalfband(N, Fpass/(Fs/2)); % HB FIR calculation
bw_coe = 18; % Filter coefficients bit width
coe = coe/max(abs(coe)); % Normalization
coe = round( coe*(2^(bw_coe - 1) - 1) ); % Quantization
I can calculate minimum and maximum output I can get at filter #0 output with this function, where X - vector for minimum and maximum value at input:
function [Y, bw] = get_fout_max( coe, X )
X_min = X(1, 1);
X_max = X(1, 2);
G_p = sum(coe(coe > 0));
G_m = sum(coe(coe < 0));
Y_min = G_p * X_min + G_m * X_max;
Y_max = G_p * X_max + G_m * X_min;
Y = [Y_min, Y_max];
bw = log2(max(abs(Y))) + 1;
For filter #0 X = [-2^(f_0_inp - 1), 2^(f_0_inp - 1) - 1], filter maximum is Y = [-13290137560, 13289875415], or round(Y/2^(f_0_full - f_0_out - f_0_scale)) = [-50698, 50697] after truncation, where f_0_out = 17 - target signal bit width after truncation, f_0_scale = 0 - shift for most significant (MSB) bit after truncation from the MSB of full length f_0_full = 35 (process shown at picture below)
Or I can get this absolute maximum if I will use next signal s as input for FIR #0 (bw_s = 16 - input signal s bit width)
s(coe > 0) = -(2^(bw_s - 1) - 0);
s(coe < 0) = (2^(bw_s - 1) - 1);
But how can I calculate the maximum output for multiple (two, three, etc) filters in cascade?
I cant use same method, becouse I think it is impossible to get signal at decimation FIR Filter #0 output for Filter #1 such:
s_for_filter_1(coe > 0) = -50698;
s_for_filter_1(coe < 0) = 50697;
I want to learn how to calculate scalers to avoid sign duplicate at filters outputs (with scaler shift) or how to calculate absolute maximum at each filter considering the previous filters.
Full code for cascade:
Fs = 450e+6/4; % Sampling Frequency
Fpass = 20e+6; % Passband
N = 42; % Filter order
coe = firhalfband(N, Fpass/(Fs/2)); % HB FIR calculation
bw_coe = 18; % Filter coefficients bit width
coe = coe/max(abs(coe)); % Normalization
coe = round( coe*(2^(bw_coe - 1) - 1) ); % Quantization
N = 2^20; % N samples in signal
bw_s = 16; % Signal Bit Width
s_gen_mode = 0;
switch s_gen_mode
case 1
% Max Amplitude at #0 FIR output
s = zeros(N, 1);
s(coe > 0) = -(2^(bw_s - 1) - 0);
s(coe < 0) = (2^(bw_s - 1) - 1);
case 2
% Static Max level in signal
s = -(2^(bw_s - 1) - 0) * ones(N, 1);
otherwise
% Sinusoidal signal
freqs = 0.4e+6 ; % carrier frequency
s = exp( 1i*2*pi * (0:N - 1)' * (freqs/Fs) );
% s = awgn(s, 1.761+6.02*16);
s = s / max([max(real(s), [], 'all'); max(imag(s), [], 'all')]);
s = round(s * (2^(bw_s - 1) - 1));
end
%% Filter #0
f_0_inp = 16; % Filter Input Bit Width (BW)
f_0_full = ceil(log2(sum(abs(coe)))) + f_0_inp; % Full output BW
f_0_out = f_0_inp + ceil(log2(2)/2); % 'Rule of thumb' bit growth (according to f_0_inp)
f_0_scale = 0; % Scale shift
s_f_0 = filter(coe, 1, s); % Filtering
s_f_0 = s_f_0(1:2:end, 1); % Decimation
s_f_0_trunc = round(s_f_0 / 2^(f_0_full - f_0_out - f_0_scale));% Truncation
%% Filter #1
f_1_inp = f_0_out;
f_1_full = ceil(log2(sum(abs(coe)))) + f_1_inp;
f_1_out = f_0_inp + ceil(log2(4)/2);
f_1_scale = 0;
s_f_1 = filter(coe, 1, s_f_0_trunc);
s_f_1 = s_f_1(1:2:end, 1);
s_f_1_trunc = round(s_f_1 / 2^(f_1_full - f_1_out - f_1_scale));
%% Filter #2
f_2_inp = f_1_out;
f_2_full = ceil(log2(sum(abs(coe)))) + f_2_inp;
f_2_out = f_0_inp + ceil(log2(8)/2);
f_2_scale = 0;
s_f_2 = filter(coe, 1, s_f_1_trunc);
s_f_2 = s_f_2(1:2:end, 1);
s_f_2_trunc = round(s_f_2 / 2^(f_2_full - f_2_out - f_2_scale));
%% Filter #3
f_3_inp = f_2_out;
f_3_full = ceil(log2(sum(abs(coe)))) + f_3_inp;
f_3_out = f_0_inp + ceil(log2(16)/2);
f_3_scale = 0;
s_f_3 = filter(coe, 1, s_f_2_trunc);
s_f_3 = s_f_3(1:2:end, 1);
s_f_3_trunc = round(s_f_3 / 2^(f_3_full - f_3_out - f_3_scale));
%% Filter #4
f_4_inp = f_3_out;
f_4_full = ceil(log2(sum(abs(coe)))) + f_4_inp;
f_4_out = f_0_inp + ceil(log2(32)/2);
f_4_scale = 0;
s_f_4 = filter(coe, 1, s_f_3_trunc);
s_f_4 = s_f_4(1:2:end, 1);
s_f_4_trunc = round(s_f_4 / 2^(f_4_full - f_4_out - f_4_scale));
%% Filter #5
f_5_inp = f_4_out;
f_5_full = ceil(log2(sum(abs(coe)))) + f_5_inp;
f_5_out = 16; % Required output bit width
f_5_scale = 5;
s_f_5 = filter(coe, 1, s_f_4_trunc);
s_f_5 = s_f_5(1:2:end, 1);
s_f_5_trunc = round(s_f_5 / 2^(f_5_full - f_5_out - f_5_scale));
s_f_5_bw = log2(max([abs(real(s_f_5_trunc)); abs(imag(s_f_5_trunc))], [], 'all')) + 1;
fprintf('Output signal bit width %.3f
', s_f_5_bw)
fprintf('Recommended scale for last filter %i
', f_5_scale + (f_5_out - ceil(s_f_5_bw)))
%% Spectrum
data = s_f_5_trunc;
figure(1); periodogram(data, blackmanharris(size(data, 1)), 'centered', 2^(nextpow2(size(data, 1))+1), Fs/(2^6));
grid minor
ylim([-100, 70])
question from:
https://stackoverflow.com/questions/65941002/maximum-possible-output-for-fir-filters-cascade 
