移動平均

参考文献及びWEB

参考文献

(1)「例題で学ぶ過渡現象」大重力、森本義広、神田一伸　共著　森北出版
(2)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著日経BP社
(3)「ディジタル信号処理の基礎」三上直樹著 CQ出版
(4)「C言語によるディジタル信号処理入門」三上直樹著 CQ出版
(5)「アナログ&ディジタルフィルタ入門」小野浩司著日刊工業
(6)「フーリエの冒険」ヒッポファミリークラブ
(7)「C言語ではじめる医用情報処理」小高和宏著　Ohmsha

あわせて読む

1次フィルタ
フーリエ級数展開
畳込みについて
時間領域と周波数領域
入力信号
平均@分散@標準偏差

移動平均

3点移動平均を作る

番目	元データ	移動平均の値
0	10	-
1	14	-
2	11	(10+14+11)/3=11.67
3	15	(14+11+15)/3=13.33
4	13	(11+15+13)/3=13
5	17	(15+13+17)/3=15

■定義
$y(i) = \frac{1}{N} \sum_{j=i-(N-1)}^{i} x(j)$

データが0番目から始まる場合の3点移動平均

$y(i) = \frac{1}{3} \sum_{j=i-2}^{i} x(j)$

i=0のとき
$y(0) = \frac{1}{3} \sum_{j=-2}^{0} x(j)=\frac{1}{3}(x(-2)+x(-1)+x(0))$

i=1のとき
$y(1) = \frac{1}{3} \sum_{j=-1}^{1} x(j)=\frac{1}{3}(x(-1)+x(0)+x(1))$

i=2のとき
$y(2) = \frac{1}{3} \sum_{j=0}^{2} x(j)=\frac{1}{3}(x(0)+x(1)+x(2))$

i=3のとき
$y(3) = \frac{1}{3} \sum_{j=1}^{3} x(j)=\frac{1}{3}(x(1)+x(2)+x(3))$

■気付いた点
1. 現在の点と過去の2点を使用して移動平均としている。
2. 3点移動平均の場合、各データに1/3の重みがかかっているとも判断できる。
3. データが3つたまる前は、過度現象と判断できる。

「平均」とはデータのでこぼこをならす操作であるので、移動平均すると図のように
データは滑らかになる。
つまり、移動平均は「ローパスフィルタ」であるといえる。

式から考えると、移動平均はインパルス応答(フィルタ係数)
$h[0]=h[1]=h[2]=\frac{1}{3}$
で表される式であると考えることができる。

つまり、フィルタ係数は、

$h[0]=\frac{1}{3}, h[1]=\frac{1}{3}, h[2]=\frac{1}{3}$

$h[n] = \frac{1}{3} \delta[n]+\frac{1}{3}\delta[n-1]+\frac{1}{3}\delta[n-2]$

この式のZ変換は、

定義 $X(z) = \sum_{n=0}^{\infty} x[n] z^{-n}$
から
$H[z] = \sum_{n=0}^{2} h[n]z^{-n} =h[0]+h[1]z^{-1}+h[2]z^{-2}$
$= \frac{1}{3} +\frac{1}{3}z^{-1}+\frac{1}{3}z^{-2}$

周波数特性を求める

$H[z] = \frac{1}{3} +\frac{1}{3}z^{-1}+\frac{1}{3}z^{-2}$

の周波数特性は、

$H[e^{j\omega\Delta t_{s}}] = \frac{1}{3} +\frac{1}{3}e^{-j\omega\Delta t_{s}}+\frac{1}{3}e^{-2j\omega\Delta t_{s}}$

ここで簡単化のためサンプリング周期を1[s]とすれば

$H[e^{j\omega}] = \frac{1}{3} +\frac{1}{3}e^{-j\omega}+\frac{1}{3}e^{-2j\omega}$

$H[e^{j\omega}] = (\frac{1}{3}e^{j\omega} +\frac{1}{3}+\frac{1}{3}e^{-j\omega})e^{-j\omega}$

$= \frac{1+e^{j\omega}+e^{-j\omega}}{3}e^{-j\omega}$

オイラーの公式

$2\cos(\omega)= e^{j\omega}+e^{-j\omega}$
を使って

$= \frac{1+2\cos(\omega)}{3}e^{-j\omega}$ (1)

そしてフィルタの振幅特性と位相特性は、

