関数filterと自作のプログラムの不一致

Question

MATLABの関数にあるconvを使わないで以下の式を実現したくてプログラムを作成しているのですが出力結果がうまくいかず解決方法に悩んでいます。教えていただけないでしょうか。1つ前にも同様の質問をしていますが、プログラムを修正しています。
具体的にはフィルタ係数a,bを指定する際にa,bともに1より大きな値で指定するとfilter(b,a,x)でプロットした結果と一致するのですが、フィルタ係数を1より小さい値にしたときに一致しません。

prompt = 'What is the value of "a" ? ';
a = input(prompt);
prompt = 'What is the value of "b" ? ';
b = input(prompt);

N=49;
x = [1; zeros(N,1)];


    
b_length=length(b);
a_length=length(a);
k1=length(x)+length(b)-1;
y1=zeros(1,N);
y2=zeros(1,N);

x1=length(x);
y=zeros(1,N);


 for n=1:k1
 for k=max(1,n+1-x1):min(n,b_length)  
      
      y1(n)=y1(n)+x(n-k+1)*b(k);
end


%for m=1:a_length
for m=max(1,n+1-x1):min(n,a_length)  
      y2(n)=y2(n)+y(n-m+1)*(-a(m));
        

end
y(n)=y1(n)+y2(n);
end
plot(y)
xlim([0 50])

左が自作のプログラムの結果で右がfilter関数を使った場合です。

Comments

[jbpb0]

コードを詳しく読んでませんが、とりあえず気づいたことだけ書きます

[その1]
https://jp.mathworks.com/help/matlab/ref/filter.html#buagwwg-2
の「詳細」の「有理伝達関数」の、
「有理伝達関数は、次の差分方程式で表すこともできます。」
と書かれてるところのすぐ下の式を見てください
「a(1)」は左辺の係数であり、右辺で使われているのは「a(2)」以降です

質問に記載の式に合わせて「y(n)=...」とするには、右辺全体を「a(1)」で割ります
a(1)を使うのはそこだけで、a(2)以降の使い方の個別のy()との掛け算には使いません

a(1)の扱いが異なるので、質問者さんが書いたコードと
filter(b,a,x)
の結果は同じにはならないと思います

for m=max(1,n+1-x1):min(n,a_length)
を、a(1)を使わないように
for m=max(2,n+2-x1):min(n,a_length)
に変えて、
y(n)=y1(n)+y2(n);
を
y(n)=(y1(n)+y2(n))/a(1);
に変えるのではないですかね

[その2]
質問のコードを実行してもエラーになります

y1, y2, yのサイズは
y1=zeros(1,N);
y2=zeros(1,N);
y=zeros(1,N);
なので、どれもN=49です

k1は
k1=length(x)+length(b)-1;
なので、length(x)=50以上ですから、
for n=1:k1
でnは50以上になり、そのループ内での
y1(n)=...
y2(n)=...
y(n)=...
で、サイズが49しかないy1, y2, yに、それよりも大きなnでアクセスしようとして、エラーになります

y1, y2, yのサイズを、
y1=zeros(1,x1);
y2=zeros(1,x1);
y=zeros(1,x1);
でxと同じx1にして、
for n=1:k1
を
for n=1:x1
に変えるのではないですかね

Answer 1

filter
の「詳細」の「有理伝達関数」の、
「有理伝達関数は、次の差分方程式で表すこともできます。」
と書かれてるところのすぐ下の式を見てください
「a(1)」は左辺の係数であり、右辺で使われているのは「a(2)」以降です

質問に記載の式に合わせて「y(n)=...」とするには、右辺全体を「a(1)」で割ります
a(1)を使うのはそこだけで、a(2)以降の使い方の個別のy()との掛け算には使いません

a(1)の扱いが異なるので、質問者さんが書いたコードと
filter(b, a, x)
の結果は同じにはならないと思いますので、そこを直しました

N=49;
x = [1; zeros(N,1)];

abl = randi(N);
a = rand(1, abl);
b = rand(1, abl);

b_length=length(b);
a_length=length(a);
x1=length(x);
y1=zeros(1,x1);
y2=zeros(1,x1);
y=zeros(1,x1);

for n=1:x1
  for k=max(1,n+1-x1):min(n,b_length)
    y1(n)=y1(n)+x(n-k+1)*b(k);
  end

  for m=max(2,n+2-x1):min(n,a_length)
    y2(n)=y2(n)+y(n-m+1)*(-a(m));
  end

  y(n)=(y1(n)+y2(n))/a(1);
end

g = filter(b, a, x);

plot(1:x1, y, 1:x1, g)

err = (y - g') ./ g';
disp(max(err))
disp(min(err))

a, bを乱数で決めているので実行する度に結果が変わりますが、常にグラフがピッタリ重なります
また、相対誤差を数値で表示してますが、とても小さな値(数値計算誤差)にしかなりません

比較のため、質問に記載のコードから、同様にa, bを乱数で決めるように変えたものを、以下に掲載します
実行してみると、ほとんどの場合でグラフが重なりません

N=49;
x = [1; zeros(N,1)];

abl = randi(N)
a = rand(1, abl)
b = rand(1, abl)

b_length=length(b);
a_length=length(a);
k1=length(x)+length(b)-1;
y1=zeros(1,N);
y2=zeros(1,N);
x1=length(x);
y=zeros(1,N);

for n=1:k1
  for k=max(1,n+1-x1):min(n,b_length)
    y1(n)=y1(n)+x(n-k+1)*b(k);
  end

  for m=max(1,n+1-x1):min(n,a_length)
    y2(n)=y2(n)+y(n-m+1)*(-a(m));
  end

  y(n)=y1(n)+y2(n);
end

g = filter(b, a, x);

plot(1:N, y, 1:x1, g)