Main Content

신호 평활화

이 예제에서는 이동평균 필터와 리샘플링을 사용하여 시간별 온도 측정값에 대한 태양이 떠 있는 주간 시간의 주기적 성분의 영향을 분리하고 개루프 전압 측정값에서 바람직하지 않은 회선 잡음을 제거하는 방법을 보여줍니다. 또한, 중앙값 필터를 사용하여 경계값을 보존하면서 클록 신호의 레벨을 평활화하는 방법도 보여줍니다. 추가적으로, 햄펄 필터(Hampel Filter)를 사용하여 큰 이상값을 제거하는 방법도 보여줍니다.

목적

평활화는 매끄럽게 함으로써 데이터에서 중요하지 않은 것(예: 잡음)을 제거하고 그 패턴을 알아내는 방법입니다. 필터링을 사용하여 이 평활화를 수행합니다. 평활화의 목적은 값의 변화를 매끄럽게 하여 데이터의 추세(경향)를 알기 쉽게 하는 것입니다.

입력 데이터를 검증할 때 신호의 추세를 파악하고자 데이터의 평활화를 수행할 수도 있습니다. 이 예제에서는 2011년 1월 한 달 동안 보스톤의 로간 공항(Logan Airport)에서 한 시간마다 측정한 온도 측정값(단위: 섭씨)을 사용합니다.

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

Figure contains an axes object. The axes object with title Logan Airport Dry Bulb Temperature (source: NOAA), xlabel Time elapsed from Jan 1, 2011 (days), ylabel Temperature ( degree C) contains an object of type line.

태양이 떠 있는 주간 시간이 온도 측정값에 미친 영향을 시각적으로 알 수 있습니다. 한 달간의 일일 온도 변화에만 관심이 있는 경우 시간별 변동은 잡음을 일으키고 일일 변화를 파악하기 어렵게 만들 수 있습니다. 태양이 떠 있는 주간 시간의 영향을 제거하기 위해 이제 이동평균 필터를 사용하여 데이터를 평활화합니다.

이동평균 필터

가장 단순한 형태로, 길이 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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")

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

필터 지연

필터링된 출력값이 약 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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")

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

평균 차이 추출하기

또는 이동평균 필터를 사용하여 태양이 떠 있는 주간 시간이 전반적인 온도에 어떻게 영향을 미치는지에 대한 더 나은 추정값을 얻을 수도 있습니다. 이를 위해 먼저 시간별 온도 측정값에서 평활화된 데이터를 뺍니다. 그런 다음, 차이 데이터를 일수로 분할하고 한 달 전체 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)")

Figure contains an axes object. The axes object with title Mean temperature differential from 24 hour average, xlabel Hour of day (since midnight), ylabel Temperature difference ( degree C) contains an object of type line.

피크 포락선 추출하기

경우에 따라 온도 신호가 매일 얼마나 높이 또는 낮게 변하는지에 대한 매끄러운 변동 추정을 얻고자 할 수 있습니다. 이를 위해 envelope 함수를 사용하여 24시간 한 주기에서 검출한 극댓값과 극솟값을 연결할 수 있습니다. 이 예제에서는 각각의 극댓값과 극솟값 사이의 기간이 최소 16시간이 되도록 합니다. 또한, 두 극값 사이의 평균을 구하여 높은 값과 낮은 값의 추세(경향)를 파악할 수도 있습니다.

