Low pass filtering in software

[Simon]

I'm currently low pass filtering my ADC signals with the current formula:

Code: [Select]

average = (average * (n-1) + x ) / n    ;

Does this approximate an RC filter (1st order 3db) ? if it does how would I calculate the frequency response. The values of "x" are delivered at regular intervals.

The ADC itself does permit me to accumulate values and right shift the result to perform averaging. this is not the same calculation but does it amount the same thing? I need 2 differently filtered streams of my readings so a hardware filter it out (I actually just took it out!).

[SiliconWizard]

This is an exponential moving average, where alpha = 1/n.
You would get much lower exec time if n is a power of 2, as you can just right-shift instead of dividing, which is always expensive for arbitrary values (although the compiler may change it to a multiplication if 'n' is known at compile time).
Alternatively, you could use floating point if you have access to a FPU and using fp otherwise makes sense.

If you want to know more: https://en.wikipedia.org/wiki/Exponential_smoothing

EMA:
y = alpha*x + (1 - alpha)*y

Your case:
y = (1/n)*x + ((n-1)/n)*y = (1/n)*x + (1 - 1/n)*y
so it's the same thing with alpha = 1/n

This article gives you the frequency and impulse responses: https://blog.mbedded.ninja/programming/signal-processing/digital-filters/exponential-moving-average-ema-filter/

[Simon]

aha, thank you. So I am on the right track. Unfortunately no FPU so the best I can do is take my 16 bit ADC value and put it into a 32 bit variable shifted up by 8 bits. so I have 8 bits on the right to act as digits after the decimal place and 7 bits on the left that the multiplication can happen into. unfortunately the values represent a signed quantity with an offset so I have to subtract that on each reading. I don't know how long an abs() in C takes as I really don't need the sign, that would mean that yes bit shifting unsigned values may work but I guess it's all the same with the hardware divider.

I do have a hardware divider that will take no more than 32 cycles and I'm running at 48MHz / 21ns cycles.

I don't know

[Siwastaja]

Also, if using integers and not floating point, remember to multiply x by value similar to (or same as) n; and divide down the average before using with the same factor. Otherwise you are getting a lot of quantization loss.

