2
$\begingroup$

I am using a DSP processor to sample signals at 108kHz and I only want to get the DC part out of it. Such a high sampling rate is used because I want to use oversampling to reduce quantization noises.

Two filters are to be designed:

  • One has 10Hz stop-band frequency; new data are read every 0.1 second;
  • One has 0.1Hz stop-band frequency; new data are read every 10 second;

What I did

I used several stages of decimation and downsample the signals to 300Hz. Then I use FIRs to filter the signal:

  • For the 10Hz bandwidth: 128 stage filters are used
  • For the 0.1Hz bandwidth: 768 stage filters are used.

My question

Performance is the most important to me: The sharper low-pass filters are, the better.

  • Could I ask for suggestions in implementing such a filter?
  • Is it possible to design a 0.01Hz low pass filter? If so, the new data can be read every 100 seconds?
  • I heard that using IIRs could be a solution, is it?

Thank you very much!

Is averaging the best DC filter?

I tried to compare with matlab:

Setting 1: %Use Matlab fir1 function

size = 1080000; h = fir1(size, 0.00000000001); fvtool(h) FIR1

Setting 2: %Build an averaging FIR

size = 1080000; h = ones(1,size); h = h *1/size; fvtool(h) average

The theory should be both two settings have similar performance. However, it seems that fir1 has narrower stopband than the simple average filter. It seems that the average DC filter will include more noises.

|improve this question
$\endgroup$
  • $\begingroup$ is your DSP fixed point? how wide are the words? is it that you're worried about the behavior of your DC filter (one with $H(e^{j0})=1$) if the cutoff frequency is far, far below Nyquist? like $$ H(z) = (1-p)\frac{z}{z-p} $$ for $ 0 < 1-p \ll 1$? (or $p$ very, very close to 1.) there are some known numerical problems sometimes people get with them kinda filters. especially IIR and fixed-point. $\endgroup$ – robert bristow-johnson Mar 27 '15 at 3:57
  • $\begingroup$ Yes, the Nyquist is around 54kHz and cut-off is around 0.1Hz.. The DSP has a FPU actually meaning it has floating point functions. And 32bit double numbers are used. $\endgroup$ – richieqianle Mar 27 '15 at 4:10
  • $\begingroup$ 32 bit is not double floats. for floats, that's single-precision. 32-bit fixed-point numbers might be considered "double", since lotsa DSP is done with 16-bit parts. so you gotta worry a little about the "cosine problem" because $$ \cos(\omega_0) = \cos\left( \pi \frac{0.1}{54000} \right) \approx 1$$ and the whole information about $\omega_0$ is contained in the difference that is from 1. so you lose your accuracy in the coefficient unless you replace every occurrence of $\cos(\omega_0)$ with $$ \cos(\omega_0) = 1 - 2 \sin^2\left( \frac{\omega_0}{2} \right) $$ . $\endgroup$ – robert bristow-johnson Mar 27 '15 at 4:37
  • $\begingroup$ You are right. It is 32 bit float with a 32-bit processor.. $\endgroup$ – richieqianle Mar 27 '15 at 5:15
  • $\begingroup$ then, if you do this with an IIR (could be 2nd-order, like a Butterworth), then you really need to worry about the numerical issues (the "cosine problem") i alluded to above. $\endgroup$ – robert bristow-johnson Mar 27 '15 at 5:55
1
$\begingroup$

The definition of the Fourier Transform of a function $x(t)$ is

$$ X(f) \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t $$

Evaluating for $f=0$ you see

$$ X(0) = \int\limits_{-\infty}^{\infty} x(t) \, \mathrm{d}t $$

As you can see the value at DC equals the area of the function.

The definite integral of a function from $-\infty$ to $\infty$ is equal to the transform at the origin $X(0)$.

This may help you if you don't have streaming data. Instead of a sharp lowpass filter or taking the fourier tranform and extracting $X(0)$, just integrate your data and the result is equivalent to obtaining $X(0)$.

|improve this answer
$\endgroup$
  • 1
    $\begingroup$ hay Amnon, would you take some time to learn how to do $\LaTeX$ markup for your mathematical expressions? $\endgroup$ – robert bristow-johnson Jul 10 '18 at 19:17
  • 1
    $\begingroup$ okay, so i did it for you. see how easy that was? also had to change the name of "$f(t)$" and "$F(f)$" to "$x(t)$" and "$X(f)$" because of the conflicting use of "$f$". Welcome to the DSP stack exchange. We are looking forward to your contributions, but please use $\LaTeX$ for your math. Wikipedia has a pretty good reference for the math markup. $\endgroup$ – robert bristow-johnson Jul 13 '18 at 17:31
2
$\begingroup$

If all you want to do is get an estimate of the DC component of your data then you can simply obtain the average value of your signal over a unit of time. This can be easily shown by looking at the expression for calculating the Fourier coefficients for 'f=0' (DC component).