[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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")

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

가중 이동평균 필터

다른 종류의 이동평균 필터는 각 샘플에 균일하게 가중치를 적용하지 않습니다.

또 다른 일반적인 필터는 [1/2,1/2]n의 이항 전개를 따릅니다. 이 유형의 필터는 n이 큰 값인 경우 정규 곡선을 근사합니다. 이 필터는 n이 작은 경우 고주파 잡음을 필터링하는 데 유용합니다. 이항 필터의 계수를 구하려면 [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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")

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

가우스 확장 필터와 다소 유사한 또 다른 필터는 지수 이동평균 필터입니다. 이 유형의 가중 이동평균 필터는 생성하기가 쉬우며 큰 윈도우 크기를 필요로 하지 않습니다.

지수적으로 가중치가 적용된 이동평균 필터를 0과 1 사이의 알파 파라미터를 사용하여 조정합니다. 알파 값이 클수록 덜 매끄러워집니다.

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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")

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

특정 일의 측정값을 확대합니다.

axis([3 4 -5 2])

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

사비츠키-골레이 필터(Savitzky-Golay Filter)

데이터를 평활화하면 극값이 다소 잘리는 것을 알 수 있습니다.

신호를 좀 더 면밀하게 추적하려면 가중 이동평균 필터를 사용하여 지정된 차수의 다항식을 지정된 수의 샘플에 최소제곱 피팅으로 근사해 볼 수 있습니다.

편의상, 함수 sgolayfilt를 사용하여 사비츠키-골레이 평활화 필터를 구현할 수 있습니다. sgolayfilt를 사용하려면 홀수 길이의 데이터 세그먼트와 함께, 세그먼트 길이보다 작은 다항식 차수를 지정해야 합니다. sgolayfilt 함수는 평활화 다항식 계수를 내부적으로 계산하여 지연을 조정하고, 데이터 레코드의 시작 부분과 끝부분에 나타나는 과도(Transient) 현상을 처리합니다.

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("Temperature (\circC)")
xlabel("Time elapsed from Jan 1, 2011 (days)")
title("Logan Airport Dry Bulb Temperature (source: NOAA)")
axis([3 5 -5 2])

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

리샘플링

경우에 따라 이동평균을 올바르게 적용하기 위해 신호를 리샘플링하는 것이 유용할 수 있습니다.

다음 예제에서는 60Hz AC 전력 공급선으로부터 잡음의 영향을 받고 있는 아날로그 기기의 입력에서 개루프 전압을 샘플링했습니다. 1kHz 샘플링 레이트로 전압을 샘플링했습니다.

load openloop60hertz
fs = 1000;
t = (0:numel(openLoopVoltage)-1) / fs;
plot(t,openLoopVoltage)
ylabel("Voltage (V)")
xlabel("Time (s)")
title("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.

이동평균 필터를 사용하여 회선 잡음 영향을 제거해 보겠습니다.

균일 가중 이동평균 필터를 생성하면 필터 지속 시간과 관련하여 주기적인 성분이 모두 제거됩니다.

1000Hz로 샘플링했을 때 60Hz의 1주기에는 대략적으로 1000 / 60 = 16.667개의 샘플이 있습니다. 이를 "올림"하여 17점 필터를 사용해 보겠습니다. 그러면 기본주파수 1000Hz / 17 = 58.82Hz에서 필터링을 극대화할 수 있습니다.

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")

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 Moving average filter operating at 58.82 Hz.

전압이 매우 매끄러워졌지만 여전히 60Hz의 작은 리플이 남아 있는 것을 알 수 있습니다.

이동평균 필터로 60Hz 신호의 1주기 전체를 캡처하도록 신호를 리샘플링하면 리플을 상당히 줄일 수 있습니다.

17 * 60Hz = 1020Hz로 신호를 리샘플링하면 17점 이동평균 필터를 사용하여 60Hz의 회선 잡음을 제거할 수 있습니다.

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")

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 Moving average filter operating at 60 Hz.

중앙값 필터

이동평균 필터, 가중 이동평균 필터, 사비츠키-골레이 필터(Savitzky-Golay Filter)는 필터 처리한 모든 데이터를 평활화 처리하여 매끄럽게 만듭니다. 그러나 데이터 평활화가 항상 필요한 것은 아닙니다. 예를 들어, 클록 신호에서 데이터를 가져오는 경우 이 데이터에 평활화를 원하지 않는 예리한 경계값이 있다면 어떻게 해야 할까요? 이런 경우, 지금까지 설명한 필터는 그다지 도움이 되지 않습니다.

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")

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent original signal, moving average, Savitzky-Golay.

이동평균 필터와 사비츠키-골레이 필터는 각각 클록 신호의 경계값 인근에서 부족보정과 과잉보정을 수행합니다.

경계값을 보존하면서 레벨을 평활화하는 간단한 방법은 다음과 같이 중앙값 필터를 사용하는 것입니다.

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

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent original signal, median filter.

햄펄 필터(Hampel Filter)를 통해 이상값 제거하기

필터의 대부분은 이상값의 영향을 받습니다. 중앙값 필터와 밀접한 관련이 있는 필터는 햄펄 필터입니다. 이 필터는 데이터를 과도하게 평활화하지 않으면서 신호에서 이상값을 제거하는 데 도움이 됩니다.

이를 확인하기 위해 기차 기적의 오디오 녹음을 불러오고 몇 가지 인공적인 잡음 스파이크를 추가합니다.

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

Figure contains an axes object. The axes object contains an object of type line.

추가한 각 스파이크의 기간이 샘플 한 개에 해당하기 때문에 3개 요소에 대한 중앙값 필터를 사용하여 스파이크를 제거할 수 있습니다.

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

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent original signal, filtered signal.

이 필터 사용을 통해 스파이크가 제거되었지만 원래 신호의 데이터 점도 많이 제거되었습니다. 햄펄 필터는 중앙값 필터와 유사하게 동작하지만, 국소 중앙값에서 표준편차의 몇 배 떨어진 값만 대체합니다.

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

Figure contains an axes object. The axes object 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.

원래 신호에서 이상값만 제거되었습니다.

추가 참고 자료

필터링과 리샘플링에 대한 자세한 내용은 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.

참고 항목

| | | |