信号の平滑化

ライブスクリプトを開く

この例では、移動平均フィルターとリサンプリングを使用して、1 時間ごとの温度の測定値について時刻の周期的成分の影響を分離したり、開ループ電圧測定値から望ましくない回線ノイズを取り除いたりする方法を示します。この例では、メディアンフィルターを使用してエッジを保存しながら、クロック信号のレベルを平滑化する方法も説明します。さらに、例では、Hampel フィルターを使用して大きな外れ値を取り除く方法についても説明します。

目的

平滑化とは、データにおける重要パターンを、ノイズなど、重要性の低いものを除去しながら見つけ出す方法です。この平滑化のためにフィルター処理を使用します。平滑化の目的は、値の変化をなだらかにしてデータのトレンドをわかりやすくすることです。

入力データを検証するとき、信号のトレンドを調べるためにデータの平滑化を行うことがあります。ここに、2011 年 1 月の 1 か月間、ボストンのローガン空港で 1 時間間隔で測定された温度が摂氏単位で表示されています。

load bostemp
days = (1:31*24)/24;
plot(days, tempC)
axis tight
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

温度の読み取りにおいて 1 日の時間が与える影響を目視できることに注目してください。その月にわたり日々の温度変化だけに関心がある場合、1 時間ごとの変動はノイズになるだけで毎日の変化を見定めにくくします。時刻の影響を取り除くために、移動平均フィルターを使用してデータを平滑化します。

移動平均フィルター

最も簡単な形では、長さ N の移動平均フィルターは連続した N 個ごとの波形サンプルの平均をとります。

各データ点に移動平均フィルターをかけるには、各点が均等に重み付けされ、合計平均の 1/24 の割合を示すようにフィルター係数を作成します。これにより、24 時間ごとの平均温度が得られます。

hoursPerDay = 24;
coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;

avg24hTempC = filter(coeff24hMA, 1, tempC);
plot(days,[tempC avg24hTempC])
legend('Hourly Temp','24 Hour Average (delayed)','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

フィルター遅延

フィルター出力が約 12 時間分遅れていることに注意します。これは、移動平均フィルターが遅れをもつことによるものです。

長さ N の対称フィルターはどれも (N-1)/2 サンプルの遅れをもちます。この遅れは手作業で調整できます。

fDelay = (length(coeff24hMA)-1)/2;
plot(days,tempC, ...
     days-fDelay/24,avg24hTempC)
axis tight
legend('Hourly Temp','24 Hour Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

平均の差の抽出

また、移動平均フィルターを使用して、時刻が全体の温度にどのように影響しているかをより正確に推定できます。これを行うにはまず、1 時間ごとの温度の測定値から平滑化したデータを減算します。次に、差分データを日ごとに分け、その月の 31 日全体について平均をとります。

figure
deltaTempC = tempC - avg24hTempC;
deltaTempC = reshape(deltaTempC, 24, 31).';

plot(1:24, mean(deltaTempC))
axis tight
title('Mean temperature differential from 24 hour average')
xlabel('Hour of day (since midnight)')
ylabel('Temperature difference (\circC)')

ピーク包絡線の抽出

温度信号の高値と低値が、毎日どのように変化するかについて、なだらかな変動推定を得ることが必要な場合もあります。これを行うには、関数 envelope を使用して、24 時間周期のサブセットにおいて検出された極端に高い値と極端に低い値をつなぎます。この例では、極端に高い値と極端に低い値のそれぞれの間に少なくとも 16 時間の間隔が空くようにしています。極端な 2 つの値の間の平均を取得することによって、高値と低値がどのようなトレンドで変化しているかを把握することもできます。

[envHigh, envLow] = envelope(tempC,16,'peak');
envMean = (envHigh+envLow)/2;

plot(days,tempC, ...
     days,envHigh, ...
     days,envMean, ...
     days,envLow)
   
axis tight
legend('Hourly Temp','High','Mean','Low','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

重み付けされた移動平均フィルター

他の種類の移動平均フィルターでは各サンプルの均等な重み付けは行われません。

別の一般的なフィルターは ${[1 / 2, 1 / 2]}^{n}$ の 2 項展開式に従います。n が大きな値の場合、このタイプのフィルターは正規曲線を近似します。n が小さい場合、これは高周波ノイズの除去に役立ちます。2 項フィルターの係数を見つけるには、 $[1 / 2, 1 / 2]$ とそれ自体の畳み込みを行い、出力と $[1 / 2, 1 / 2]$ の畳み込みを指定された回数繰り返します。この例では、総反復回数は 5 です。

h = [1/2 1/2];
binomialCoeff = conv(h,h);
for n = 1:4
    binomialCoeff = conv(binomialCoeff,h);
end

figure
fDelay = (length(binomialCoeff)-1)/2;
binomialMA = filter(binomialCoeff, 1, tempC);
plot(days,tempC, ...
     days-fDelay/24,binomialMA)
axis tight
legend('Hourly Temp','Binomial Weighted Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

ガウスフィルターとやや似た別のフィルターとして指数移動平均フィルターがあります。このタイプの重み付けされた移動平均フィルターは作成が容易で、大きなウィンドウサイズを必要としません。

指数関数的に重み付けされた移動平均フィルターを 0 から 1 までの alpha パラメーターで調整します。alpha の値が大きいほど滑らかさが少なくなります。

alpha = 0.45;
exponentialMA = filter(alpha, [1 alpha-1], tempC);
plot(days,tempC, ...
     days-fDelay/24,binomialMA, ...
     days-1/24,exponentialMA)

axis tight
legend('Hourly Temp', ...
       'Binomial Weighted Average', ...
       'Exponential Weighted Average','location','best')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')

1 日の測定値を拡大します。

axis([3 4 -5 2])

Savitzky-Golay フィルター

データの平滑化によって、極端な値がいくらか切り捨てられています。

信号をもう少し厳密に追跡するには、指定の次数の多項式を、最小二乗により指定したサンプル数に当てはめようとする重み付けされた移動平均フィルターを使用できます。

便宜上、関数 sgolayfilt を使用して Savitzky-Golay 平滑化フィルターを実装できます。sgolayfilt を使用するには、奇数長のデータセグメントと、セグメント長より厳密に小さい多項式の次数を指定します。関数 sgolayfilt は平滑化多項式係数を内部で計算し、遅延調整して、データの最初と最後で生じる過渡的な影響の処理を行います。

cubicMA   = sgolayfilt(tempC, 3, 7);
quarticMA = sgolayfilt(tempC, 4, 7);
quinticMA = sgolayfilt(tempC, 5, 9);
plot(days,[tempC cubicMA quarticMA quinticMA])
legend('Hourly Temp','Cubic-Weighted MA', 'Quartic-Weighted MA', ...
       'Quintic-Weighted MA','location','southeast')
ylabel('Temp (\circC)')
xlabel('Time elapsed from Jan 1, 2011 (days)')
title('Logan Airport Dry Bulb Temperature (source: NOAA)')
axis([3 5 -5 2])

リサンプリング

移動平均を適切に行うために、信号のリサンプリングをした方がよい場合もあります。

次の例は、ノイズの影響を受けているアナログ計器の、60 Hz 電力供給ラインからの入力における開ループ電圧のサンプリング結果です。1 kHz のサンプリングレートで電圧をサンプリングしました。

load openloop60hertz
fs = 1000;
t = (0:numel(openLoopVoltage)-1) / fs;
plot(t,openLoopVoltage)
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')

移動平均フィルターを使用して回線ノイズの影響を取り除いてみます。

一律に重み付けされた移動平均フィルターを作成した場合、フィルターの継続時間における周期性成分がいずれも取り除かれます。

1000 Hz でサンプリングした場合、60 Hz の 1 サイクル中のサンプル数は約 1000 / 60 = 16.667 です。この数を丸めて、17 点フィルターを使用してみます。これにより、基本周波数 1000 Hz / 17 = 58.82 Hz においてフィルター処理が最大になります。

plot(t,sgolayfilt(openLoopVoltage,1,17))
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')
legend('Moving average filter operating at 58.82 Hz', ...
       'Location','southeast')

電圧はかなり平滑されていますが、まだ微小な 60 Hz リップルが残っていることがわかります。

60 Hz の信号 1 サイクル分を完全に移動平均フィルターで捕捉できるように信号をリサンプリングすれば、リップルを大幅に減らすことができます。

17 × 60 Hz = 1020 Hz で信号をリサンプリングすれば、17 点移動平均フィルターを使用して 60 Hz 回線ノイズを取り除くことができます。

fsResamp = 1020;
vResamp = resample(openLoopVoltage, fsResamp, fs);
tResamp = (0:numel(vResamp)-1) / fsResamp;
vAvgResamp = sgolayfilt(vResamp,1,17);
plot(tResamp,vAvgResamp)
ylabel('Voltage (V)')
xlabel('Time (s)')
title('Open-loop Voltage Measurement')
legend('Moving average filter operating at 60 Hz', ...
    'Location','southeast')

Median Filter

移動平均フィルター、加重移動平均フィルターおよび Savitzky-Golay フィルターは、フィルター処理するすべてのデータを平滑化します。ただし、これが常に必要であるとは限りません。たとえば、クロック信号から取得されたデータについて、鋭いエッジの平滑化が望ましくない場合があるとします。これまでに説明したフィルターは、このような処理には不向きです。

load clockex

yMovingAverage = conv(x,ones(5,1)/5,'same');
ySavitzkyGolay = sgolayfilt(x,3,5);
plot(t,x, ...
     t,yMovingAverage, ...
     t,ySavitzkyGolay)

legend('original signal','moving average','Savitzky-Golay')

移動平均フィルターと Savitzky-Golay フィルターでは、クロック信号のエッジ付近で、それぞれ補正不足と過剰補正が生じます。

エッジを保存しながらレベルを平滑化する簡単な方法は、メディアンフィルターを使用することです。

yMedFilt = medfilt1(x,5,'truncate');
plot(t,x, ...
     t,yMedFilt)
legend('original signal','median filter')

Hampel フィルターによる外れ値の除去

フィルターの多くは外れ値の影響を受けます。メディアンフィルターと近い関係にあるのが Hampel フィルターです。このフィルターは、データを過剰に平滑化することなく信号から外れ値を除去するのに役立ちます。

これを確認するため、列車の警笛のオーディオ録音を読み込み、人工的なノイズスパイクをいくらか加えます。

load train
y(1:400:end) = 2.1;
plot(y)

導入した各スパイクの持続時間は 1 サンプル分のみであるため、このスパイクの除去には要素が 3 つのみのメディアンフィルターを使用します。

hold on
plot(medfilt1(y,3))
hold off
legend('original signal','filtered signal')

フィルターによりスパイクは除去されましたが、同時に元の信号から多数のデータ点が除去されています。Hampel フィルターは、メディアンフィルターと同様に動作しますが、局所的な中央値から標準偏差の数倍離れた値のみを置き換えます。

hampel(y,13)
legend('location','best')

外れ値のみが元の信号から除去されています。

参考情報

フィルター処理およびリサンプリングの詳細については、Signal Processing Toolbox を参照してください。

参考文献: Kendall, Maurice G., Alan Stuart, and J. Keith Ord. The Advanced Theory of Statistics, Vol. 3: Design and Analysis, and Time-Series. 4th Ed. London: Macmillan, 1983.

参考

envelope | hampel | medfilt1 | resample | sgolayfilt