DSPRelated.com
Blogs

Linear-phase DC Removal Filter

Rick LyonsMarch 30, 200826 comments

This blog describes several DC removal networks that might be of interest to the dsprelated.com readers.

Back in August 2007 there was a thread on the comp.dsp newsgroup concerning the process of removing the DC (zero Hz) component from a time-domain sequence [1]. Discussed in that thread was the notion of removing a signal's DC bias by subtracting the signal's moving average from that signal, as shown in Figure 1(a).

DC Removal Figure 1

This article is available in PDF format for easy printing

Figure 1.

At first I thought such a DC removal scheme would have the very poor passband peak-peak ripple performance of roughly 13 dB because of the 13 dB stopband ripple of a moving averager. But upon modeling this DC removal scheme I found that its passband behavior was better than I had expected. My modeling involved a computationally efficient D-point moving averager (MA), shown in Figure 1(b). The bottom path, in Figure 1(b), is a simple delay line having a length equal to the averager's group delay, (D-1)/2 samples. This enables us to time-synchronize the averager's v[n] output with the x[n] input in preparation for the subtraction operation.

The D-point moving averager (MA) in Figure 1(b) has a transfer function defined by

DC Removal Equation 1

The second ratio in Eq. (1) is merely a recursive running sum comprising a D-length comb filter followed by a digital integrator. There are two delay lines in Figure 1(b); the bottom path's (D-1)/2 length delay, and the comb subfilter's delay line of length D within the MA.

This DC removal network's passband performance, when D = 31, is shown in Figure 1(c). (The frequency axis value of 0.5 corresponds to a cyclic frequency of half the input signal's Fs sample rate.) While the network has the desired infinite attenuation at zero Hz, its passband peak-peak ripple is unpleasantly large at 2.9 dB. We can do better, as we shall see.

If D is an integer power of two the 1/D scaling in (1) can be performed using a binary right shift by log2(D) bits, making Figure 1(b) a multiplier-free network. However, in that scenario the MA's group delay is not an integer number of samples, making it difficult to synchronize the delayed x[n] and the v[n] sequences. To solve this problem we can use two cascaded D-point MAs as shown in Figure 2(a). Because the cascaded MAs have an integer group delay of D-1 samples, we can be clever and tap off the first moving averager's comb delay line eliminating the bottom-path delay line in Figure 1(b). This way we still only need implement two delay lines in Figure 2(a), one delay line for the comb subfilter in each MA.

DC Removal Figure 2

Figure 2.

The magnitude response of the Figure 2(a) dual-MA DC removal network, for D = 32, is shown in Figure 2(b). In that figure we show the DC removal filter's passband with its narrower transition region width and a much improved peak-peak ripple of 0.42 dB. What we've created then is a linear-phase, multiplierless, DC removal network having a narrow transition region near zero Hz.

Happily, it's worth noting that standard tapped delay-line, linear phase, highpass FIR filter designs using least-squares error minimization, or the Parks-McClellan method, require more than one hundred taps to approximate our D = 32 DC removal filter's performance.