振幅特性: $| H(e^{j\omega\Delta t_s}) |= \sqrt{ H_{r}^{2}(e^{j\omega\Delta t_s}) + j H_{i}^{2}(e^{j\omega\Delta t_s})}$
位相特性: $\angle (e^{j\omega\Delta t_s}) = \arctan\frac{ H_{i}(e^{j\omega\Delta t_s}) }{ H_{r}(e^{j\omega\Delta t_s})}$

より

振幅特性は
$| H[e^{j\omega}] | = | \frac{1+2\cos(\omega)}{3} | |e^{-j\omega}|$
$| H[e^{j\omega}] | = | \frac{1+2\cos(\omega)}{3} | \times 1$
$| H[e^{j\omega}] | = | \frac{1+2\cos(\omega)}{3} |$

位相特性は
$\angle (e^{j\omega}) = -\omega$

Matlabコード moving_average.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% moving_average sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生データを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);

x=[10,14,11,15,13,17];
[gyo,max]=size(x);
iti=[0:1:max-1];

% 生データを表示する
plot(iti,x,'ro-','linewidth',2);
hold on;

%3点移動平均
data=0;
y=[];
N=3;

for i=[N-1:1:max-1] %0番目から始めるといちいちめんどくさい
  data=0;
  for j=(i-(N-1)):1:i
   data=data+1/N*x(j+1); %x(j+1) x(0)ができないためこういった処理も生じる
  end
  y=[y,data];
end

y
iti2=[N-1:1:max-1]

%3点移動平均をプロットする
plot(iti2,y,'bs-','linewidth',2);

axis([0 5 8 18]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

grid on;

% キャプション
h=legend('元信号',...,
'3点移動平均した信号',...,
2);
set(h,'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

print -djpeg moving_average_sample_t.jpg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 移動平均の周波数特性、位相特性を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
subplot(2,1,1)

%サンプリング周期 1[s]
dt=1;

%サンプリング周波数1Hz
df=1/dt;

%周波数軸を作る
fre = [0:df*2/1000:df*2];

%FFTの絶対値を求める
H=abs((1+2*cos(2*pi*fre*dt))/3);

plot(fre,H,'r-','linewidth',2);
axis([0 df*2  0 Inf]);
grid on;

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2)

%サンプリング周期 1e-3[s]
dt=1e-3;

%サンプリング周波数1kHz
df=1/dt;

%周波数軸を作る
fre = [0:df*2/1000:df*2];

%FFTの絶対値を求める
H=abs((1+2*cos(2*pi*fre*dt))/3);

plot(fre,H,'b-','linewidth',2);
axis([0 df*2  0 Inf]);
grid on;

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Amplitude','Fontsize',15);

print -djpeg moving_average_sinpuku_tokusei.jpg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 移動平均の位相特性を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)

%サンプリング周期 1[s]
dt=1;

%サンプリング周波数1Hz
df=1/dt;

%周波数軸を作る
fre = [0:df*2/1000:df*2];

%FFTの絶対値を求める
H_iso=-(2*pi*fre);

plot(fre,H_iso,'r-','linewidth',2);
axis([0 df*2  -Inf Inf]);
grid on;

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Phase[rad]','Fontsize',15);

print -djpeg moving_average_iso_tokusei.jpg

FFTで移動平均をもとめる

Z変換で移動平均が求められるならば、ある一定量のデータを移動平均する
場合は周波数領域で計算したほうが速い。

Matlabソース

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fft moving_average sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生データを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); 

%サンプリング周期 1[s]
dt=1;

%サンプリング周波数1Hz
df=1/dt;

x=[10,14,11,15,13,17];
[gyo,max]=size(x);
iti=[0:1:max-1]; 

% 生データを表示する
plot(iti,x,'ro-','linewidth',2);
hold on;

%3点移動平均
data=0;
y=[];
N=3;

for i=[N-1:1:max-1] %0番目から始めるといちいちめんどくさい
  data=0;
  for j=(i-(N-1)):1:i
   data=data+1/N*x(j+1); %x(j+1) x(0)ができないためこういった処理も生じる
  end
  y=[y,data];
end

y
iti2=[N-1:1:max-1]

%3点移動平均をプロットする
plot(iti2,y,'bs-','linewidth',2);

axis([0 5 8 18]);
grid on;

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3点FFTで移動平均する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%元信号をフーリエ変換する
fftx=fft(x)

%インパルス応答(フィルタ係数)をフーリエ変換する
hn=[1/3,1/3,1/3]
ffthn=fft(hn,6) %このように書くと6点で足りない点つまりこの場合は3点を0補間できる

yfft=fftx.*ffthn

%逆フーリエ変換
yifft=ifft(yfft)
y
plot(iti,yifft,'m-','linewidth',2);

