신호 평활화

라이브 스크립트 열기

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

목적

평활화는 매끄럽게 함으로써 데이터에서 중요하지 않은 것(잡음)을 제거하고 그 패턴을 알아내는 방법입니다. 필터링을 사용하여 이 평활화를 수행합니다. 평활화의 목적은 값이 천천히 변화하도록 만들어, 데이터의 추세(경향)를 더 쉽게 파악할 수 있게 하는 것입니다.

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

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

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

이동평균 필터

가장 단순한 형태에서, 길이 $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)"])

평균 차이 추출하기

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

[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}$ 의 이항 전개를 따릅니다. 이 유형의 필터는 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(tempH,tempC,tempH-fDelay/24,binomialMA)
addTitleAndAxisLabelsBosttempHelperFcn
legend(["Hourly Temperature" "Binomial Weighted Average"])

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

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

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

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

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

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

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.

리샘플링

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

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

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.

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

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

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

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.

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

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

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

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

중앙값 필터

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

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.

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

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

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 Filter)를 통해 이상값 제거하기

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

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

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.

추가한 각 스파이크의 기간이 샘플 한 개에 해당하기 때문에 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(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 헬퍼 함수는 이 예제에 표시된 온도 플롯이 포함된 Figure에 제목과 축 레이블을 추가합니다.

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 헬퍼 함수는 이 예제에 표시된 전압 플롯이 포함된 Figure에 제목과 축 레이블을 추가합니다.

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