# Linear-phase DC Removal Filter

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

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

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 *F*_{s} 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 log_{2}(*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.

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:

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 2*D*-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.

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

**Previous post by Rick Lyons:**

Free DSP Books on the Internet

**Next post by Rick Lyons:**

Computing Large DFTs Using Small FFTs

- Comments
- Write a Comment Select to add a comment

can you please send me a matlab code of above article

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

Thanks a lot Mr.Rick Lyons.

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.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.