remove DC component of real time signal

I am using a interface for acquiring a signla in matlab.

The signal is sampled at 500Hz. I have to remove the DC component in real time.

How can I do that?

I know that simulink has the DC blocker block, but I need to accomplish this inside a matlab routine!

Answer 1

[Star Strider]

The D-C component is the mean of the signal, so for a recorded signal, simply subtract the mean.

Otherwise, since the D-C offsset is defined as having a frequency of zero (rad/s, Hz, or anything else), a highpass filter (or bandpass filter) with its passband (or lower passband) slightly above zero should remove the D-C offset. The resulting signal should then have a zero D-C offset.

[L]

@Star Strider My signal is not recorded, so subtracting the mean is not a good alternative. It is a real time signal coming in . I can record up to only 30 samples. I tried to subtract the mean from that, but the signal gets higly distorted and 30 is a very low number of samples to do a high pass filter.

[Star Strider]

Not really. An elliptic filter can be of order 6 or 8 and still be quite effective, if implemented as a second-order-section (‘sos’) realisation, even using filtfilt. I can design it here if you provide the sampling frequency.

Otherwise, implementing it in hardware as a Bessel filter (combing anti-aliasing and D-C offset removal) can be relatively straightforward if you have the power supply available. (I’ve not done any hardware filtering in a while, however there should be references for it online, if not already implemented and requiring only some ‘outboard’ components.)

[L]

sampling is 500Hz

My signal of interest is actually between 0.25 and 4 hz.

So I could band pass it.. but if I could just remove the DC will also help.

[L]

@Star Strider

sampling is 500Hz (2ms). Matlab receives a vector of 10 points at each 20ms. I can save up to 30 samples, if I want.

My signal of interest is actually between 0.25 and 4 hz.

So I could band pass it.. but if I could just remove the DC will also help.

What is important is to keep the signal without any introduced delays.

[Star Strider]

This is a difficult filter to realise because the lower transition region is too steep. I can make it work with a lower attenuation (35 dB rather than 50 dB) for ‘Rs’ and then it distorts the filtered signal significantly less. I encourage you to experiment with various passband, stopband, and atttenuation values to get the result you want.

Fs = 500;                                                           % Sampling Frequency (Hz)
Fn = Fs/2;                                                          % Nyquist Frequency (Hz)
Wp = [0.25 4.0]/Fn;                                                 % Passband Frequency (Normalised)
Ws = [0.01 5.0]/Fn;                                                 % Stopband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 35;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 

[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^16, Fs)                                                % Filter Bode Plot
set(subplot(2,1,1), 'XLim',[0 10])                                  % Zoom To Region-Of-Interest
set(subplot(2,1,2), 'XLim',[0 10])

L = 10*30;
t = linspace(0, Fs*L, Fs*L+1).'/Fs;
signal = sum(sin(2*pi*t.*[0.1:0.1:5 10:10:240]),2);

filtered_signal = filtfilt(sos,g,signal);

figure
plot(t, signal, 'DisplayName','Original')
hold on
plot(t, filtered_signal, 'DisplayName','Filtered')
hold off
xlabel('Time (s)')
ylabel('Amplitude')
legend('Location','best')

[FTsignal,Fv] = FFT1(signal, t);
[FTfiltered_signal,Fv] = FFT1(filtered_signal,t);

figure
plot(Fv, abs(FTsignal)*2, 'DisplayName','Original')
hold on
plot(Fv, abs(FTfiltered_signal)*2, 'DisplayName','Filtered')
hold off
grid
xlim([0 5])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('Location','best')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

.