Negative result from filtering a positive signal with a band-pass filter

Question

I have used scipy.signal.remez to calculate the coefficients for a band-pass filter and when I use it to filter a sinusoidal signal that goes between 0 and some positive number (e.g. $2^{16}$) I get that the filtered signal is negative ( it is attenuated, as I expected ) but it has a negative offset, is this expected? I need my outputs to be between 0 and $2^{16}$ (since I am implementing this filter in an FPGA and send the filtered signal to a DAC which accepts positive numbers between 0 and $2^{16}$ Should I just deal with this by adding a positive offset to re-centre the attenuated signal?

Code used to generate the FIR filter:

fs = 2500000         # Sample rate, Hz
band = [95000, 105000]  # Desired pass band, Hz
trans_width = 5000    # Width of transition from pass band to stop band, Hz
numtaps = 200       # Size of the FIR filter.

edges = [0, band[0] - trans_width,
     band[0], band[1],
     band[1] + trans_width, 0.5*fs]
taps = scipy.signal.remez(numtaps, edges, [0, 1, 0], Hz=fs)

The taps are real and look like this:

Here is an example of the result of using this FIR filter to filter a 200KHz sine wave in python:

Answer 1

What you see is what one would expect. As pointed out in Marcus Müller's answer, your band pass filter has a relatively poor stop band attenuation, and since DC is in the filter's stop band, DC is not sufficiently attenuated. You can easily predict what is going to happen:

1. compute the filter's DC gain, which is just the sum over all filter coefficients: $H(0)=\sum_nh[n]$ This will be a not so small negative number.

2. Take the DC offset of your input signal ($2^{15}$), and multiply it with $H(0)$. This must be the DC offset of your output signal.

So the DC offset of the output signal can be easily predicted, even without computing the output!

What can you do to get rid the DC offset of the output signal? Since you can perfectly predict it, you can just subtract the DC offset $=2^{15}\cdot H(0)$ from your output signal. Another fun solution is to increase the filter length by $1$ sample, and set that new last sample to the value $-H(0)$ (as computed above). This will make sure that $H(0)$ of the new filter is exactly $0$, and it will only insignificantly change the overall frequency response of the filter.

For the ones who wonder if it is normal that the first and last sample of the impulse response are much larger than all other samples: yes, this is normal, especially for narrow band pass filters, and it is a direct consequence of the equi-ripple optimality criterion implied by the Remez algorithm. These two impulses (first and last sample) are echos generating the sinusoidal (i.e., equi-ripple) stop band behavior in the frequency domain. The ripple in the impulse response generates the (impulse-like) narrow pass band in the frequency domain. So everything can be explained by the basic Fourier relationship "sinusoid $\Longleftrightarrow$ impulse".

Answer 2

Ideally, the gain of "strict" (real) band-pass filter should be close to $0$ near the DC at $0$ Hz, or at least be sensibly lower than $1$. Otherwise, it would be more a low-pass. Hence the DC or "average" value of a positive signal is lowered. This is apparent in the example from the doc:

In other words, if you apply a real band-pass filter to your, by definition the filtered signal should be zero-average (on the long run), because you initially "add" an initial offset to the sine (which is normally positive and negative) to cast it between $0$ and $2^{16}−1$.

Another potential issue is quantization: if you design a filter with "real" coefficients, and carelessly quantize its coefficients, the resulting filter might have a frequency behavior different from that of the designed one. This however is not so problematic with FIR filters, if the quantization is fine enough.

Answer 3

You are essentially asking for a band-pass filter that also preserves the DC offset. Literally, that would be a combination of a band-pass filter and a low-pass filter, such that the gain is 1 both at DC and at the frequency you are filtering for.

However, if you know what your DC offset is (e.g. $2^{15}$), it is probably much better to subtract this offset before filtering and add it after filtering. That way, you are not as dependent on the exact filter characteristics at DC.

Answer 4

I ran your code and made a magnitude response plot:

import numpy
import scipy.signal
from matplotlib import pyplot
try:
    import seaborn
except: 
    pass #graphs won't be as pretty *shrug*

fs = 2500000         # Sample rate, Hz
band = [95000, 105000]  # Desired pass band, Hz
trans_width = 5000    # Width of transition from pass band to stop band, Hz
numtaps = 200       # Size of the FIR filter.

edges = [0, band[0] - trans_width,
     band[0], band[1],
     band[1] + trans_width, 0.5*fs]
taps = scipy.signal.remez(numtaps, edges, [0, 1, 0], Hz=fs)

def plotme(taps, filename):
    pyplot.figure(figsize=(10,5), dpi=90)
    omega, h = scipy.signal.freqz(taps)
    pyplot.plot(omega, numpy.abs(h))
    pyplot.xlim( ( 0, numpy.pi ) )
    pyplot.ylim( ( 0, max(numpy.abs(h)) ) )
    pyplot.tight_layout()
    pyplot.savefig(filename)

plotme(taps, "/tmp/freq.png")

As you can see, at DC ($f=0$), your filter has a suppression of but 0.2 – that is but -14dB (power)!

Now, I noticed that these first and last taps look strange. So what happens when we modify your taps vector to exclude these?

plotme(taps[1:-1], "/tmp/freq2.png")

Now the overall frequency response looks much better, but the DC component got worse. However, you might solve that by cascading with a very relaxed high-pass filter with a cut-off somewhere between DC and your passband.

Hm, what about we set the zeroth tap to - the last tap?

Things look slightly more benign for the DC component, but we see the overall suppression has greatly decreased! My impression here is that Remez tries to find a minimum for overall attenuation error, and the DC leakage is the price you pay for acceptable overall attenuation.

My advice here is not to use the specific Remez algorithm; a "normal" windowing FIR designer usually can produce single-passband bandpasses just as well by designing an equivalen low-pass first and then applying spectral shift magic.

Let's do a rough analysis of how much headroom you have for proper attenuation in you 200 taps:

So, Belanger says:

$$N\approx \frac 23 \log_{10} \left[\frac1{10 \delta_1\delta_2}\right]\,\frac{f_s}{\Delta f}\tag{*}$$

with

$$\begin{align} f_s &\text{ the sampling rate}\\ \Delta f& \text{ the transition width,}\\ \delta_1 &\text{ the ripple in passband,}\\ &\text{ ie. "how much of the original amplitude can you afford to vary"}\\ \delta_2 &\text{ the suppresion in the stop band}. \end{align}$$

Let's find $\delta_2$. $(*)$ implies that, assuming $\delta_1=\delta_2$ is appropriate,

$$\begin{align} 10\cdot 10^{\frac 32 N\frac{\Delta f}{f_s}} &\approx \delta_2^{-2}\\ &\text{plugging in your numbers}\\ \delta_2^{-2} &\approx 10\cdot10^{\frac 32 200\frac{5\cdot10^3}{2.5\cdot 10^6}}\\ &=10\cdot 10^{0.6}\\ &\approx 40\\ &\implies \delta_2 \approx 0.15 \end{align}$$

So, your real limiting factor is really just that your filter is too short to get something with a $\frac{f_s}{500}$ transition width done with nice attenuation.

The for that is probably a decimating stage before you do your bandpassing – if the highes frequency you're interested in is much much lower than half your sampling rate, you should just use a low pass filter to "roughly" cut out a lower bandwidth first, decimate in the process (decimate by $N$: throw away $N-1$ samples of every $N$ samples; filtering is necessary, because that would usually imply aliasing), and then design your band pass on the resulting lower sampling rate.