deconv

多項式 $2{\mathit{x}}^3+7{\mathit{x}}^2+4\mathit{x}+9$ および ${\mathit{x}}^2+1$ の係数をそれぞれ含む 2 つのベクトル y および h を作成します。y から h の逆畳み込みを計算することで、最初の多項式を 2 つ目の多項式で除算します。逆畳み込みにより、多項式 $$2x+7$$ に対応する商の係数と $$2x+2$$ に対応する剰余の係数が得られます。

y = [2 7 4 9];
h = [1 0 1];
[x,r] = deconv(y,h)

x = 1×2

     2     7

r = 1×4

     0     0     2     2

R2023b 以降

ガウス分布の形状の信号 x を作成します。ランダム ノイズから成るインパルス応答 h によってこの信号の畳み込みを行います。

N = 200;
n = 0.1*(1:N);

rng("default")
x = 2*exp(-0.5*((n-10)).^2);
h = 0.1*randn(1,length(x));
y = conv(x,h);

元の信号、インパルス応答、および畳み込みが行われた信号をプロットします。

figure
tiledlayout(3,1)
nexttile
plot(n,x)
title("Original Signal")
nexttile
plot(n,h)
title("Impulse Response")
nexttile
plot(0.1*(1:length(y)),y)
title("Convolved Signal")

次に、既定の多項式の長除法を使用して、インパルス応答 h について信号 y の逆畳み込みを求めます。この手法では、逆畳み込みの計算は不安定であり、結果が急激に増加することがあります。

[x1,r1] = deconv(y,h);
x1(end)

ans = 7.5992e+90

代わりに、数値的に安定した計算を行うために、最小二乗法を使用して逆畳み込みを求めます。

[x2,r2] = deconv(y,h,Method="least-squares");

両方の逆畳み込みの結果をプロットします。ここでは、最小二乗法により、ガウス分布の形状をもつ元の信号が正しく返されています。

figure
tiledlayout(2,1)
nexttile
plot(n,x1)
title("Deconvolved Signal Using ""long-division"" Method")
nexttile
plot(n,x2)
title("Deconvolved Signal Using ""least-squares"" Method")

R2023b 以降

2 つのベクトルを作成します。xin と h の畳み込みについて、xin と同じサイズの中央部分を求めます。畳み込みが行われた信号 y の中央部分の長さは、完全な長さ length(xin)+length(h)-1、つまり 10 ではなく、7 です。

xin = [-1 2 3 -2 0 1 2];
h = [2 4 -1 1];
y = conv(xin,h,"same")

y = 1×7

    15     5    -9     7     6     7    -1

インパルス応答 h について信号 y の最小二乗法による逆畳み込みを求めます。"same" オプションを使用して、畳み込みが行われた信号 y が中央部分のみであることを指定します (y = conv(x,h,"same") + r)。deconv により、丸め誤差内で x に元の信号が復元されることを示します。

[x,r] = deconv(y,h,"same",Method="least-squares")

x = 1×7

   -1.0000    2.0000    3.0000   -2.0000    0.0000    1.0000    2.0000

r = 1×7
10-14 ×

         0    0.0888    0.1776         0         0         0         0

R2023b 以降

それぞれ 2 つの要素をもつ 2 つのベクトルを作成し、"valid" オプションを使用してそれらの畳み込みを行います。このオプションを使用すると、ゼロが加えられたエッジを含まずに計算された畳み込みの部分のみが返されます。この場合、畳み込みが行われた信号 y の要素は 1 つのみです。

xin = [-1 2];
h = [2 5];
y = conv(xin,h,"valid")

y = -1

インパルス応答 h について、畳み込みが行われた信号 y の最小二乗法による逆畳み込みを求めます。"valid" オプションを使用した場合、deconv では常に元の信号が x に返されるとは限りません。そうではなく、norm(x) を最小化する、逆畳み込み問題の解が返されます。

[x,r] = deconv(y,h,"valid",Method="least-squares")

x = 1×2

   -0.1724   -0.0690

r = -3.3307e-16

解を確認するために、h を使用して、計算された信号 x の完全な畳み込みを求めることができます。畳み込みが行われた信号の中央部分は、逆畳み込み問題を定義した元の y と同じです。

yfull = conv(x,h,"full")

yfull = 1×3

   -0.3448   -1.0000   -0.3448

この問題では、deconv は元の信号とは異なる信号を返します。それは、2 変数をもつ 1 つの方程式 ($$-1=5\cdot x(1)+2\cdot x(2)$$) の解を求めているからです。この系は "劣決定" です。つまり、この系では、変数の数が方程式の数より多くなっています。この系では、最小二乗法を使用して残差ノルム (norm(y - conv(x,h,"valid"))) を 0 に最小化すると、無限の解が存在します。このため、deconv では norm(x) を最小化する解も得られます。

次の図は、この劣決定の問題の状況を示しています。青のラインは、方程式 $$x(2)=-1/2-5/2\cdot x(1)$$ の解の数が無限であることを示しています。オレンジ色の丸は、原点から解のラインまでの最小距離を表します。deconv が返す解はこのラインと円の接点と一致します。これは、原点に最も近い解を表します。

R2023b 以降

2 つの信号 x と h を作成し、それらの畳み込みを行います。y の畳み込みが行われた信号にランダム ノイズを追加します。

N = 200;
n = 0.1*(1:N);

rng("default")
x = 2*exp(-0.8*(n - 8).^2) - 4*exp(-2*(n - 10).^2);
h = 2.*exp(-1*(n - 5).^2).*cos(4*n);
y = conv(x,h);
y = y + max(y)*0.05*randn(1,length(y));

元の信号、インパルス応答、および畳み込みが行われた信号をプロットします。

figure
tiledlayout(3,1)
nexttile
plot(n,x)
title("Original Signal")
nexttile
plot(n,h)
title("Impulse Response")
nexttile
plot(0.1*(1:length(y)),y)
title("Convolved Signal with Added Noise")

次に、正則化係数を指定せずに最小二乗法を使用して、インパルス応答 h についてノイズを含む信号 y の逆畳み込みを求めます。既定では、正則化係数は 0 です。

[x1,r1] = deconv(y,h,Method="least-squares");

元の信号と逆畳み込みが行われた信号をプロットします。ここでは、正則化係数が指定されていない関数 deconv では、ノイズを含む信号から元の信号を復元できていません。

figure;
tiledlayout(3,1);
nexttile
plot(n,x)
title("Original Signal")
nexttile
plot(n,x1)
title("Deconvolved Signal Without Regularization");

代わりに、正則化係数 1 を指定して最小二乗法を使用して、h について y の逆畳み込みを求めます。ノイズを含む信号が関係する問題など、悪条件の逆畳み込み問題では、最小二乗解で過適合が生じないように正則化係数を指定できます。

[x2,r2] = deconv(y,h,Method="least-squares",RegularizationFactor=1);

この逆畳み込みが行われた信号をプロットします。ここでは、正則化係数を指定した関数 deconv により、元の信号が復元されています。

nexttile
plot(n,x2)
title("Deconvolved Signal Using Regularization")