4
$\begingroup$

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:

taps

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

filtered_signal

$\endgroup$
  • 1
    $\begingroup$ can you add the parameters you've used for your Remez? I agree with Laurent, per definition, a real bandpass wouldn't pass DC, but maybe you're referring to a complex bandpass that has a passband between $f_{cutoff,low}<0<f_{cutoff,upper}$ ? $\endgroup$ – Marcus Müller Aug 19 '16 at 11:05
  • $\begingroup$ Yes, I will add the code generating the filter to the question $\endgroup$ – SomeRandomPhysicist Aug 19 '16 at 11:15
  • $\begingroup$ @Marcus Müller D* right, I too often forget about complex filters $\endgroup$ – Laurent Duval Aug 19 '16 at 11:22
  • 1
    $\begingroup$ If I sum the answers up, your question seems more or less answered. Offset change is normal, and you can compensate. However, I feel like we do not know enough to reach the actual goal $\endgroup$ – Laurent Duval Aug 20 '16 at 16:58
  • 1
    $\begingroup$ Indeed my question is answered, as you surmised, from these answers I now know that offset change is normal and predictable from the coefficients and can be adjusted for in filtering so that I produce an attenuated signal centered the same as the input signal. Thank you everyone who answered my question and discussed it for your help. $\endgroup$ – SomeRandomPhysicist Aug 20 '16 at 17:44
5
$\begingroup$

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".

$\endgroup$
4
$\begingroup$

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:

Remez band-pass filter design

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.

$\endgroup$
  • $\begingroup$ Ok, so, as I understand it, scipy.signal.remez shouldn't be producing a band-pass filter with a negative offset? $\endgroup$ – SomeRandomPhysicist Aug 19 '16 at 10:47
  • 1
    $\begingroup$ If you apply a real band-pass filter, by definition the filtered signal should be zero-average, because you "add" an initial offset to the sine (which is normalize positive and negative) to cast it between $0$ and $2^{16}-1$ $\endgroup$ – Laurent Duval Aug 19 '16 at 11:27
  • $\begingroup$ Laurent, @MarcusMüller, please see my answer below for an explanation of the first and last taps; it's normal (under these circumstances). $\endgroup$ – Matt L. Aug 19 '16 at 18:29
  • $\begingroup$ @Matt L. I need to fully understand it first. A single-sample (symmetric) echo is still not clear to me, yet my faculties in FIR design are very poor. $\endgroup$ – Laurent Duval Aug 19 '16 at 18:36
  • 1
    $\begingroup$ @LaurentDuval: Just compute the frequency response corresponding to the first and last sample: $\delta[n]+\delta[n-N]$ The corresponding frequency domain magnitude is (the magnitude of) a cosine, and this is what we see in the (stop band) of the frequency response. $\endgroup$ – Matt L. Aug 19 '16 at 18:39
4
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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

frequency response

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

not getting better

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?

slightly more benign

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.

$\endgroup$
  • 1
    $\begingroup$ Thanks for your help, I was planning to do as you suggest and use a low-pass filter to remove the high frequency components and then use a slower clock frequency for the band-pass filter such that I can get much better attenuation in the region of interest. What would you suggest to use as a "normal" windowing FIR designer? I've been using scipy.signal.remez and matlab's fdatool which seem to produce very similar filters. (as found also by @Laurent Duval in his comment above) $\endgroup$ – SomeRandomPhysicist Aug 19 '16 at 13:21
  • $\begingroup$ @SomeRandomPhysicist I often find myself in a situation where I can work with a complex baseband (=lowpass) signal just as well as with a real passband signal, so what I often do is: 1. mutliply the real signal by an $e^{j2\pi f_\text{shift}}$, rotating everything in frequency domain so that the center of my passband ends up at $f=0$, and then 2. filtering with a simple real low pass filter, eg. with a Kaiser Window, while simultaneously reducing the sample rate (note that you can now reduce the sample rate twice as much, because you're having a complex signal!). $\endgroup$ – Marcus Müller Aug 19 '16 at 13:29
  • $\begingroup$ @SomeRandomPhysicist the nice thing about that is that there's often filter implementations that include the option to rotate (and otherwise, multiplying with $e^{j...}$ isn't that bad, either, especially when you can put that in a Look-up-table) while filtering and allow for simultaneous decimation – a decimating FIR is especially nice, because when implemented correctly, it can run at the output sampling rate, instead of the input rate, which makes much longer filters feasible! $\endgroup$ – Marcus Müller Aug 19 '16 at 13:32
  • $\begingroup$ By the way, the same works for bandpass filters, too: You can reduce the sampling rate (==decimate), as long as the resulting rate still fulfills Nyquist's demands. The resulting signal will simply be aliased into your first Nyquist zone. So don't be shy to implement longer filters. $\endgroup$ – Marcus Müller Aug 19 '16 at 13:39

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.