4
$\begingroup$

I'm trying to model a sensor system that has an averaging behaviour. The frequency response is almost identical to a box filter and looks roughly like this:

Desired Frequency Response

Transferring this into a discrete time model would require a box filter of non-integer length - e.g. $N=2.5$ samples. Now I am looking for ways to model this system.

Here are my attempts and why they failed for me:

1. Ordinary Lowpass

As the desired frequency response has a lowpass characteristic, it would seem logical to try a lowpass filter first. However, they fail to reproduce the zero found in the desired frequency response. Also, they end in a zero at nyquist, which is not wanted.

2. Interpolated box filter

Using the impulse response $h[i] = [1, 1, f]$ where $0 < f < 1$ allows me to approximte a box filter with $N$ somewhere between 2 and 3. Here are the frequency responses of these filters for $Fs = 24kHz$ and $f = 0, 0.1, 0.2, ... , 1$:

Interpolated Box Filter

The problem is that the attenuation only approaches zero for $N=2$ and $N=3$. For anything in between it becomes way less with the worst being $N=2.5$ where the attenuation is only about -16dB.

3. Downsampled Box Filter:

I designed the desired box filter for a higher samplerate, e.g. oversampled by a factor of $S=32$. Then I lowpass-filtered it with a windowed-sinc and got these impulse responses:

Lowpass filtered box filters

I downsampled this to my original samplerate by keeping only the samples $S/2 + i*S$ and got these impulse responses: Downsampled lowpass filtered box filters

However, the frequency responses of this look very similar to the simple "interpolated" filters from attempt #2. They are so similar, that it doesn't even make sense to add another picture here. The major difference is a significantly higher computational load and an additional processing delay. Increasing the size of the windowed sinc lowpass kernel doesn't actually improve things much, it only adds additional delay due to the pre-ringing.

4. Crude oversampling

The idea was to interpolate $S$ samples for each actual sample and apply the box filter to these. I used 4-point interpolation that accounts for samples $i-1, i, i+1, i+2$ for each output sample at a position between $i$ and $i+1$. I can then re-arrange the formula to calculate the specific contribution of each input sample to the final output value like this:

h = zeros(ceil(N) + 2)
totalNumOversampledSamples = S * N
for i = 0 .. totalNumOversampledSamples:
    samplePosition = i / S
    intSamplePosition = floor(samplePosition)
    fractional = samplePosition - floor(samplePosition)
    // get interpolation coefficients for a 4pt interpolation
    a,b,c,d = getInterpolationCoefficients(fractional)
    // add those to the impulse response
    h[intSamplePosition - 1] += a
    h[intSamplePosition]     += b
    h[intSamplePosition + 1] += c
    h[intSamplePosition + 2] += d
// normalize
h /= sum(h)

(I assumed the first $S$ samples to not be interpolated to avoid adding another coefficient to the front of my impulse response) The resulting filter is quite efficient, but unfortunately, the resulting frequency response is pretty bad - probably due to the poor interpolation scheme used: Crudely interpolated box filter

5. Additional thoughts

I though of upsampling my input data, then applying an ordinary box filter to it before downsampling again. With this method, I could actually realise a "fractional length" box filter because in the upsampled domain, the box filter can be of integer length. However, this operation is entirely linear, so it should be possible to transform the same operation to an ordinary FIR filter and skip the upsampling step - which I did attempt in my 3rd approach. I am not sure why it didn't work.

Here's the actual question:

How could I model this system to fulfill these criteria:

  1. Keep the characteristic shape, especially the "zero" of the desired transfer function, or at least a high attenuation.
  2. Be able to "sweep" the zero(s) across the frequency spectrum much like it would be possible with a "moving average" filter in a continuous-time system.
  3. Keep computational load within reason (this must be able to run in real time)
  4. Phase response is not important
$\endgroup$
2
$\begingroup$

The problem might be solved already by the existing answers, but I thought I'd add my solution, which adds another degree of freedom resulting in a much closer match of the filters' magnitude responses. What I came up with is a simple system of four linear equations with the following conditions:

  1. unity gain at DC
  2. gain of the continuous-time (CT) filter at Nyquist
  3. zero at the same frequency as the CT filter

This is similar to the existing answers, but with the additional condition that the responses at Nyquist are also identical. This makes the resulting magnitude responses match each other very closely (see figure below).

As an example, I chose the width of the CT boxfilter as $T=6e-5$, and used a sampling frequency $f_s=48 \textrm{ kHz}$. The discrete-time (DT) filter has four samples because there are $4$ degrees of freedom (note that two degrees are taken by the zero at positive and negative frequencies). The result looks like this (top: magnitude responses, bottom: impulse response of DT filter):

enter image description here

Note that there's virtually no difference between the magnitude responses of the CT and DT filters.

EDIT: With this method one can incorporate an arbitrary number of zeros, which is necessary if the width of the CT impulse response becomes larger compared to the sampling period. In that case we naturally end up with a longer filter.

