MATLABで「たたみこみ」を語る

## はじめに
　「たたみこみ」とは「合成積」のことです（と言われるとよけい分からない）。
・和室とは関係ありません
・柔道の寝技でもありません。
東北工業大学の中川先生は、[「人生、たたみこみ」](https://www.ice.tohtech.ac.jp/nakagawa/laplacetrans/convolution3.htm)と書いておられますが、私は正直、そこまで達観した「たたみこみ観」は持っていないので、実践的な「たたみこみ」の説明を書いてみます。
なお、ＣＮＮ（たたみこみニューラルネットワーク）については、本稿では触れません。

## １．なぜ「たたみこみ」が必要なのか
　「たたみこみ」の例としては、音響のエコー（反射）が良く知られています。
下図のように人がいて、１カ所に反射物があったとします。
![図１.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/ae4d8c9f-e64d-a63f-fc73-f8ec896da01a.jpeg)
人と反射物の距離がＬで音速がVcだとすると、人が発生した音声は、反射物で反射し、2L/Vc 秒後に人に戻ってきます（２は往復なので）。
また、図２のように２カ所に反射物があり、人からの距離が２つの反射物までの距離が、それぞれLa，Lbだとしますと、それぞれの反射物からの音は、2La/Vc、2Lb/Vc 秒後に返ってきます。
![図2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/a077d960-0a4a-9540-233a-e82493c5bbab.jpeg)
それぞれから返ってくる音の大きさが同じで、波形に崩れがない（線形）とします。
このとき、２つの反射物からの音は、時間をずらして加算されます。
これをMATLAB的に書くと、

```matlab:sum2point.m
wave=load('wave.dat'); %波形データを読み込む
wavelen=length(wave); %波形データの点数

fs=10e3; % サンプリング周波数[Hz]
Vc=340; % 音速[m/sec]

La=4; % A点までの距離[m]
Lb=6; % B点までの距離[m]

Pa=round(2*La/Vc*fs); % A点までの往復の時間をサンプリング点数に換算
Pb=round(2*Lb/Vc*fs); % B点までの往復の時間をサンプリング点数に換算

WaveAll=zeros(1,600); % 合成波形を収納するためのゼロ配列

figure(1);

subplot(3,1,1);
plot(WaveAll);
text(20,0.5,'ゼロ配列');

WaveAll(Pa:(Pa+wavelen-1))=WaveAll(Pa:(Pa+wavelen-1))+wave;
subplot(3,1,2);
plot(WaveAll);
text(20,0.5,'Ａ点の反射を加える');

WaveAll(Pb:(Pb+wavelen-1))=WaveAll(Pb:(Pb+wavelen-1))+wave;
subplot(3,1,3);
plot(WaveAll);
text(20,0.5,'さらにＢ点の反射を加える');
```
実行すると、下記のようなグラフが表示されます。
![figure1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/9183ecb4-1a98-bfc3-aa08-f348d79958a5.jpeg)
また、A点からの距離と、B点からの距離が近い場合は、波形が重なりますが、重なったところは両者の波形が加算された形になります。
![figure2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/7e3059d8-7b0c-57c8-4b55-24c854b9ec06.jpeg)
ところで、上で説明した演算は、けっこう面倒くさいものがあります。例えて言うなら、
　　　y = ( 2 * a ) + ( 4.5 * a ) + ( 11 * a )
でａの値が与えられたときに、カッコの中から計算しちゃうような面倒くささがあります。
　　　y = ( 2 + 4.5 + 11 ) * a
で計算すべきでしょう。
たたみこみ演算は、これに近い感じです（あくまでも「感じ」ですが）。

##　２．MATLABの「たたみこみ関数」convを使う
　前章の２点からの反射の例では、波形そのものは同じで、返ってくる時間だけが異なっていました（反射信号の振幅も同じ）。
そこで、伝搬時間Pa, Pb（正確には伝搬時間をサンプリング周波数で割ったもの）と振幅の１だけを用いて、これら（遅延と振幅）を表現するインパルス列を作成します。
![図３.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/99d1e4ba-ee8c-7755-6169-e040da5aa29b.jpeg)

このインパルス列と信号波形の畳み込みを行ないます。
　　WaveOut = conv ( WaveIn , Inp ) ; 
というような使い方で、スクリプトは、例えば、

```matlab:ex_conv.m
wave=load('wave.dat'); %波形データを読み込む

fs=10e3; %サンプリング周波数[Hz]
Vc=340; %音速[m/sec]

La=4; %A点までの距離[m]
Lb=6; %B点までの距離[m]

Pa=round(2*La/Vc*fs); %A点までの往復の時間をサンプリング点数に換算
Pb=round(2*Lb/Vc*fs); %A点までの往復の時間をサンプリング点数に換算

WaveAll=zeros(1,600); %合成波形を収納するためのゼロ配列

figure(2);

WaveAll(Pa)=1;
subplot(4,1,1);
plot(WaveAll);
ylim([-0.2 1.2]);
text(20,0.5,'Ａ点反射のインパルス列');

WaveAllA=conv(wave,WaveAll); %たたみこみ
subplot(4,1,2);
plot(WaveAllA);
xlim([0 600]);
text(20,0.5,'Ａ点の反射のみ');

WaveAll(Pb)=1;
subplot(4,1,3);
plot(WaveAll);
ylim([-0.2 1.2]);
text(20,0.5,'Ａ点Ｂ点反射のインパルス列');

WaveAll=conv(wave,WaveAll); %たたみこみ
subplot(4,1,4);
plot(WaveAll);
xlim([0 600]);
text(20,0.5,'Ａ点、Ｂ点の反射');
```
実行すると下のようなグラフが出ます。
![figure3.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/32a2734d-dbf8-0f5e-f070-a59f503149eb.jpeg)
また、Ａ点からの反射がＢ点からの２倍という場合は、インパルスの振幅を変えることで対応できます。
![figure3a.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/0a09979b-526d-0dcf-8f18-9807e010a760.jpeg)
ここで、インパルス列とのたたみこみということで、説明しましたが、離散時間系においては、どのような波形（関数）もインパルスの組み合わせです。
よって、よく出てくる離散値のたたみこみの式、
![式１.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/09f869ea-f2e9-472f-8bed-b9b9c803883d.jpeg)
と同じことになります。
ということで、「たたみこみ」が何か、どうやるかのざっくりとした説明を終わります。


##　３．多数点からの反射をいかに計算するか

前章までで、シンプルな解説は終わったので、ここからは書きたいように書いてみましょう。
反射物として板みたいなものを考えます。板はＸＹ平面上にあるとします。
![図4.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/e417c1bc-de64-858a-ce6a-db7b70e4b71b.jpeg)
板の大きさ、観測点の位置関係は、次のとおりとします。
![図5.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/6ffc2a85-b0b5-5b4b-5818-352372f02445.jpeg)
こういった場合には、セオリーとして、まず反射板をメッシュに切ります（シミュレーションの基本？）。
信号の波長とか演算精度と相談ですが、ここは0.1m格子でメッシュを切ってみましょう。
meshgridという関数が便利です。

```matlab
x=-2.5:0.1:2.5;% X方向メッシュ設定
y=-1.5:0.1:1.5;% Y方向メッシュ設定

[X,Y]=meshgrid(x,y);% メッシュ生成
```
これで、メッシュ各点のＸ座標、Ｙ座標の２次元配列ができました。
メッシュ各点から観測点Ｐまでの距離を求めます。
abs関数を用います。
![図6.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/ecdb9a97-279a-e36a-c9b7-ec24b211cc45.jpeg)

```matlab
P=6; %観測点のＺ座標[m]
Lxyz=abs(abs(X+1i*Y)+1i*P);

%%% 距離をサンプリング点に換算する
fs=10e3; %サンプリング周波数[Hz]
Vc=340; %音速[m/sec]
Pxyz=round(2*Lxyz/Vc*fs); % ２は往復
```
abs関数は絶対値計算と、複素数の大きさ計算に使えます。pythonなどでも使えるようですが。
![figure8.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/ac65999f-5717-8bc4-0262-f49fe8c061db.jpeg)
これで、たたみこみの元になる時間遅延データが完成しました。
問題はこの後、インパルス列をどうやって生成するかです。
forループを回して「集計」するのが一般的かと思いますが、ここはヒストグラム（正しくはhistcounts)を使います。入力の振幅が揃っていないと使えないとか制約はありますが、個人的には「推し」です。

```matlab
Pedges=0:600; %ヒストグラムの境界の位置
[N,edges]=histcounts(Pxyz,Pedges);

subplot(2,1,1);
histogram(Pxyz,Pedges);
xlim([0 600]);
text(20,60,'インパルス列');

WaveAll=conv(wave,N); %たたみこみ
subplot(2,1,2);
plot(WaveAll);
xlim([0 600]);
text(20,200,'たたみこみ波形');
```
下図のようなグラフが表示されます。

![figure4.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/434947/1901638e-d2fc-c2ef-9082-51708777a747.jpeg)


###終わりに

「たまにはマジメなことを書くか」と言っておられる人もいたので、私もマジメに書いてみましたが、いかがでしたでしょうか。
ということで、終わります。
お疲れさまでした。