The easiest way to do this is via a moving average filter (http://en.wikipedia.org/wiki/Moving_average) which would work over a buffer.

In terms of an FIR filter, this looks like 'h = ones(1,Nbuff)./Nbuff' (so, 'h=[0.5,0.5]' for a 2 sample buffer). 'Nbuff' stands for the order of the filter. This can be applied directly on the buffer you obtain from your CODEC through a convolution function as provided by your processor's command set.

You could also use an 'online algorithm' (an IIR filter essentially) which maintains an updated estimate of the mean of the signal on a per-sample basis. For this purpose please see: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm You could also design a "typical" IIR filter that can be as strict as you like, provided that your processor and sampling period can accommodate the computations required.

(Of course you can adapt between working in terms of frames and working in terms of individual samples)

Just a note, I have a feeling that you are struggling with the concept of Group Delay. The cut off frequency of your filter does not exactly mean that "new data can be read every 100 seconds". New data will be fed to your filter constantly, at the sampling rate of the system, but the phase delay caused on each harmonic component of your signal between the input and the output of the system might vary. How this varies, depends on the type of the system (whether it is FIR / IIR) and its form. Here is a starting point: http://en.wikipedia.org/wiki/Group_delay_and_phase_delay

Hope this helps.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks A_A, your answer helps a lot! It seems that both IIR and FIR are performing average essentially. "new data can be read every 100 seconds", I mean: the DC signal is varying slowly actually with a defined pattern. For example the DC signal varies every 10 seconds and I know the pattern. I want to wait long enough to ensure that I get the new DC instead of the old one. Something like step response.. $\endgroup$ – richieqianle Mar 26 '15 at 3:21
  • $\begingroup$ I have edited the post, could you please have a look at the section of "is averaging the best dc filter"? $\endgroup$ – richieqianle Mar 26 '15 at 9:11
2
$\begingroup$

any low-pass filter with DC gain of 1 ( $H(1)=1$ ) is a candidate for DC extraction filter. this kind of "get DC" filter problem is pretty much equivalent and complementary to the "DC blocking filter" problem.

i don't see why you wouldn't use a far lower order IIR filter for your LPF. you can get sharp without 768 taps (and, presumably 768 multiply-accumulate instructions) using an IIR filter.

excluding (for the moment) the cost of the filter, the tradeoff is between a solid DC output and responsiveness when the DC value has actually changed. how quickly do you need your DC filter to adapt when the actual DC in the input has changed?

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks Robot! I am new to DSP and I chose FIR because I read a tutorial about it. My stated algorithm was built half year ago and now I am re-evaluating it. 1. How quick should the filter: I can decide the speed of change of DC input. It is decided in this way. To get a satisfied result, the noise bandwidth should be controlled within 0.1Hz; To get DC extraction filter with bandwith narrower than 0.1Hz I need a A stage IIR/FIR; For the IIR/FIR to work properly, I need to wait B seconds. Then the DC in the input will change every b seconds, which is controlled by myself. $\endgroup$ – richieqianle Mar 27 '15 at 2:03
  • $\begingroup$ Could you also please comment on comparison between average and IIR? parkorn and A_A seem to be correct in saying that the DC essentially is the average of incoming signal and it can also be proved by calculation of FFT. Thanks a lot! $\endgroup$ – richieqianle Mar 27 '15 at 2:05
  • $\begingroup$ from my POV, essentially by definition, DC is the average of the incoming signal. now if your signal is like a soundfile, with a beginning and an end, then the DC component is literally the average value of the soundfile. if it's real time, then DC is a moving and weighted average of some definition (whatever the impulse response is) of the most current samples of the input. if the coefficients of the weighted moving average, which are the coefficients of the impulse response $h[n]$, if they add to $1$, then $H(e^{j0})=1$. dunno who parkorn or A_A are. $\endgroup$ – robert bristow-johnson Mar 27 '15 at 3:46
  • $\begingroup$ anyway, there are lotsa different $h[n]$ such that $$\sum\limits_{n=-\infty}^{\infty} h[n] = 1$$ some of them will be for filters with faster response time (narrower impulse response) than others (fatter impulse response). again, you can make lotsa different skinny or fat impulse responses with $$H(e^{j0})=1$$ with either FIR or IIR filters. depends on what you need your DC value for. $$ $$ just learned who pakornosky and A_A are. $\endgroup$ – robert bristow-johnson Mar 27 '15 at 3:49
  • $\begingroup$ Thank you very much Robert! That is true in defining "What DC is". In my senario, there is a true DC buried in noises. So the question becomes what is the best h[n] or IIR to filter all the noises and find out the real DC... $\endgroup$ – richieqianle Mar 27 '15 at 4:18
2
$\begingroup$

so, what's wrong with a simple 1st-order, 1-pole, IIR filter, where you're a little careful with the coefficient that soooo close to 1?

$$ y[n] = (1-p) x[n] + p y[n-1] $$

where $ p = \cos(\omega_0) = \cos\left(\pi \frac{0.1}{54000} \right) $

but express the whole thing in terms of $1-p$ instead of $p$.

$$\begin{align} y[n] & = (1-p) x[n] + p y[n-1] \\ & = y[n-1] + (1-p) \left( x[n] - y[n-1] \right) \\ & = y[n-1] + 2\sin^2\left( \frac{\omega_0}{2} \right) (x[n] - y[n-1]) \\ \end{align}$$

where $ \omega_0 \ll 1 $ or $\pi$ or whatever in the same neighborhood.

|improve this answer
$\endgroup$
  • $\begingroup$ Hi Robert, thanks for your help. I tested fir1, moving average, ellip IIR filter and block average algorithm on matlab. Block average I mean take a block of data and average it without previous states. It seems that block average gives the best result. I think I will just use the averaging filter. $\endgroup$ – richieqianle Mar 27 '15 at 10:14
  • $\begingroup$ Hi Robert, I have been reading your answers again recently. Shall not IIR filters be in form like: y[k] = b(1)*x[k] + b(2)*x[k-1] - a(2)*y[k-1]. ? In your equation above, there seems to be no x[k-1]. Could you explain on that? $\endgroup$ – richieqianle Apr 1 '16 at 6:16
  • 1
    $\begingroup$ problem is that your $a(2)$ is very very close to $-1$. maybe too close for your arithmetic precision. so instead of $$ y[k] = b(1)*x[k] + b(2)*x[k-1] - a(2)*y[k-1]$$, why not $$ y[k] = b(1)*x[k] + b(2)*x[k-1] - (a(2)+1)*y[k-1] + y[k-1] $$ then your $a(2)+1$ is very close to zero and you can tweek that number accurately. oh, for a DC LPF filter, it's likely that $b(2)=0$. $\endgroup$ – robert bristow-johnson Apr 1 '16 at 6:30
  • $\begingroup$ Thanks for your reply. I agree with you on that. However, there is still $b(2)*x[k-1]$ term in the formula. But in your suggested formula there is no such term. Could you explain on this? $\endgroup$ – richieqianle Apr 1 '16 at 6:33
  • $\begingroup$ for a DC LPF? i doubt that it $b(2)$ is anything other than zero. my first equation in my answer (this is a year old, so i am not very familiar) shows a 1-pole LPF. NO zeros. if there was a zero somewhere, then $b(2)$ would not be zero. $\endgroup$ – robert bristow-johnson Apr 1 '16 at 6:35
1
$\begingroup$

The easiest way (and actually make sense) is to find the average of your input signal.
If you only need the DC component, there is nothing to do with Passband or Stopband, sharp enough or not, FIR or IIR, ...

Just remind yourself again, what do you actually want?

Extended Ans

  • Time Domain

    DC component = constant (0 Hz) = not vary in amplitude.
    Other components (freq > 0 Hz) = go positive and negative over times.

    If you calculate the mean value of the waveform, the positive and negative part will cancel each other and leave you the DC component.

    • Start with sinusoidal signal, Have you got this point?
    • Any signal can be realized as the sum of sinusoidal signals which vary in frequency and amplitude.

    So, you can get DC component from averaging the signal.

  • Frequency Domain

    Use FFT to find spectrum of signal, the first point of output is the DC component (0 Hz). This is theoretically supported by Fourier.

Evaluation
To test your approach, you need experiments.

  • Generate your known signal with DC offset.

  • Get DC part using each approach (1. filtering with LPF and 2. averaging)

  • Calculate the accuracy of DC output with a priori offset.

Try more waveform as well as the random signals.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for your reply! I just want to get the DC and get rid of the noises which are huge. I talked about Passband/Stopband/FIR/IIR because I want to minimize the noises. Do you mean: averaging the sampled signal is the best low pass filter for DC? And the only factor deciding how much noises will enter the result is number of samples? $\endgroup$ – richieqianle Mar 26 '15 at 3:26
  • 1
    $\begingroup$ Refer to DC_component "the DC bias, DC component, DC offset, or DC coefficient is the mean value of the waveform" Or this is not your point, Pls let me know. $\endgroup$ – pakornosky Mar 26 '15 at 8:30
  • 1
    $\begingroup$ I have already added the further explanation above. :) $\endgroup$ – pakornosky Mar 26 '15 at 8:51
  • $\begingroup$ Thanks for your reply. I have edited the post, could you please have a look at the section of "is averaging the best dc filter"? $\endgroup$ – richieqianle Mar 26 '15 at 9:11
  • $\begingroup$ What I actually mean is "The mean of signal = DC". This is fact and go directly to your point. What you trying to do now is to filter out others and leave the DC. So, calculating DC directly or filtering the other freqs. Hope this help. $\endgroup$ – pakornosky Mar 26 '15 at 9:36

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.