Here is an example for the same sampling rate as before ($f_s=48 \textrm{ kHz}$), but with a longer CT impulse response with $T=15e-5$:

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ I tried to cheat and impose unity gain at DC and Nyquist, while keeping the 3 taps. It's working, but the coefficients are complex (the middle is always 0), and the lobes go overboard inbetween the notch and the farther end of the spectrum compared to the notch, except when the notch is at half: i.stack.imgur.com/YFifD.png . I suppose it's only fair, since otherwise it would have to have been a nonlinear phase. Also, a 1st order IIR with complex xoefficients doesn't seem to work, either: the notch is ok, the gains are not. Most probably I'm doing it wrong. $\endgroup$ – a concerned citizen Apr 16 at 17:48
  • $\begingroup$ @aconcernedcitizen: I think we should keep things real :) At Nyquist I think the DT filter should have the same gain as the CT filter, then their magnitudes match best. And yes, the filter I suggest has a non-linear phase, but linear phase was not a requirement. $\endgroup$ – Matt L. Apr 16 at 17:56
  • $\begingroup$ Oh, I was just goofin' around. That part with the nonlinear phase was a simple conclusion, that's all. I suppose that's what I get for cheating. :-) Still, I don't see anywhere the attenuation at Nyquist for OP. Am I missing it? $\endgroup$ – a concerned citizen Apr 16 at 18:02
  • 1
    $\begingroup$ @LoveDaOOP: You could use the same method if you have more than one zero in the frequency interval $[0,F_s]$. For each zero you get two more equations, i.e., two more filter taps. So if necessary you could easily force the DT filter to have the exact same notches as the CT filter. As for your other attempts, I don't know, I didn't take a close look (mainly because I thought the solution I had in mind was more straightforward and would yield a better result anyway). $\endgroup$ – Matt L. Apr 17 at 10:22
  • 1
    $\begingroup$ @LoveDaOOP: I added an example with a longer CT impulse response and exact implementation of the notches. $\endgroup$ – Matt L. Apr 17 at 16:16
2
$\begingroup$

I don't know if this is what you really want but, inspired by your 2nd attempt, I thought about the triangular window, which doesn't have to have the ends null (like Bartlett), but which has variable zeroes on the unit circle, depending on the taper, and I simply tried using h=[f, 1, f], with f=[0.5 : 0.1 : 10], and here are the results (normalized gain):

test

It approaches Nyquist/2 asymptotically and the lobe at Nyquist keeps on rising, while if $f\leq0.5$, the roots are real, and you have no notch. But, since it's a simple, three taps filter, the definig equation is a quadratic, and since you need a clear notch, then you can impose the roots to have a specific transfer function. It also looks like you don't need it to be "too lowpassy", that is, the attenuation towards Nyquist is pretty awful, so, for example, if you need a notch at $\frac34$ Nyquist, then all you have to do is:

$z=-\exp\left(\pm j\pi(1-\frac{3}{4})\right) \space => \space (x-z_1)(x-z_2)=x^2+\sqrt2x+1 => [1, \sqrt2, 1]$

3/4

Which can be normalized, or not. If you need $0.8\frac{f_0}{2}$, thenL

$$z=-\exp\left(\pm j\pi(1-0.8)\right)=>[1,1.618,1]$$

0.8

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ It's possible you missed on this, but the last two examples show how precise the method is: first you impose a zero at 3/4 normalized frequency, then the graph has the X-axis divided in 4, and the zero falls exactly on the 3rd mark. The 2nd example shows a zero at 0.8, and the X-axis is divided into 10 parts, the zero falling right on the 8th. I didn't show a log axis for the magnitude because that wouldn't have made it clear where the zero was at (it was meant to be a visual representation, after all, I presumed you know that the log magnitude would have been the same as what you have). $\endgroup$ – a concerned citizen Apr 15 at 18:48
  • 1
    $\begingroup$ @DanBoschen No! Your approach may have ended in the same point as mine, but you had another way of viewing it from the beginning! Please leave the answer, it's just another way of reaching the goal. $\endgroup$ – a concerned citizen Apr 15 at 18:56
  • $\begingroup$ @LoveDaOOP turns out my answer is actually identical to this one that citizen had already given-- so I deleted my solution. A comment for this post is the coefficients can be solved using c1,c3=1, c2=−2cos(ω), or divide all by 2(1-cos(ω)) to normalize the scaling independent of frequency. $\endgroup$ – Dan Boschen Apr 15 at 18:56
  • 1
    $\begingroup$ @DanBoschen Please leave the answer back, now I feel guilty for saying that. $\endgroup$ – a concerned citizen Apr 15 at 18:59
  • $\begingroup$ @citizen It looks like you are doing the exact same thing (placing the zeros). And you showed how to normalize so it's all there! $\endgroup$ – Dan Boschen Apr 15 at 18:59
1
$\begingroup$

that has an averaging behaviour.

this looks more like a notch filter to me, designed to notch out 16 kHz or so. (maybe to block CRT line line frequency ?)

I would try to model this as a notch filter at the notch frequency with a Q of about $sqrt(2)$ or thereabouts. this appears to be an analog filter, so I'd use a fairly high sample rate to avoid any bilinear distortion at high frequencies or do an impulse invariant transform from an analog prototype.

that has an averaging behaviour.

Could also be reflection. The way to model this would not be a box filter but an impulse response with a value of 1 at time 0 and 1 at a time that corresponds to a half a period at 16 kHz, i.e. 31.25 $\mu s$. That's 3 samples at 96 kHz or 1.5 samples at 48 kHz sample rate. The latter one would require a fractional delay.

So which one is it ?

The transfer function for a single reflection at 31 $\mu s$ or a notch of Q=0.7 at 16 kHz look pretty much identical in the range that you show in your graph. Even the phase response is more or less identical in your frequency range of interest.

Fractional delays are hard to do, so I would start with a notch filter.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hilmar, thanks for the input. It is actually an averaging behaviour (magnetic pickup with a finite gap size that averages the "magnetic information". The Notch is actually a great idea, as typically only a single zero will be in the frequency response. I could combine this with a 1 pole lowpass to model the part above the zero. $\endgroup$ – LoveDaOOP Apr 15 at 13:04
  • $\begingroup$ Average to me means low pass. This clearly comes up again after the notch, so I think it's a stretch to call it averaging $\endgroup$ – Hilmar Apr 15 at 14:27

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.