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.
- M. Vetterli, P. Prandoni. Signal Processing for Communications. EPFL press.
- https://en.wikipedia.org/wiki/Spectral_flatness
scipy
functionscipy.signal.iirfilter
. But the problem would be the same with a FIR designed with firwin. $\endgroup$ – g6kxjv1ozn Sep 26 '18 at 12:31