% キャプション
h=legend('元信号',...,
'3点移動平均した信号',...,
'3点FFT移動平均した信号',...,
2);
set(h,'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

print -djpeg moving_average_fft.jpg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

信号が連続で送信されるときについて考えてみる

2048ポイントだけを計算してみる

filterTest@20091104a.zip

普通に移動平均をかける場合と、FFT移動平均の場合とで比較してみる

上を拡大したのが下図である。

図により最初の10ポイントが両者一致しない(フィルタの過度状態のため)

1200ポイントごと計算してみる(移動平均と移動平均FFT)

filterTest@20091105a.zip

結果は以下のとおり、最初の9つのデータに違いが現れる(過度現象のため)。
そのため、9つの過去データを入れたうえでFFTしなければならない。

 1.2000000e+003 -1.1000000e-002 -6.4000000e-003
 1.2010000e+003 -8.7000000e-003 -5.8000000e-003
 1.2020000e+003 -1.2900000e-002 -9.5000000e-003
 1.2030000e+003 -1.5600000e-002 -1.2700000e-002
 1.2040000e+003 -1.7500000e-002 -1.3200000e-002
 1.2050000e+003 -1.8100000e-002 -1.7000000e-002
 1.2060000e+003 -2.1500000e-002 -1.8100000e-002
 1.2070000e+003 -2.0100000e-002 -1.7500000e-002
 1.2080000e+003 -1.6900000e-002 -1.9600000e-002
 2.4000000e+003 -6.5000000e-003  1.5000000e-003
 2.4010000e+003 -7.5000000e-003  1.0000000e-003
 2.4020000e+003 -5.3000000e-003  3.0000000e-003
 2.4030000e+003 -6.6000000e-003  2.5000000e-003
 2.4040000e+003 -7.2000000e-003 -1.0000000e-003
 2.4050000e+003 -6.2000000e-003  5.0000000e-004
 2.4060000e+003 -5.4000000e-003 -1.6000000e-003
 2.4070000e+003 -5.2000000e-003 -1.9000000e-003
 2.4080000e+003 -3.2000000e-003 -1.1000000e-003
 3.6000000e+003 -6.7000000e-003  1.8000000e-003
 3.6010000e+003 -3.2000000e-003  3.3000000e-003
 3.6020000e+003  1.8000000e-003  5.4000000e-003
 3.6030000e+003 -1.0000000e-004  4.0000000e-003
 3.6040000e+003  2.0000000e-004  2.5000000e-003
 3.6050000e+003  1.5000000e-003  1.4000000e-003
 3.6060000e+003  2.2000000e-003  2.6000000e-003
 3.6070000e+003 -1.7000000e-003 -5.0000000e-004
 3.6080000e+003  1.5000000e-003  0.0000000e+000
 4.8000000e+003 -1.8900000e-002 -3.2000000e-003
 4.8010000e+003 -1.9100000e-002 -4.9000000e-003
 4.8020000e+003 -1.9400000e-002 -7.8000000e-003

■コード

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 生データを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(2,1,1);
%//250Hz
dt=0.004;
load 'input1.txt'
load 'output_moveave.txt'
load 'output_fft_moveave.txt'

%row 行 column 列
[row,column]=size(output_fft_moveave)

time=[0:row-1];

input1=input1(1:row);
output_moveave=output_moveave(1:row);

% 生データを表示する
plot(time,input1,'r-','linewidth',2);
hold on;

axis([-Inf Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on;

% キャプション
h=legend('元のデータ',...,
'3',...,
2);
set(h,'FontSize',15,'FontName','MSゴシック');

xlabel('Time[s]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2); 

% 移動平均データを表示する
plot(time,output_moveave,'ro-','linewidth',2);
hold on;
% 移動平均データを表示する
plot(time,output_fft_moveave,'b+-','linewidth',2);
hold on;

axis([-Inf Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15,'FontName','Century');

grid on;

% キャプション
h=legend('移動平均',...,
'移動平均(FFT)',...,
'3',...,
4);
set(h,'FontSize',15,'FontName','MSゴシック');

xlabel('Time[s]','Fontsize',15,'FontName','Century');
ylabel('Amplitude','Fontsize',15,'FontName','Century')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
result2=[];
for i=1:1:row
   
  if(output_moveave(i) ~= output_fft_moveave(i))
      result=cat(2,i-1,output_moveave(i),output_fft_moveave(i));
      result2=cat(1,result2,result);
  end
end
  
save kentoudata.txt result2 -ascii

print -djpeg move_ave_discuss2.jpg