On a practical note, the MAs in Figure 2(a) contain integrators that can experience data overflow. (An integrator's gain is infinite at DC.) Using two's complement fixed-point arithmetic avoids integrator overflow errors if we ensure that the integrator (accumulator) bit width is at least:

DC Removal Equation 2

where q[n] is the input sequence to an accumulator, and ékù means: if k is not an integer, round it up to the next larger integer.

For an even narrower transition region width, in the vicinity of zero Hz, than that shown in Figure 2(b), we can set D to a larger integer power of two, however this will not reduce the DC removal network's passband ripple.

At the expense of three additional delay lines, and four new addition operations per output sample, we can implement the linear phase DC removal filter shown in Figure 3(a). That quad-MA implementation, having a group delay of 2D-2 samples, yields an improved passband peak-peak ripple of only 0.02 dB, as shown in Figure 3(b), as well as a reduced-width transition region relative to the dual-MA implementation.

DC Removal Figure 3

Figure 3.

The frequency magnitude responses of the Figure 3(a) network, for several values of D, are shown in Figure 4.

Figure 4.

So the bottom line to this blog is: At the expense of multiple delay lines, it is possible to perform efficient (multiplier-free) linear-phase DC removal.

[1] Posts during July 31 - August 8, 2007, Subject: "Regarding DC removal filter".



[ - ]
Comment by lishenghu365April 14, 2008
Thanks ,it is a terrific job!!!
[ - ]
Comment by Ahmad1990April 30, 2018
Rick Lyons

can you please send me a matlab code of above article 

[ - ]
Comment by Rick LyonsApril 30, 2018

Ahmad1990, here is some MATLAB code for you to work with:

% Filename: DC_Removal_for_distribution.m
%
%    Models the " DC Removal" scheme using cascaded
%    recursive running sum (RRS) filters and subtracting  
%    the filter's output from a delayed version of the input.
%    Models single-stage, 2-stage, and 4-stage RRS filters.
%     Richard Lyons, April, 2018

clear, clc

% Set the desired delay line length
D = 31 % WARNING!! For a 1-stage DC blocker 'D' must be an odd number.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Design 1-, 2-, & 4-stage lowpass RRS (moving averaging) filters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_1Stage = (1/D)*[1, zeros(1,D-1), -1]; A_1Stage = [1, -1];
B_2Stage = conv(B_1Stage, B_1Stage); A_2Stage = conv(A_1Stage, A_1Stage);
B_4Stage = conv(B_2Stage, B_2Stage); A_4Stage = conv(A_2Stage, A_2Stage);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Add or delete the appropriate percent signs (%) in front of
%  'B_filter' variables to choose the number of deired RSS
%  stages (single-, two-, or four-)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implement single-stage RSS
%B_filter = B_1Stage; A_filter = A_1Stage; B_delay_line = [zeros(1, (D-1)/2), 1]'; Order = 1;

% Implement two-stage RSS
B_filter = B_2Stage; A_filter = A_2Stage; B_delay_line = [zeros(1, D-1), 1]'; Order = 2;

% Implement four-stage RSS
%B_filter = B_4Stage; A_filter = A_4Stage; B_delay_line = [zeros(1, 2*D-2), 1]'; Order = 4;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analyze freq-domain performance of the DC-removal filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ImpRespFilter,T] = impz(B_filter, A_filter, 4*D+10);
[ImpRespDelay,T] = impz(B_delay_line, 1, 4*D+10);
ImpResp_DC_Remove_Filter = ImpRespDelay -ImpRespFilter;

H_Spec_mag = abs(fft(ImpResp_DC_Remove_Filter, 1024));
H_Spec_mag_dB = 20*log10(H_Spec_mag/max(H_Spec_mag));
%H_Spec_mag_dB = thresh(H_Spec_mag_dB, Threshold); % For plotting
FreqAxis = (0:1023)/(1024);
    
    figure(1), clf
    subplot(3,1,1), plot(FreqAxis, H_Spec_mag,'-r')
    axis([0, 0.5, 0, 1.1])
    title('FFT mag. of cascaded imp. responses')
    ylabel('Linear'), grid on, zoom on
    subplot(3,1,2), plot(FreqAxis, H_Spec_mag_dB,'-r')
    axis([0, 0.5, -30, 2])
    ylabel('dB'), grid on, zoom on
    subplot(3,1,3), plot(FreqAxis, H_Spec_mag_dB,'-r')
    axis([0, 0.2, -0.5, 0.1])
    ylabel('dB'), xlabel('Freq (times Fs)'), grid on, zoom on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute & plot group delay to show phase linearity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Turn off warning when plotting group delay
warning off MATLAB:divideByZero,
[GD,W] = grpdelay(B_filter, A_filter, 1024, 'whole');

    figure(2), clf
    plot(W, GD,'-r', 'linewidth', 1)
    ylabel('Group delay (samples)'), xlabel('Freq in rad./sample')
    grid on, zoom on

[ - ]
Comment by Ahmad1990May 1, 2018

Thanks a lot Mr.Rick Lyons.

[ - ]
Comment by DKDiveDudeDecember 30, 2023

If anyone have the Matlab C Converter, I would very much appreciate if they would post the C code here, thanks!

[ - ]
Comment by an2orApril 2, 2008
Rick As always, an excellent and informative article! Regards, Andor
[ - ]
Comment by bluewhaleJuly 22, 2008
Thats whats expected from an ace.
[ - ]
Comment by xie.qiangApril 8, 2009
Hi Rick, It's perfect. Since there is attenuation at around the very low frequency for the DC revomal filter. Is there any way to remove the low frequency DC by using this filter?
[ - ]
Comment by Rick LyonsApril 9, 2009
Hello xie.qiang, Sorry, I'm not sure what you mean by "low frequency DC". The filter has infinite attenuation at DC (zero Hz). Regards, [-Rick-]
[ - ]
Comment by xie.qiangApril 14, 2009
Hi Rick, I means that when the signal frequency is very low compare to the sampling rate, say 20 Hz with the sample rate of 192KHz, and the DC removal filter will also attenuate the 20Hz signal but what I really want to do is need to removal the DC in the 20Hz signal. Hope this could be clear. :-)
[ - ]
Comment by Rick LyonsApril 14, 2009
Hello xie.qiang, I am still having trouble understanding you. What do you mean by "the DC in the 20Hz signal"? Do you mean that you have a signal that is a 20 Hz sinewave riding on a DC bias (a 20 Hz sinewave added to a constant)? Regards, [-Rick-]
[ - ]
Comment by xie.qiangApril 15, 2009
Hi Rick, Sorry for the mis-understanding. Yes, it's correct, I have a signal riding on a DC bias, and my sampling rate is 192K, by using the DC remvoal filter, it will do remove the DC bias, but it will also attenuate the 20 Hz Signal, is it? Actually, it is so interested to me is that I have such an application that I want to get the DC bias of my 20 Hz signal. For example, I want to calculate the RMS value of 10Hz signal of x(n), first multiplied x(n) by it's self get x2(n), [x2(n) = x(n) * x(n) ], and then x2(n) is a 20Hz signal ride on a DC bias, the DC is the RMS value. One solution is that I can apply a IIR filter to the signal of x2(n), however, it will wait so long time to get the stabe output of IIR. Since if I can use your DC removal filter, then I can remove the DC bias of signal x2(n), by using the output of DC removal filter y2(n) minus the x2(n), I will finally get the DC (or RMS of 10Hz signal).. Do you think is it possible?? :-)
[ - ]
Comment by Rick LyonsApril 15, 2009
Hello xie.qiang, before you could use the above DC removal filter you would have to decimate your 192 kHz signal by 400 or 500. (Your 20 Hz signal is far too low in frequency, relative to Fs sample rate, to use the DC removal filter.) You might consider using a high-order decimating CIC filter to extract your DC signal from x2(n). Just a thought.) Good Luck, [-Rick-]
[ - ]
Comment by xie.qiangApril 15, 2009
Hi Rick, Thanks your suggestion. Actually, my singnal to be measured is from 10Hz to 48KHz, it is ture only that the low frequency bother me so much.
[ - ]
Comment by gopalJSeptember 24, 2009
hii Rick and all, I want to remove DC and very low frequency components from my signal. The sampling frequency is 40 Hz. I want to remove DC and frequencies upto 0.3Hz. The transition band should be less than 0.2 Hz. Please help me on this. The delay should not be more than 15 samples. I would be really grateful for this. Thanking all of you in advance.
[ - ]
Comment by Rick LyonsJanuary 12, 2010
To gopalJ, For some reason I did not see your comment until today. Do you still need help? [Rick Lyons-]
[ - ]
Comment by AfinkoJanuary 12, 2010
Hi Rick, it is a great solution with very less computation requirement. I want to implement it in DSP. Is there a way how to change the "cut-off" frequency of this filter? Thanks. Afi
[ - ]
Comment by Rick LyonsJanuary 12, 2010
To Afi, If I understand your question, then my answer is, "No, not in the usual sense." The cutoff freq (the beginning of the highpass filter's passband) depends on how many moving averager (MA) operations exist in the top path of the filter. See how Figure 3's cutoff freq is a little lower than the Figure 2 filter's cutoff freq. The cutoff freqs of those two filters are fixed (uncontrollable). Good Luck, [-Rick Lyons-]
[ - ]
Comment by AfinkoJanuary 13, 2010
Thnaks Rick, You understand my question correct. Do you have any idea how to modify the schematic in order to obtain cutoff freq at approx. 0.1 Fs? The passband ripple should be less then 0.2dB. The FIR filter with such steep attenuation require huge number of coefficients...
[ - ]
Comment by Rick LyonsJanuary 13, 2010
To Afi, No, I know of no way to change the DC-removal filter's structure (block diagram) to accommodate your design parameters. However, your two design parameters are insufficient to specify a highpass filter. To completely specify a highpass filter you must also specify the filter's transition region width and stopband attenuation (which you did not do). I played around, a little, with MATLAB and learned that you can build your highpass filter having a transition region width of 25 Hz and a stopband attenuation of 63 dB using a 100-tap tapped-delay line FIR filter. If you use the "folded" structure of such a filter, you can reduce it's computational workload to roughly 50 multiplies per filter output sample. That's not too bad, is it? [-Rick-]
[ - ]
Comment by AfinkoJanuary 13, 2010
Thanks Rick, Now I am using FIR filter with response shoen in figure: http://img705.imageshack.us/img705/4079/97967298.jpg This filter I generated by MATLAB FDAtool and it has 500 coefficients. Is there any another way how to obtain filter like this with small computation requirements? Afi
[ - ]
Comment by Rick LyonsJanuary 13, 2010
To Afi, Please send me an E-mail at: [-Rick-]
[ - ]
Comment by Leo2011February 15, 2011
Hi Rick, I've a signal from 10Hz to 100kHz in the band and sample with 300kHz, and would like to remove the DC from those signal in band. May I know is there a method to do so? I did try the DC bias removal method from streamlining DSP, but the roll-off in passband is very huge.
[ - ]
Comment by igorlkOctober 14, 2011
As far as I know, recursive moving average converges to gaussian filter.
[ - ]
Comment by aditijoisherNovember 23, 2020

Hey I am doing a project on dc blocker and i need help!!

Project states:

Demonstrate any filter (FIR/IIR) with all parameter of filter .

DC blocker:

Steps:

  1. Choose a suitable Transfer Function suitable for the application
  2. Pass signal (x(n))
  3. Find output discrete /analog
  4. State which kind of filter (eg.FIR/IIR)
  5. Find out stability,PZ,side lobes ,transition band,order of the filter,attenuation

Please help!

[ - ]
Comment by Rick LyonsNovember 24, 2020

Hello aditijoisher. I suggest your go to the following web page:

https://www.dsprelated.com/freebooks/filters/DC_Blocker.html#:~:text=The%20dc%20blocker%20is%20a%20small%20recursive%20filter,%28%29%20and%20a%20pole%20near%20dc%20at%20.

Implement the IIR filter described by the web page's Equation (B.10).

Analyze that IIR filter to answer all the questions in your above "Step# 5.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Please login (on the right) if you already have an account on this platform.

Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: