信号の平滑化

ライブスクリプトを開く

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

目的

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

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

load bostemp
tempH = (1:31*24)'/24;
plot(tempH,tempC)
addTitleAndAxisLabelsBosttempHelperFcn

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

移動平均フィルター

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

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

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

avg24hTempC = filter(coeff24hMA,1,tempC);
plot(tempH,[tempC avg24hTempC])
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly Temperature" "24-Hour Average (Delayed)"])

Figure contains an axes object. The axes object with title Dry Bulb Temperature at Boston Logan International Airport (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains 2 objects of type line. These objects represent Hourly Temperature, 24-Hour Average (Delayed).

フィルター処理された出力は約 12 時間分遅れています。出力に遅延が生じるのは、移動平均フィルターに固有の群遅延があるためです。詳細については、grpdelayを参照してください。

長さ $N$ の対称フィルターは $(N - 1) / 2$ サンプルの遅れをもちます。24 時間平均プロットにおける遅延を手動で補正するには、フィルター遅延を計算し、信号サンプルを "x" 軸に沿って移動させます。

fDelay = (length(coeff24hMA)-1)/2;
plot(tempH,tempC,tempH-fDelay/24,avg24hTempC)
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly Temperature" "24-Hour Average (Delay Compensated)"])

平均の差の抽出

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

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

plot(1:24,mean(deltaTempC))
axis tight
title("Mean Temperature Difference from 24-Hour Average")
xlabel("Hour of Day (Since Midnight)")
ylabel("Temperature Difference (\circC)")

Figure contains an axes object. The axes object with title Mean Temperature Difference from 24-Hour Average, xlabel Hour of Day (Since Midnight), ylabel Temperature Difference ( degree C) contains an object of type line.

ピーク包絡線の抽出

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

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

plot(tempH,[tempC envHigh envMean envLow])
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly" "High" "Mean" "Low"],Location="best")

Figure contains an axes object. The axes object with title Dry Bulb Temperature at Boston Logan International Airport (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains 4 objects of type line. These objects represent Hourly, High, Mean, Low.

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

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

別の一般的なフィルターは ${[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(tempH,tempC,tempH-fDelay/24,binomialMA)
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly Temperature" "Binomial Weighted Average"])

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

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

alpha = 0.45;
exponentialMA = filter(alpha,[1 alpha-1],tempC);
plot(tempH,tempC, ...
     tempH-fDelay/24,binomialMA, ...
     tempH-1/24,exponentialMA)
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly Temperature" "Binomial Weighted Average" ...
       "Exponential Weighted Average"],Location="northwest")

Figure contains an axes object. The axes object with title Dry Bulb Temperature at Boston Logan International Airport (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains 3 objects of type line. These objects represent Hourly Temperature, Binomial Weighted Average, Exponential Weighted Average.

フィルター処理された出力の違いをより詳細に確認するために、"x" 軸の範囲を 1 日にしてビューをズームインします。

xlim([3 4])

Figure contains an axes object. The axes object with title Dry Bulb Temperature at Boston Logan International Airport (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains 3 objects of type line. These objects represent Hourly Temperature, Binomial Weighted Average, Exponential Weighted Average.

Savitzky-Golay フィルター

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

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

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

cubicMA   = sgolayfilt(tempC,3,7);
quarticMA = sgolayfilt(tempC,4,7);
quinticMA = sgolayfilt(tempC,5,9);
plot(tempH,[tempC cubicMA quarticMA quinticMA])
addTitleAndAxisLabelsBosttempHelperFcn
axis([3 5 -5 2])
legend(["Hourly Temperature" "Cubic-Weighted MA" ...
    "Quartic-Weighted MA" "Quintic-Weighted MA"],Location="southeast")

Figure contains an axes object. The axes object with title Dry Bulb Temperature at Boston Logan International Airport (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains 4 objects of type line. These objects represent Hourly Temperature, Cubic-Weighted MA, Quartic-Weighted MA, Quintic-Weighted MA.

リサンプリング

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

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

load openloop60hertz
Fs = 1000;
t = (0:numel(openLoopVoltage)-1)/Fs;
plot(t,openLoopVoltage)
addTitleAndAxisLabelsVoltHelperFcn("Open-Loop Voltage Measurement")

Figure contains an axes object. The axes object with title Open-Loop Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains an object of type line.

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

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

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

Vavg5882 = sgolayfilt(openLoopVoltage,1,17);
plot(t,Vavg5882)
addTitleAndAxisLabelsVoltHelperFcn("Open-Loop Voltage Measurement")
legend("With Moving-Average Filter at 58.82 Hz",Location="southeast")

Figure contains an axes object. The axes object with title Open-Loop Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains an object of type line. This object represents With Moving-Average Filter at 58.82 Hz.

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

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

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

Frs = 1020;
Vrs = resample(openLoopVoltage,Frs,Fs);
trs = (0:numel(Vrs)-1)/Frs;
VrsAvg6000 = sgolayfilt(Vrs,1,17);
plot(t,Vavg5882,trs,VrsAvg6000)
addTitleAndAxisLabelsVoltHelperFcn("Open-Loop Voltage Measurement")
legend("With Moving-Average Filter at " + [58.82 60] +  " Hz", ...
    Location="southeast")

Figure contains an axes object. The axes object with title Open-Loop Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains 2 objects of type line. These objects represent With Moving-Average Filter at 58.82 Hz, With Moving-Average Filter at 60 Hz.

フィルター処理された出力の違いをより詳細に確認するために、"x" 軸の範囲を 0.5 秒にしてビューをズームインします。

xlim([0.75 1.25])

Median Filter

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

load clockex
yMA = conv(x,ones(5,1)/5,"same");
ySG = sgolayfilt(x,3,5);
plot(t,x,t,yMA,t,ySG)
addTitleAndAxisLabelsVoltHelperFcn("Clock-Signal Voltage Measurement")
legend(["Original Signal" ...
    "With Moving-Average Filter" "With Savitzky-Golay Filter"])

Figure contains an axes object. The axes object with title Clock-Signal Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains 3 objects of type line. These objects represent Original Signal, With Moving-Average Filter, With Savitzky-Golay Filter.

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

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

yMedFilt = medfilt1(x,5,"truncate");
plot(t,x,t,yMedFilt)
addTitleAndAxisLabelsVoltHelperFcn("Clock-Signal Voltage Measurement")
legend(["Original Signal" "With Median Filter"])

Figure contains an axes object. The axes object with title Clock-Signal Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains 2 objects of type line. These objects represent Original Signal, With Median Filter.

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

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

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

load train
tTrain = 1/Fs*(0:numel(y)-1);
y(1:400:end) = 2.1;
plot(tTrain,y)
addTitleAndAxisLabelsVoltHelperFcn("Train-Whistle Voltage Measurement")

Figure contains an axes object. The axes object with title Train-Whistle Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains an object of type line.

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

yFilt = medfilt1(y,3);
tFilt = 1/Fs*(0:numel(yFilt)-1);
hold on
plot(tFilt,yFilt)
hold off
legend(["Original" "Filtered"] + "Signal")

Figure contains an axes object. The axes object with title Train-Whistle Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains 2 objects of type line. These objects represent OriginalSignal, FilteredSignal.

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

hampel(y,13)
legend(Location="best")
addTitleAndAxisLabelsVoltHelperFcn("Train-Whistle Voltage Measurement")

Figure contains an axes object. The axes object with title Train-Whistle Voltage Measurement, xlabel Time (s), ylabel Voltage (V) contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent original signal, filtered signal, outliers.

付録: 補助関数

addTitleAndAxisLabelsBosttempHelperFcn 補助関数は、この例で表示される温度プロットを含む図にタイトルと軸のラベルを追加します。

function addTitleAndAxisLabelsBosttempHelperFcn
    title({"Dry Bulb Temperature" ...
        "at Boston Logan International Airport (source: NOAA)"})
    xlabel("Time elapsed from Jan 1, 2011 (days)")
    ylabel("Temperature (\circC)")
    axis tight
end

addTitleAndAxisLabelsVoltHelperFcn 補助関数は、この例で表示される電圧プロットを含む図にタイトルと軸のラベルを追加します。

function addTitleAndAxisLabelsVoltHelperFcn(titleString)
    title(titleString)
    xlabel("Time (s)")
    ylabel("Voltage (V)")
end

参考文献

[1] 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