# Removing transients in highpass filtering with MATLAB

I've made a simple first order IIR highpass filter with a zero at z = 1 and a pole at z = 0.9. Its frequency response looks like this:

Now, I filter a DC signal using this filter. Here's the MATLAB code I use to do it:

b = [1 -1]; % Zero at z = 1
a = [1 -0.9]; %Pole at z = 0.9

figure(1)
freqz(b, a)

t = 1:100;
x(1:length(t)) = 1; % Constant function

y = filter(b, a, x);

figure(2)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(3)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');


As my filter is highpass, I expect the DC to become zero, or atleast become severely attenuated. However, the output signal I get looks like this:

From my understanding, this exponential output is a transient produced because I haven't set the initial conditions correctly. Sure enough, setting x[-1] = 1 solves the problem. However, this works only for this particular input DC signal. For any general input signal, how do I set the initial conditions so that transients aren't produced?

Edit : I'm aware that the filtfilt() function does forward-backward filtering with transient minimization, but I really want to port the filter to an embedded platform, so I need to understand how transient removal works. Thanks in advance for the help!

Edit 2 : As suggested by Kuba Ober, I tried setting x[-1] as the value that it actually should have been. It works fine for a DC input, but here's what happened for a sinusoidal input:

clc; clear all;
p = 0.9

a = [1 -p]
b = [1 -1]

n = 1:100; % Samples
f = 0.2; % Frequency in Hz
Fs = 10; % Sampling rate in samples per second
t = n/Fs; % Time axis

x = sin(2*pi*f*t);

% Filter with the appropriate initial conditions
y = filter(b, a, x, filtic(b, a, [], [sin(2*pi*f*0)]));

figure(1)
plot(t, x)
xlabel('Time');
ylabel('Input Signal');

figure(2)
plot(t, y)
xlabel('Time');
ylabel('Output Signal');


Here's the input signal :

And here's the output :

The first peak is visibly smaller than the second, which indicates some transients being present. I'm not entirely sure about this, but I think the reason it doesn't work is because just setting x[-1] is not enough, I also need to set y[-1]. The problem here is that there's no way to find out what y[-1] actually should be.

Edit 3 : Let me provide a little more info on the problem I'm working on. I'm trying to use filters to remove noise from ECG (Electrocardiogram) signals in an embedded platform. Here's a typical ECG signal, after filtering:

Here's what an ECG signal looks before filtering:

Note the DC offset in the signal before filtering. For filtering, I need a notch filter to remove high frequency power line noise and a highpass filter to remove the DC and the low frequency "drifting" of the signal.

The filters I use need to be linear phase, since the time domain morphology of an ECG signal is very important for diagnosis. However, my filter doesn't need to be real-time, as I'm doing the processing offline after acquiring the ECG signal from the patient. So, for implementing nonlinear phase IIR filters, I'm currently using forward-backward zero phase filtering.

One opinion that's shared by @Matt L. and @Drazick is that transients are unavoidable in real-time filtering and that I should use a longer input signal and crop off the first few seconds of the filtered output instead. This is something I'd like to avoid, as acquiring long ECGs from a living patient is somewhat difficult. Also, I do not have to filter in real-time, so any technique of transient removal that hinges on knowing the entire signal in advance is perfectly all right. Any help is appreciated!

• Can't you "set" x[-1] to whatever value it actually had? As in - take one more sample before the samples that you'll filter, and use its value? – Reinstate Monica Jun 16 '14 at 16:42
• Hi, isn't it that you need to set x[-1] = -1? Actually, filter command in MATLAB has provision to set the initial value. – learner Jun 17 '14 at 1:38
• In my opinion if the focus is to remove DC then initial conditions must be set to -1. In that case there might be a penalty if input signal has no DC. – learner Jun 17 '14 at 1:46
• @Kuba Ober : That's a great idea! Unfortunately, it doesn't quite seem to work. I've tried out your suggestion and detailed the results in an edit. – Azura Jun 17 '14 at 5:38
• @learner : I assume that you mean that I should call the filter function using filter(b, a, x, -1). That is absolutely correct, but it's not the same as setting x[-1] = -1. The fourth argument is actually the initial conditions for the delays in the direct-form II implementation. To actually set x[-1] = -1, you need to use filter(b, a, x, filtic(b, a, [], [-1]). This doesn't quite work, though - only setting x[-1] = 1 works for removing a DC with amplitude 1. link – Azura Jun 17 '14 at 5:46

Your first order filter recursion for some real constants $a,b,c$ is

$$y[n] = a x[n] + b x[n-1] - c y[n-1]$$

with the two initial memory states $x[-1]$ and $y[-1]$ at $n=0$.

Your "no transient" condition can be translated to $y[0]=0$ and a necessary second condition so that you can solve for both of your memory states. The second condition could be, that the discrete derivative of $y$ also vanishes at $n=0$, so $y[0]-y[-1]=0$. You can also take any other condition that seems sensible to you.

The two equations give you a unique solution for the two unknown memories, namely: $$y[-1] = 0$$ and $$x[-1] = -\frac{a}{b}x[0]$$

Alternatively, your conditions may better be chosen as $y[0]=0$ and $y[0]-y[-1]=x[0]-x[-1]$ in order to capture the initial slope of the input. The resulting recursion equation at $n=0$ is then $$0=a x[0]+b x[-1]+c(x[0]-x[-1])$$ giving you the solution $$x[-1]=-\frac{a+c}{b-c}x[0]$$ and $$y[-1]= -\left(1+\frac{a+c}{b-c}\right)x[0]$$

• @Drazick, I explicitly wrote in my answer that setting initial conditions like these is not a general solution. Of course this is making assumptions, but I disagree with your understanding of "work". This always works, it might just not give the result you expected. The long term output signal will always be zero mean because the filter of the OP rejects $\omega=0$. Only the transient response will interfere with that, and that is why enforcing this kind of response helps reducing the transient response. Also, the second method does not even give $y[-1]=0$ for the general case. – Jazzmaniac Jun 20 '14 at 12:24