2
$\begingroup$

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):

enter image description here

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.

$\endgroup$
  • $\begingroup$ It is confusing to me to know how you obtained the red curve. Did you apply the filter in the frequency domain by multiplying the spectrum by a shifted and scaled rectangular window or you did you convolve the time signal by a sinc? $\endgroup$ – JFonseca Sep 25 '18 at 19:11
  • $\begingroup$ @WaltzForZizi I just applied an IIR filter (such as Butterworth) to the input, with ready-to-use scipy function scipy.signal.iirfilter. But the problem would be the same with a FIR designed with firwin. $\endgroup$ – g6kxjv1ozn Sep 26 '18 at 12:31
  • $\begingroup$ Can we limit your application to the detection of 1000 Hz sinusoidal bursts? Is an overall delay of the filtered result acceptable? (By delay, I mean that the sharp transient would be intact but the overall waveform is delayed in time?). What is the primary purpose for you of needing to maintain the transient (in case there are other solutions to your purpose)? $\endgroup$ – Dan Boschen Sep 27 '18 at 13:03
  • $\begingroup$ @DanBoschen Yes to limit it to 1000 Hz sinusodial content (mixed with some background noise). Other sinusoidal contents are present too, but at 1050 Hz, 1200 Hz, etc. Sharp transient intact + overall delay is ok (if I can reshift the signal to compensate the delay). It's a sound for which I need to isolate some harmonics. I would like to maintain transients to keep the "fast" attack of the sound. $\endgroup$ – g6kxjv1ozn Sep 27 '18 at 17:15
  • $\begingroup$ I would not necessarily refer to this as the optimum solution but would be tempted to use this bandpass filter (since the bandwidth is easily tunable): dsp.stackexchange.com/questions/40482/… and combine that with a correlator to optimally detect the presence of the 1KHz signal. The greater the sensitivity for the correlator/detector, the longer the delay- so as long as the delay is not an issue I think this could work quite well for you. Specifically I would use the tuneable bandpass with widest bandwidth in default mode, and then .... $\endgroup$ – Dan Boschen Sep 27 '18 at 17:24
2
$\begingroup$

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: enter image description here

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:

enter image description here

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
$\endgroup$
  • $\begingroup$ Yes, the goal is to 1) filter out what is outside of the range [950hz, 1050hz] 2) except when there is a transient (sudden envelope change), because in this case all frequencies are required to create the sharp attack (if the transient is bandpass filtered too, it will be smoothened which is not wanted!). Would you have a code example for your proposal? Would be great to test on a real signal how it works! $\endgroup$ – g6kxjv1ozn Sep 26 '18 at 19:35
  • $\begingroup$ @g6kxjv1ozn, please see the new answer, I added a resonator filter that is also tunable like the LPF filter proposed by DanBoschen. $\endgroup$ – JFonseca Oct 4 '18 at 3:08
  • $\begingroup$ Thank you @JFonseca. How would you implement it on the example signal (blue curve in the question)? By small blocks of 100 samples each (re-computing the lambda for each block)? Also what is the easiest way to find the actual coefficients of the filter (FIR or IIR) from this H(z) with Python or Matlab? $\endgroup$ – g6kxjv1ozn Oct 4 '18 at 6:52
  • $\begingroup$ @g6kxjv1ozn, I added a sample code, with a plot in the format you requested, to illustrate the algorithm and the filtering procedure. By nature, the filter is IIR because the transfer function has coefficients in the denominator. Please mark as answer if it works for you. $\endgroup$ – JFonseca Oct 7 '18 at 2:56
2
$\begingroup$

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()

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.