2
$\begingroup$

I need to smooth some data segments and one of the techniques I am experimenting is a low-pass filter. So I designed a low-pass filter in matlab and applied it to my signals. Here is what happened with two of them.

enter image description here

Panel A is one of the data segments plotted with its filtered signal. It looks good. This is what I expect and Panel B is its signed residual which is the filtered signal minus the original signal. I see there is no DC offset and no linear trend so I am happy.

Panel C is another data segment plotted with its filtered result and there is an obvious offset between the two. The filtered signal "doesn't run through" the signal like how I think it should. Panel D is its residual and we confirm that the mean is not zero. And sometimes (not shown here) I see a linear trend as well in the residual.

So my question is what is going on here? How should one interpret this? I am using low-pass filtering so the DC component should be included so the mean of the residual should be close to zero. Does this say something about the power at high frequencies like is there "too much" power at higher frequencies? Something with leakage or something like power leaking from the zero channel to other neighboring frequencies? Maybe I did something wrong? Is there anyway to fix this? I need to be sure that all of the residuals are like Panel B with no linear trends and the mean being around zero.

The bigger focus is to look at both the smoothed signal (the low-pass filtered signal) AND the residual with the emphasis being on the residual. If this can be fixed then good otherwise I am thinking after filtering like this, then subtracting a least squares fit straight line from the residual. Then adding the signed residual to the original signal to get the smoothed-version that goes "through" the signal.

Yep, I know guys that generally it isn't a good idea to low-pass and subtract. But long story and I have no choice. As you can see in the code, the data segments all have length 41 with 30 seconds cadence so the frequency step size is 0.000813 which is of course also the first bin after the zero channel. So what I am trying to do is to get rid of the power in the zero channel and then look at all other (higher) frequency bins so in other words, detrending (hence my cutoff frequency at 9mHz). The problem is that I need a smoothed version of the signal as well as the detrended signal. So as far as I can tell, either low pass and then subtract...OR...high pass and then add.

Thanks everyone.

Here is my code


% Filter Design
Fpass = 1e-09;   % Passband Frequency
Fstop = 0.009;   % Stopband Frequency
Apass = 1;       % Passband Ripple (dB)
Astop = 60;      % Stopband Attenuation (dB)
Fs    = 1/30;    % Sampling Frequency
h = fdesign.lowpass('fp,fst,ap,ast', Fpass, Fstop, Apass, Astop, Fs);
Hd = design(h, 'equiripple', ...
        'MinOrder', 'any', ...
        'StopbandShape', 'flat');
xavg = filtfilt(Hd.numerator,1,x);
yavg = filtfilt(Hd.numerator,1,y);
xdelta = xavg - x;
ydelta = yavg - y;

t=1:41;

figure
subplot(2,2,1)
plot(t,x,'o',t,xavg)
legend('x','Filtered x')
title('Panel A - Signal x and Filtered Signal x')
subplot(2,2,3)
plot(t,y,'o',t,yavg)
legend('y','Filtered y')
title('Panel C - Signal y and Filtered Signal y')

subplot(2,2,2)
plot(t,xdelta)
title('Panel B - Signed Residual - Filtered x minus x')
subplot(2,2,4)
plot(t,ydelta)
title('Panel B - Signed Residual - Filtered y minus y')

Edit: I didn't see the comments before. I just saw them right now and checked out what geometrikal pointed out. And here is the pic enter image description here different data sets but same idea. I couldn't find the original numbers I used so just took another set of numbers. So here I am only plotting two signals, no residuals. Call the top X which gives a bigger offset and the bottom Y. So for both of the signals, I plot the signals with results of the various filtering methods. Here is the code to show what I am doing.


% Just normal filtering
xfilt = filtfilt(Hd.numerator,1,x);
yfilt = filtfilt(Hd.numerator,1,y);

% Filtering by subtracting the mean and then adding it
xfiltmean = filtfilt(Hd.numerator,1,x-mean(x))+mean(x);
yfiltmean = filtfilt(Hd.numerator,1,y-mean(y))+mean(y);

% The scaled filter
sc = mean(filtfilt(Hd.numerator,1,ones(1,41)));
xfiltsc = xfilt/sc;
yfiltsc = yfilt/sc;

and yep the filter doesn't sum to one as the picture clearly shows. sum(Hd.numerator)=0.9546 so which one of the two is the problem is this? What does this mean? Is it the zeropadding or this scaling thing? Because applying both of these fixes individually seems to work or at least they give roughly the same answers. Both of them together like this


(filtfilt(Hd.numerator,1,x-mean(x))+mean(x))/sc;

won't work because now the filtered signal will clearly be too large (and verified this with my data). And the magnitude of the data doesn't seem to matter even if I do


x = x+1000;

I get the same results. Both filters are almost identical. So which fix should I apply? Removing the mean or scaling?

