Bandpass filter that automatically adapts its bandwidth when a transient is detected (to avoid to smoothen the transient)

Question

Let's say we want to isolate a band 1000 hz +/- 50 hz.

Obviously, limiting the bandwidth by applying a passband filter will always destroy a bit the sharp transients (a Dirac or a rectangular envelope / Heavyside step function requires all frequencies, so if limit the bandwidth, we lose a part of it, and it becomes smoother).

Question: are there some adaptative band-pass filters that would auto-extend their bandwidth for a short time when a transient is detected, in order to not lose the sharp transients?

In this example (1000 hz sinusoid modulated by a rectangle envelope, input in blue):

the filter would of course still focus on the 1000hz +/- 50hz band, but it would extend its bandwidth near the transient so that the transient is not smoothened like with a normal filter (signal in red).

Does such an adaptative bandpass exist, and is it available easily in most languages (Matlab, Python, etc.)?

NB: on this graph there is nothing else except the 1000hz sinusoid, so you may wonder "why bandpass filtering?", but it's just an example, in the general case, it would be a broadband signal.

Answer 1

Ok, I think I got your point. You want a BPF, $H(z)$, that auto extends its bandwidth accordingly to the energy distribution in the magnitude spectrum. If you have a pure 1k Hz sinusoidal tone (that corresponds, in the frequency domain, to a dirac delta located at $\omega_0=\pm2\pi 1$k rad/s), you want to pass only frequencies in the 1k$\pm 50$ Hz range, and if you have a transient event with a white noise-like distribution, you want an all-pass filter to preserve the sharp attack.

What you need is a resonator filter [1]: $$ H(z)=\frac{(1-\lambda)\sqrt{1+\lambda^2-2\lambda\cos(2\omega_0)}}{1-(2\lambda\cos(\omega_0))z^{-1}+\lambda^2 z^{-2}}, $$

its behavior for different values of $\lambda \in [0,1]$ is like so:

so for $\lambda\to 0$ you will get a flat response to catch transient events, and for $\lambda\to 1$ you will have a localized filter at the desired frequency. Here, for illustration purposes, I set $w_0=\pi/2$, but you can change the desired frequency using the formula $w_0=2\pi F_0/F_s$.

For setting $\lambda$ automatically you can use the spectral flatness estimator [2]:

$$ f = \frac{\left(\prod_{n=0}^{N-1}{x[n]}\right)^{1/N}}{\frac{1}{N}\sum_{n=0}^{N-1}{x[n]}}, $$

which is $f=1$, when the magnitude spectrum is completely flat, and $f=0$, when the magnitude spectrum is completely localized. Therefore, you can make $\lambda=1-f$. I wrote the following code to exemplify how you can apply this control:

Fs=16e3;
F0=1e3;
w0 = 2*pi*F0/Fs;
x1 = [zeros(1,50),2*rand(1,50)-1];
x2 = 0.7*sin(w0.*[1:100])+0.3*rand(1,100);
x3 = 0.7*sin(3.5*w0.*[1:100])+0.3*rand(1,100);
y = [adaptiveResonatorFilter(x1,w0), adaptiveResonatorFilter(x2,w0), adaptiveResonatorFilter(x3,w0)];
plot([x1,x2,x3],'linewidth',2)
hold on
plot(y,'linewidth',2)
xlabel('Samples')
ylabel('Amplitude')
legend('Original','Filtered')

function y = adaptiveResonatorFilter(x,w0)
    X = fft(x);
    mX = abs(X);
    mX = mX/max(mX);
    sf = mean(mX,'g')/mean(mX,'a')
    lambda = ifelse(0.5<1-sf, 0.99, 0.0)
    B = (1-lambda)*sqrt(1+lambda^2-2*lambda*cos(2*w0));
    A = [1,-2*lambda*cos(w0), lambda^2];
    [H,W] = freqz(B,A,linspace(-pi,pi,length(mX)));
    Y = X .* fftshift(H);
    y = real(ifft(Y));
end

which gives the following output:

where you can see that transient part is kept untouched, the 1k Hz pure tone contaminated with noise has been cleared and the 3.5k Hz pure tone has been attenuated, as you wanted.

Note: I am taking this as the definition of "transient attack". Please correct me if I misunderstood.

---

1. M. Vetterli, P. Prandoni. Signal Processing for Communications. EPFL press.
2. https://en.wikipedia.org/wiki/Spectral_flatness

Answer 2

Additional resource: this is JFonseca's code translated in Python:

import numpy as np
import matplotlib.pyplot as plt
import scipy, scipy.stats, scipy.signal

def adaptiveResonatorFilter(x,w0):
    X = np.fft.fft(x)
    mX = abs(X)
    mX = mX / max(mX)
    sf = scipy.stats.gmean(mX) / np.mean(mX)
    l = 0.99 if 0.5<1-sf else 0.0
    B = [(1-l)*np.sqrt(1+l**2-2*l*np.cos(2*w0))]
    A = [1, -2*l*np.cos(w0), l**2]
    return scipy.signal.lfilter(B, A, x)

Fs = 16000
F0 = 1000.0
w0 = 2*np.pi*F0/Fs

x1 = np.concatenate((np.zeros(50), 2*np.random.rand(50)-1))
x2 = 0.7*np.sin(w0*np.arange(100))+0.3*np.random.rand(100)
x3 = 0.7*np.sin(3.5*w0*np.arange(100))+0.3*np.random.rand(100)
x = np.concatenate((x1, x2, x3))
y = np.concatenate((adaptiveResonatorFilter(x1,w0), adaptiveResonatorFilter(x2,w0), adaptiveResonatorFilter(x3,w0)))

#x = np.concatenate((np.zeros(int(0.01*Fs)), np.sin(w0*np.arange(int(0.01*Fs)))))
#y = adaptiveResonatorFilter(x, w0)

plt.plot(x)
plt.plot(y)
plt.show()