CODEWITHSUNDEEP
matlabsignal-processingconvolution
yurnero
yurnero

Reputation: 315

MATLAB's filter: how to setup the initial conditions

I have a very simple problem: generate y_1, ..., y_n such that

y_t = x_t + a * y_{t-1} + b * y_{t-2}

where a, b, x_1,...,x_n are given, as are y_0 and y_{-1}. I can do this with MATLAB in 2 ways (with a loop and with some manual adjustments before using filter). How do I do this using the built-in allowance for initial conditions of filter?

My code is as follows:

clc; clear; tic;

n        = 20;
x        = (1:1:n)';  %//' this is here so code looks good in SO.com
yzero    = -3;
yminus1  = -4;
a        = 1;
b        = 1;
yinitial = [yminus1;yzero];

% alternative 1
y = zeros(n,1);
for i = 1:n
    if i == 1
        y(1) = x(1) + a * yzero + b * yminus1;
    elseif i == 2
        y(2) = x(2) + a * y(1) + b * yzero;
    else
        y(i) = x(i) + a * y(i-1) + b * y(i-2);
    end
end

% alternative 2
z    = x;
z(1) = z(1) + a * yzero + b * yminus1;
z(2) = z(2)             + b * yzero;
s    = filter(1,[1,-a,-b],z);

% alternative 3
r    = filter(1,[1,-a,-b],x,yinitial);


toc;

I was hoping to get the same answers for y, s, and r. But only y and s are the same and r is very different. What am I doing wrong please?

Upvotes: 0

Views: 502

Answers (2)

yurnero
yurnero

Reputation: 315

For alternative 3, using yinitial as I have done is incorrect. Here are 2 ways to fix it. First, I need to redefine yinitial:

yinitial = [yzero;yminus1];

Then, I get xinitial either from

xinitial = filter([b,a],1,yinitial);
xinitial = xinitial(end:-1:1);

or from

xinitial = filtic(1,[1,-a,-b],yinitial);

Then, I call

r        = filter(1,[1,-a,-b],x,xinitial);

Upvotes: 0

dfrib
dfrib

Reputation: 73196

The first element of the initial conditions are slightly off in the third alternative.

You use initial conditions:

yinitial =

    -4
    -3

But the correct ones are

yinitialFix =

    -7
    -3

Why? the 4th argument in call to filter in alternative 3 contains the initial conditions: in so that these initial steps will be performed:

x(1) = x(1) + yinitialFix(1)
x(2) = x(2) + yinitialFix(2)

Hence, you've actually (implicitly) computed the correct values to yinitialFix yourself in alternative 2, in your initial step:

  • initial step shifts z(1) by a * yzero + b * yminus1 = -7,
  • in the same way, z(2) is shifted by b * yzero = -3.

With the applied fix, alternative 3

%// ...
%// alternative 3;
yinitialFix = [-7; -3];
r    = filter(1,[1,-a,-b],x,yinitialFix);

yields the same result as the two previous ones.

Upvotes: 1

Related Questions