7
$\begingroup$

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:

enter image description here

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:

Exponentially decaying output signal

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 :

enter image description here

And here's the output :

enter image description here

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:

enter image description here

Here's what an ECG signal looks before filtering:

enter image description here

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 @Royi 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!

$\endgroup$
9
  • $\begingroup$ 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? $\endgroup$ Jun 16, 2014 at 16:42
  • $\begingroup$ Hi, isn't it that you need to set x[-1] = -1? Actually, filter command in MATLAB has provision to set the initial value. $\endgroup$
    – learner
    Jun 17, 2014 at 1:38
  • $\begingroup$ 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. $\endgroup$
    – learner
    Jun 17, 2014 at 1:46
  • $\begingroup$ @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. $\endgroup$
    – Azura
    Jun 17, 2014 at 5:38
  • $\begingroup$ @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 $\endgroup$
    – Azura
    Jun 17, 2014 at 5:46

1 Answer 1

9
+50
$\begingroup$

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]$$

(Please check my calculations!)

But in general you cannot expect a simple initial condition to give you the same result as knowing the signal history. So you can only take this to a certain point and in general it would probably be better if you just discarded the transient response of your output.

$\endgroup$
3
  • $\begingroup$ Setting y[-1] = 0 as initial condition for this would only work (Partially) for zero mean signals. By the properties of IIR filters you can not tell the real initial condition of the filter unless you have the history of the signal. $\endgroup$
    – Royi
    Jun 20, 2014 at 10:34
  • 2
    $\begingroup$ @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. $\endgroup$
    – Jazzmaniac
    Jun 20, 2014 at 12:24
  • 1
    $\begingroup$ You are right, great answer. $\endgroup$
    – Royi
    Jun 20, 2014 at 12:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.