Thanks guys for all your help.

$\endgroup$
  • 4
    $\begingroup$ A low-pass filter should be expected to pass DC or zero frequency, no? $\endgroup$ – Dilip Sarwate May 4 '13 at 1:39
  • $\begingroup$ @DilipSarwate: but the OP question is that the low-pass filter is not passing DC correctly. This is really a high-pass filter formed by subtracting the low-pass filter from the data. (i.e., the title is wrong.) $\endgroup$ – Wandering Logic May 7 '13 at 12:05
  • $\begingroup$ Set one of the first few samples to 1000 and see what happens, i.e. x(1) = 1000. $\endgroup$ – geometrikal May 8 '13 at 7:59
  • $\begingroup$ The fix you should apply is to divide the filter by its sum. (Normalize the filter to 1). @geometrikal was right! You should switch the right answer to geometrikal's post. $\endgroup$ – Wandering Logic May 8 '13 at 11:23
  • 1
    $\begingroup$ Done, indeed @geometrikal was right. I tried by making some outliers and indeed the scaling works much better than the removing/adding mean so that was the underlying problem. Thanks everyone! $\endgroup$ – Fixed Point May 8 '13 at 23:01
2
$\begingroup$

It could be a few things.

Zero padding of the input signal:

Some filtering operations pad the end of the signal with zeros before convolving it with the filter kernel. I don't think filtfilt does this though.

Filter doesn't sum to 1:

Lets say you had a discrete signal that was all 1's. The low pass filtering should also return all 1's. Because you are operating on discrete signals, even if the filter is designed properly the errors that come in through discretisation can cause the filter to give a different value. For example, compare the discrete Gaussian kernel to the sampled Gaussian kernel below (image taken from wikipedia, more info at http://en.wikipedia.org/wiki/Sampled_Gaussian_kernel). I think this is the most likely explanation. You could scale your filter output as a crude fix.

enter image description here

The Gaussian kernel is often used for low-pass smoothing because it generates a scale-space. In practice this means that smoothing does not generate any new features in the signal, and that more smoothed signals will always have the same or less detail as less smoothed signals. I would try that if you haven't already.

$\endgroup$
  • $\begingroup$ Please see my edit. Thanks! Also may I ask why don't you think that filtfilt zeropads? $\endgroup$ – Fixed Point May 8 '13 at 5:57
  • $\begingroup$ @FixedPoint Because there doesn't appear to be edge effects. Imagine if you added 10 zeros on each side of your data, the filter would start heading towards 0 a few samples each side of the edge. $\endgroup$ – geometrikal May 8 '13 at 8:02
  • $\begingroup$ Makes sense, but what if those edge effects are being deleted. I am thinking like zeropad, filter, then remove the edges. Could MATLAB be doing something like that? $\endgroup$ – Fixed Point May 8 '13 at 22:40
3
$\begingroup$

The DC offset is just the average of all your samples. You should calculate that first, subtract it, then run your filter.

I think what is happening here is that your data is getting padded with zeros on the left and right, so you are actually averaging your real data with a bunch of extra zeros. (I don't know matlab, so I may be very wrong about this.) So the DC value you are getting right now is way off (low because your data is all large and positive.)

$\endgroup$
  • $\begingroup$ Oh my good god, talk about Occam's razor. So obvious when you point it out and you are absolutely right. I checked this with positive and negative signals. The bias is always towards zero and it is larger for larger values. So now I subtract the mean, filter, and then add the mean back to get a smooth signal. And of course there's no bias anywhere now...at least none that I can see. $\endgroup$ – Fixed Point May 4 '13 at 0:01
  • $\begingroup$ Another unrelated question, I was also playing with high pass filters. Is the order of the filter just length(Hp.numerator)? Because MATLAB keeps telling me that the length of the data segment should at least be 3 times as the filter order. In my case the data length is 41 and length(Hp.numerator)=81 with Hp being my high-pass filter. $\endgroup$ – Fixed Point May 4 '13 at 0:03
  • $\begingroup$ I think this assumes the DC level is centred in the data. E.g. what happens when there is an outlier? Set one of your data points to something really big and that will give you a big offset! $\endgroup$ – geometrikal May 4 '13 at 11:59
  • $\begingroup$ For a crude fix: make a signal of all ones x = ones(1,41), apply the filter, xavg will be all the same of another value, lets call it sc. Now apply the filter to a real signal and multiply the result by 1/sc. $\endgroup$ – geometrikal May 4 '13 at 12:05
  • $\begingroup$ @geometrikal: I was using "DC offset" in the same sense that the OP used it: as the 0-frequency component of the input signal. That's the mean. If you subtract the mean from the input data, then the mean of the result will be 0. Outliers don't change the definition of the mean. $\endgroup$ – Wandering Logic May 5 '13 at 3:05

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.