Ccmmutty logo
Commutty IT
・4 min read

[MATLAB]加算平均処理とデジタル微分

MATLAB
https://cdn.magicode.io/media/notebox/05df1b48-7233-46ad-8b3d-03f53fa36e04.jpeg

はじめに

今回は加算平均処理、デジタル微分処理についてまとめる。

実行環境

  • Ubuntu20.04
  • MATLAB R2022a

加算平均とは

同じ計測を複数回繰り返し、その平均を求める。
→ノイズが減少し、信号成分を抽出できる。

  • 雑音が加法的でランダムかつ平均値ゼロ
  • 信号成分は繰り返し計測しても不変 上記2つの仮定が成立する必要が有る。
一度の測定では雑音成分が大きく信号成分を評価しづらく、複数回測定を繰り返しても信号成分は常に同じである場合に適している。

加算平均処理は以下の式で表すことができる。
x1(t)=S1(t)+n1(t)x2(t)=S2(t)+n2(t)x3(t)=S3(t)+n3(t)x4(t)=S4(t)+n4(t)...xi(t)=Si(t)+ni(t)x_{1}\left( t\right) =S_{1}\left( t\right) +n_{1}\left( t\right) \\ x_{2}\left( t\right) =S_{2}\left( t\right) +n_{2}\left( t\right) \\ x_{3}\left( t\right) =S_{3}\left( t\right) +n_{3}\left( t\right) \\ x_{4}\left( t\right) =S_{4}\left( t\right) +n_{4}\left( t\right) \\ ...\\ x_{i}\left( t\right) =S_{i}\left( t\right) +n_{i}\left( t\right)
これらより、M回の加算平均は
1M∑i=1mxi(t)=1M∑i=1Msi(t)+1M∑i=1Mni(t)\dfrac{1}{M}\sum ^{m}_{i=1}x_{i}\left( t\right) =\dfrac{1}{M}\sum ^{M}_{i=1}s_{i}\left( t\right) +\dfrac{1}{M}\sum ^{M}_{i=1}n_{i}\left( t\right)

加算平均処理の実装

load handel;
s = y(1:20000); % 音声データの最初の20000点を切り出し
N = length(s); % 切り出した音声データの⻑さをNに格納
n = randn(N,1); % N点分の擬似正規白色雑音(平均0,分散1)を生成
x = s + n; % 音声データにノイズを重畳
sound(s, Fs); pause % ノイズなしハレルヤの再生
sound(x, Fs); pause % ノイズに埋もれたハレルヤの再生
x_ave = s*0; % 加算平均結果を格納するx_aveを初期化
M = 10000; % 加算平均回数の指定

for i = 1:M
ni = randn(N,1); % i回目の計測時に重畳するノイズ
xi = s + ni; % i回目の計測データ(信号+雑音)
x_ave = x_ave + xi; % xiをx_aveにM回加算
end
x_ave = x_ave / M; % x_aveをMで割って平均を計算

sound(x_ave, Fs); % 加算平均結果の再生

微分

時刻tkt_kにおける微分係数y(t)y(t)は
 y(tk)=dx(t)dt∣t=tk=lim⁡Δt→0Δx(t)Δt∣t=tk\displaystyle\ y(t_k) = \frac{dx(t)}{dt}|_{t=t_k} = \lim _{\Delta t\rightarrow 0}\dfrac{\Delta x\left( t\right) }{\Delta t}| _{t=t_{k}}
これは接線の傾きになる。
微分とは、全てのtにおける微分係数を求めることである。

時刻tkt_kにおいて微分係数がy(tk)=0y(t_k) = 0の場合、接線の傾きは0でその時刻でx(t_k)は増減していないことがわかる。

要するに時刻tkt_kにピークを持つということ。
微分は波形のピーク(時刻とその値)を求めるために使える。

デジタル微分

kk点目におけるデジタル信号x(k)x(k)の微分係数y(k)y(k)は下式によって求められる。
y(k)=−cx(k−1)+cx(k+1)y(k) = -c x(k-1) +c x(k+1)
=c(x(k+1)−x(k−1))= c (x(k+1) - x(k-1))
ccは標本化周波数によって決まる定数

上式はkkの1時刻後と前の差を取ってcc倍する。
これにより、時刻kkにおけるxxの傾きが求められ、同じ処理を全ての時刻kで行う事によりx(n)x(n)のデジタル微分y(n)y(n)が求められる。

デジタル微分は高周波雑音を増幅することができる。

デジタル微分処理の実装

t= 0:1/1000:1;  % fs=1000で時刻データを1秒間生成
x = sin(2*pi*2*t);  % 2Hzの正弦波を生成し，xに代入
x(75:78) = x(75:78)+0.005;  % xの75~78点目に突発性ノイズ0.005を付加
x(575:578) = x(575:578)+0.002;  % 同575~578点目に0.002を付加

fs = 1000;
c = fs/2;
N = length(x);
y = x*0;
for k=2:N-1
    y(k) = c*(x(k+1) - x(k-1));
end
y(1) = (-x(1) + x(2))*fs;
y(N) = (-x(N-1) + x(N))*fs;


plot(t,x,'r-', t,y,'b-')

最後に

加算平均処理は信号を変化させることなく信号成分の抽出ができる。 デジタル微分処理はノイズの高周波成分を増幅することができるので、信号内に潜在している突発性ノイズの検出等に役立つ。

Discussion

コメントにはログインが必要です。