Why? Think about what happens with n = 1000 and x = 57 (let's assume average is 0 initially):

average = (0 * 999 + 57) / 1000 = 57 / 1000 = 0
Repeat this as many times as you want, but result stays 0, when it should be slowly approaching 57.

Fix this by:
average = (0 * 999 + 57*1000) / 1000 = 57000 / 1000 = 57
(use as average/1000 = 57/1000, which still rounds to 0, but hey, now average is at value 57 so next round gets bigger and bigger...

And as SiliconWizard says, use powers of two for efficient division-by-shift:

average_accumulator = (average_accumulator*255 + 256*x ) / 256;
cur_average = average_accumulator/256;

[gf]

if it does how would I calculate the frequency response. The values of "x" are delivered at regular intervals.

The EMA  y=alpha*x+(1-alpha)*y  is basically a Direct Form I implementation of the transfer function  H(z)=alpha/(1+(alpha-1)*z^-1).

In Octave (or Matlab), you can plot the frequency response in the following way:

Code: [Select]

pkg load signal
alpha = 0.25;      % put your alpha here
fs = 1000;         % put your sample rate here
b = alpha          % numerator coefficients of transfer function
a = [ 1 alpha-1 ]  % denominator coefficients of transfer function
[H,f] = freqz(b, a, 10000, fs);
plot(f, 20*log10(abs(H)))
grid on
xlabel("frequency")
ylabel("dB")

EDIT: Goggle also revealed this document
https://tttapa.github.io/Pages/Mathematics/Systems-and-Control-Theory/Digital-filters/Exponential%20Moving%20Average/Exponential-Moving-Average.html
which even gives closed-form expressions for the frequency response and cutoff frequency.

[Nominal Animal]

I don't know how long an abs() in C takes

On ARMv6-m (Cortex-M0/M0+) and ARMv7e-m (Cortex-M3/M4/M7) it's two or three instructions, no jumps.

If you have <stdlib.h> or equivalent, just use its abs().
If you use GCC or Clang, use __builtin_abs().  This should yield the fastest integer absolute value the compiler can generate.

(Unless there is an even better way, it usually ends up becoming t = x arithmetic right-shifted all the way so that sign bit is copied to all bits, then abs(x) = (x+t)^t = (x^t)-t.)

Here's the code I would probably use on a Cortex-M:

Code: [Select]

const int32_t  average_scale;  // 1 .. 32767, constexpr for C++
int32_t  average_magnitude_scaled;  // 0 .. 1073741823
uint32_t  average_magnitude;  // 0 .. 32768

uint32_t  update_average_magnitude(const int16_t sample) {
    const int32_t  sample_magnitude = abs(sample);
    average_magnitude_scaled += (int32_t)average_scale * (sample_magnitude - (int32_t)(average_magnitude_scaled >> 15));
    average_magnitude = (uint32_t)(average_magnitude_scaled + 16384) >> 15;
    return average_magnitude;
}
It compiles to less than two dozen instructions, no conditionals; and uses a single 32×32=32-bit (signed) multiplication.  The rest are binary operations, additions, and subtractions.  You could definitely use multiple in parallel, too, even on a Cortex-M0.

For the mathematical analysis, average_scale = 32768 - 32768 * alpha = 32768 - 32768 / n, or equivalently, alpha = 1 - scale/32768 = 1/n, and n = 32768 / (32768 - scale).  Essentially, alpha is the weight of the new sample, and (1-alpha) the weight of the average thus far.

Note that in my implementation, average_magnitude is rounded to average_magnitude_scaled.  This is mostly because the original sample range is signed, and I want zero samples to have the same unit range as all other magnitudes, i.e. ]-0.5,+0.5[ instead of double the others ]-1.0,+1.0[ that it would be if truncating.

(To get the scaled average magnitude to converge to the middle of the integer-valued region, we do need to truncate the scaled average magnitude when dividing it down in the average calculation.  Fortunately, on Cortex-M, this generates even better code than using pre-computed value in a separate variable, because adding or subtracting a value shifted by some bits is only one instruction.  If a separate variable were used, even average_magnitude, it would first have to be loaded into a register, thus only adding more machine code.)

The smallest alpha this can represent is 1/32768 = 0.000030517578125, corresponding to n = 32768, when average_scale = 1. Then, convergence from zero to sample magnitude 32768 takes 359605 consecutive samples of value -32768.  average_magnitude does reach 32768 after 343220 consecutive samples, but it takes an additional 16385 samples for the average_magnitude_scaled to converge and no longer change.

Note that you cannot choose an arbitrary n or alpha, because you're limited to integer values of average_scale (and therefore integer values of 32768/n).  This, the reduced choice in larger n or very small alpha, is the price we pay by converting the integer division to an integer multiplication.

[bson]

Note that the low-pass property comes from the addition, not the magnitude adjustment.  If you for example add up 256 8-bit values into a 16-bit sum, then this sum represents a low-pass filter.  Any magnitude adjustment is purely to avoid overflow, but if you're fine with a 16-bit sum then there's no reason make a magnitude adjustment.  It's the summing, not the division.

[Nominal Animal]

Yes, the division is not relevant mathematically.  At the core, at sample \$y_i\$, the exponential average is
$$A_i = y_i \alpha + y_{i-1} \alpha (1 - \alpha) + y_{i-2} \alpha (1 - \alpha)^2 + \dots + y_0 \alpha (1 - \alpha)^i = \alpha \sum_{k=0}^i y_k (1-\alpha)^{i - k}$$
and the division is only used so that we can express the above using integer operations, via \$\alpha = (n - 1) / n\$ or \$\alpha = 1 - \operatorname{scale} / 32768\$.


Another often used filter is moving average, the simple average (sum!) over N consecutive samples.  It isn't useful here, because as a low-pass filter, it has very slow roll-off; it is best for removing noise while keeping a sharp step response.  Anyway, to calculate this, you need a (circular) buffer to remember the last N samples:

    uint16_t  circular_buffer[N] = { 0 };
    uint_fast16_t  circular_buffer_index = 0;
    uint32_t  circular_buffer_sum = 0;

When a new sample is added, you subtract the old magnitude from the sum, add the new magnitude to the sum, and store the new magnitude into the buffer:

    circular_buffer_sum -= circular_buffer[circular_buffer_index];
    circular_buffer_sum += value;
    circular_buffer[circular_buffer_index] = value;
    if (++circular_buffer_index >= N)
        circular_buffer_index = 0;

It's obviously extremely fast to calculate, just one addition and one subtraction per sample –– although you do need to also keep N last samples in memory.

The downside is that circular_buffer_sum is scaled by N, i.e. the average is circular_buffer_sum/N, or (circular_buffer_sum + N/2)/N if rounded.  For power of two N this simplifies to a bit shift, but many compile-time constant N can be calculated by 32×32=high-32-bit multiplication (by reciprocal).  (Use a constexpr in C++ and C23 and later, or a preprocessor macro in ANSI C/C99/C11/C17, for N.)

The amplitude frequency response for this filter at sample wavelength \$\lambda\$, \$1 \le \lambda \le N/2\$, is
$$H(\lambda) = \left\lvert \frac{\sin(\pi \, \lambda)}{N \sin\left(\frac{\pi \, \lambda}{N}\right)} \right\rvert$$
or, in relative frequency form, \$0 \lt f \le 0.5\$,
$$H(f) = \left\lvert \frac{\sin(\pi \, f \, N)}{N \sin(\pi \, f)} \right\rvert$$

You can use my FIR analysis tool page to look at the frequency response; just use as many ones as you have N.

[radiolistener]

I have a question regarding FIR filters.

I'm using FIR LPF filter with normalized kernel coefficients but it produces overshoot peaks for square wave:


The input signal is a square wave with range -1...+1. While output signal has 1.18545 overshoot peaks.

I tried to test different test waveforms and it produces different overshoot.

Is there a reliable way to estimate max possible peak overshoot on the filter output from FIR kernel coefficients?

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

[BrianHG]

I can look for one of my old designs where I created an audio 16 bit in, 16 bit out bass and treble filter if anyone is interested.  Basically a 1st order 3db filter.  It used only 16bit adds and one 16x8 bit input multiplier, 24 bit out.  It had 2x8 bit controls, 1 for the series R value in the low pass filter for the treble effect, 1 for the resistor load across the output series capacitor for the high pass filter, or DC filter.  As a bonus, it had 6 channel mixer inputs.

Though it was written in Verilog, so a little decoding of the language would be needed.

Oh, a cheap way to see the response of your filter is to feed it random noise, like a few kb worth and run the output through an FFT.  You should see a nice response curve.  No filter should give you a flat line so long as your random number generator is truly white.

[tggzzz]

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

Simple. Ask an LLM. You know how to do that since you have posted LLM generated C code elsewhere.

Whether you can trust the response is a separate issue. I suppose you could - as some lawyers have done - ask the LLM if its response is correct 

[Simon]

I have a question regarding FIR filters.

I'm using FIR LPF filter with normalized kernel coefficients but it produces overshoot peaks for square wave:


The input signal is a square wave with range -1...+1. While output signal has 1.18545 overshoot peaks.

I tried to test different test waveforms and it produces different overshoot.

Is there a reliable way to estimate max possible peak overshoot on the filter output from FIR kernel coefficients?

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

I have seen waveforms like that before when I was playing around constructing a square wave from sine waves. My conclusion at the time was that there is no such thing as a true square wave in nature by the rules of Fourier transforms as they have a massive peak on each edge, the more sines you add the higher the peaks but the less the oscillation.

The other thing to look at is what is the sampling frequency compared to the signal frequency. I knocked up a spreadsheet to do some filtering experiments, I just put a sine wave with an offset in. When I set the sine frequency to be just off the sampling frequency the output started to go right down to 0 slowly and then came back up despite the main amplitude being 1 with a up to 0.1 sine wave added to it.

[radiolistener]

Simple. Ask an LLM. You know how to do that since you have posted LLM generated C code elsewhere.

I already did it on different AI models, but they all have no idea how to fix it.

I have seen waveforms like that before when I was playing around constructing a square wave from sine waves.

In my case, this happens because the original signal is an ideal square wave, and the FIR filter removes some high-frequency harmonics from it. The waveform is correct and acceptable, without these overshoots, the signal would be incorrect. However, the question is how to estimate this overshoot from the FIR coefficients in order to avoid sample out-of-range issues. 

As I can see, well-known audio converters like sox also struggle with this issue and require manual gain correction before applying the filter. 

[tszaboo]

I'm currently low pass filtering my ADC signals with the current formula:

Code: [Select]

average = (average * (n-1) + x ) / n    ;

Does this approximate an RC filter (1st order 3db) ? if it does how would I calculate the frequency response. The values of "x" are delivered at regular intervals.

The ADC itself does permit me to accumulate values and right shift the result to perform averaging. this is not the same calculation but does it amount the same thing? I need 2 differently filtered streams of my readings so a hardware filter it out (I actually just took it out!).

I seriously suggest coming up with the filter response first, and programming it second. I've seen it too many times, where this naiive programming method is used, with random coefficients, random length, without thinking what you are actually trying to resolve. New microcontrollers are fast enough to run very long FIR or IIR filters, and they have hardware accelerators and optimized libraries to do it. When you write it yourself, it's not optimal, and it can lead to very nasty things like excessive noise or unstable IIR filters.
And then you either don't know why it's doing it, or no way of inspecting it with standard evaluation methods.

[tggzzz]

I'm currently low pass filtering my ADC signals with the current formula:

Code: [Select]

average = (average * (n-1) + x ) / n    ;

Does this approximate an RC filter (1st order 3db) ? if it does how would I calculate the frequency response. The values of "x" are delivered at regular intervals.

The ADC itself does permit me to accumulate values and right shift the result to perform averaging. this is not the same calculation but does it amount the same thing? I need 2 differently filtered streams of my readings so a hardware filter it out (I actually just took it out!).

I seriously suggest coming up with the filter response first, and programming it second. I've seen it too many times, where this naiive programming method is used, with random coefficients, random length, without thinking what you are actually trying to resolve. New microcontrollers are fast enough to run very long FIR or IIR filters, and they have hardware accelerators and optimized libraries to do it. When you write it yourself, it's not optimal, and it can lead to very nasty things like excessive noise or unstable IIR filters.
And then you either don't know why it's doing it, or no way of inspecting it with standard evaluation methods.

Yes and no.

First determine the system's requirements.

Second determine whether a conventional analogue filter or its digital equivalent is necessary, as opposed to being one possibility.

If not, consider alternatives that are easy to implement in the digital domain, even if they have no analogue domain equivalent. Then work out whether they are sufficient.

Obvious examples are adding random noise where the noise generation mechanism is probably not important, or N-path filters for exceptionally high Q bandpass or bandstop filters or mixers.

Thus while I agree that twiddling coefficients is suboptimal, I would ask whether Simon needs to worry about how closely it needs to approximate an RC filter.

[benseno]

(i am not an expert in this field)
In analog design the overshoots are controlled through choice of filter type, order .... The same or similar means should be available in the digital counterpart:
https://www.engineersgarage.com/types-of-filter-responses/

[tggzzz]

I have a question regarding FIR filters.

I'm using FIR LPF filter with normalized kernel coefficients but it produces overshoot peaks for square wave:


The input signal is a square wave with range -1...+1. While output signal has 1.18545 overshoot peaks.

I tried to test different test waveforms and it produces different overshoot.

Is there a reliable way to estimate max possible peak overshoot on the filter output from FIR kernel coefficients?

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

I have seen waveforms like that before when I was playing around constructing a square wave from sine waves. My conclusion at the time was that there is no such thing as a true square wave in nature by the rules of Fourier transforms as they have a massive peak on each edge, the more sines you add the higher the peaks but the less the oscillation.

And that is the rather obvious way for radiolistener to determine the maximum "overshoot".

Simple Fourier formulae state that a square wave is the sum of an infinite series of sine waves. Is that reliable "in nature"? Consider that the dual question is whether it is possible to have a single frequency. If the sine wave is not infinite duration, then there will be different frequency components due to it being turned and off.

[benseno]

I have a question regarding FIR filters.

I'm using FIR LPF filter with normalized kernel coefficients but it produces overshoot peaks for square wave:


The input signal is a square wave with range -1...+1. While output signal has 1.18545 overshoot peaks.

I tried to test different test waveforms and it produces different overshoot.

Is there a reliable way to estimate max possible peak overshoot on the filter output from FIR kernel coefficients?

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

I have seen waveforms like that before when I was playing around constructing a square wave from sine waves. My conclusion at the time was that there is no such thing as a true square wave in nature by the rules of Fourier transforms as they have a massive peak on each edge, the more sines you add the higher the peaks but the less the oscillation.

And that is the rather obvious way for radiolistener to determine the maximum "overshoot".

Simple Fourier formulae state that a square wave is the sum of an infinite series of sine waves. Is that reliable "in nature"? Consider that the dual question is whether it is possible to have a single frequency. If the sine wave is not infinite duration, then there will be different frequency components due to it being turned and off.

In analog circuitry just using RC lpf to filter square wave will not result in any overshoot, aka simple averaging filter in the digital domain.

[radiogeek381]

I have a question regarding FIR filters.

I'm using FIR LPF filter with normalized kernel coefficients but it produces overshoot peaks for square wave:

Can you provide some details?

Like

1. what was the sample rate vs. the square wave frequency? (Especially, how many samples per cycle of your square wave?)
2. How many taps are in the FIR filter?
3. Was this implemented as a streaming (time domain, sequential sample) filter?
4. If you made multiple calls to the filter routine, did it save any samples between calls?
5. How were the FIR coefficients calculated?
6. Is it *really* an FIR filter?

[SiliconWizard]

If you want to avoid "overshoots", yes one kind of low-pass filtering that guarantees that is an averaging filter (like moving average or EMA).

This all boils down to the step response of the filter.
For the basics: https://en.wikipedia.org/wiki/Step_response

[Nominal Animal]

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal.

To examine the effect of Gibbs' phenomenon, just apply the step function.

The definition of output sample \$i\$ with a finite impulse response filter with \$N\$ coefficients \$c_0\$ through \$c_{n-1}\$ applied to input samples \$x\$ is
$$y_i = \sum_{k=0}^{N-1} x_{i-k} c_k$$
A step function is a single edge of a pulse, being zero for all negative arguments, and exactly 1 for all nonnegative arguments.  Thus, when \$m\$ samples of the step function are within the filter window,
$$y_i = \sum_{k=0}^{m-1} x_{i-k} c_k = \sum_{k=0}^{m-1} c_k, \quad 0 \le m \le N$$
since those \$x\$ are all 1, and all others zero (or corresponding filter coefficient would be zero).

The peak in this response is simply \$\max \lvert y_i \rvert\$ for all \$m\$ and \$i\$.

If we consider also the Gibbs' phenomenon on trailing edges, it turns out that we want to find out the maximum magnitude of the sum in any consecutive subsequence of the filter coefficients, i.e.
$$\max \lvert \sum_{k=i_0}^{i_1-1} c_k \rvert, \quad 0 \le i_0 \lt i_1 \le N$$

Thus, if you have \$N\$ coefficients in your FIR filter, you can find the maximum peak by calculating the \$(N+1)N/2\$ sums, and picking the one with the largest magnitude.  This is the peak value relative to an edge from 0 to 1 or 1 to 0 you can get.

[radiolistener]

Can you provide some details?

Like

1. what was the sample rate vs. the square wave frequency? (Especially, how many samples per cycle of your square wave?)

sample rate is 1536 kHz, square wave frequency is 800 Hz. It's about 1920 samples per square wave period.

2. How many taps are in the FIR filter?
5. How were the FIR coefficients calculated?
6. Is it *really* an FIR filter?

FIR LPF is designed in Octave using b = fir1(N, Wp, kaiser(N + 1, 6.181877));

Here is FIR response:

3. Was this implemented as a streaming (time domain, sequential sample) filter?
4. If you made multiple calls to the filter routine, did it save any samples between calls?

it is implemented as streaming filter and has internal buffer to store previous samples.

[radiolistener]

To examine the effect of Gibbs' phenomenon, just apply the step function.

The definition of output sample \$i\$ with a finite impulse response filter with \$N\$ coefficients \$c_0\$ through \$c_{n-1}\$ applied to input samples \$x\$ is
$$y_i = \sum_{k=0}^{N-1} x_{i-k} c_k$$
A step function is a single edge of a pulse, being zero for all negative arguments, and exactly 1 for all nonnegative arguments.  Thus, when \$m\$ samples of the step function are within the filter window,
$$y_i = \sum_{k=0}^{m-1} x_{i-k} c_k = \sum_{k=0}^{m-1} c_k, \quad 0 \le m \le N$$
since those \$x\$ are all 1, and all others zero (or corresponding filter coefficient would be zero).

The peak in this response is simply \$\max \lvert y_i \rvert\$ for all \$m\$ and \$i\$.

That was my first attempt to estimate it.
But when I use classic pulse response [1,0,0,0,...], I get Amax = 0.0285229364172188 which is close to zero and don't exceed -1...+1 range.

When I use pulse response for [1,-1,-1,-1,...] it gives Amax = 1.17101652962875, which is more close to 1.18 overshoot value which I get with square wave, but still smaller than it.

I also tried to experiment with different square waves on the input and found that overshoot depends on the phase and period of square wave and can be a little bit higher, up to Amax=1.21 or something like that.

[radiolistener]

Here is a small test code in Octave/Matlab which shows issue with signal overshoot on the FIR output:

Code: [Select]

pkg load signal;

Fs = 48000*32;     % Sample rate [Hz]
Fpass = 22000;     % Passband [Hz]

Wp = Fpass / (Fs / 2);  % Normalized passband
beta = 6.181877;        % Kaiser window beta (60 dB stopband attenuation, 0.01 dB passband ripple)
N = 3040;               % FIR order (estimated with ceil((A - 8) / (2.285 * pi * Wt)))

kernel = fir1(N, Wp, kaiser(N + 1, beta));

% Normalize FIR kernel
kernel = kernel / sum(kernel);
fprintf('sum(kernel):      %g\n', sum(kernel));
fprintf('sum(abs(kernel)): %g\n', sum(abs(kernel)));

% Show FIR response
freqz(kernel, 1, 2^16, Fs);


N = 2^14;               % Signal length
%signal = cos(800 * 2*pi*(0:N-1)/Fs);        % sine wave
signal = square(800 * 2*pi*(0:N-1)/Fs);     % square wave

% Apply zero padding to the signal to see details include filter phase delay
signal = [signal zeros(1, length(kernel))];
N = length(signal);

% Apply FIR filter to the signal
response = filter(kernel, 1, signal);

Amax = max(abs(response));
fprintf("Amax = %g\n", Amax);



% ===Plot waveform===
% X axis in samples
time_min = 0;
time_max = N;
time_step = (time_max-time_min) / 5;     % grid time step (5 divs per figure)
if time_step >= 1000000000
    xscale = 1/1000000000;
    xunit = 'GigaSample';
elseif time_step >= 1000000
    xscale = 1/1000000;
    xunit = 'MegaSample';
elseif time_step >= 1000
    xscale = 1/1000;
    xunit = 'kiloSample';
else
    xscale = 1;
    xunit = 'Sample';
end
xmult = 1;
amp_pkpk = range(signal);                   % Amplitude peak-to-peak
scale = 10^(ceil(-log10(amp_pkpk)) + 5);    % 5 digits for rounding
ampl_min = round((min(signal) - amp_pkpk * (sqrt(2) - 1) / 2) * scale) / scale;
ampl_max = round((max(signal) + amp_pkpk * (sqrt(2) - 1) / 2) * scale) / scale;
figure;
plot((0:N-1)*xmult, signal,  'b', 'LineWidth', 2, 'DisplayName', 'Signal');
hold on;
plot((0:N-1)*xmult, response, 'r', 'LineWidth', 2, 'DisplayName', 'Response');
xlabel(xunit);
ylabel('Amplitude');
ylim([ampl_min ampl_max]);
xlim([time_min time_max])
grid on;
grid minor;
xticks = get(gca, 'XTick');
xticklabels = arrayfun(@(x) sprintf('%.0f', x*xscale), xticks, 'UniformOutput', false);
set(gca, 'XTickLabel', xticklabels);
legend('show');
title(sprintf('signal: min=%g, max=%g\nresponse: min=%g, max=%g', min(signal), max(signal), min(response), max(response)));


[Nominal Animal]

If you want the maximum response from any input consisting only of ±1, construct one that matches the signs of the input filter in reverse order.  Then, the maximum amplitude is the sum of the magnitudes of the filter coefficients,
$$\sum_{k=0}^{N-1} \lvert c_k \rvert$$
and occurs when the input signs match the corresponding coefficient signs.

This has basically nothing to do with Gibbs' phenomena, and everything to do with maximal filter response amplitude, which for FIR filters is the sum of the magnitudes of the filter coefficients.