Main Content

Signal Processing Toolbox로 데이터 필터링하기

저역통과 FIR 필터 – 윈도우 방법

이 예제에서는 명령줄 함수 fir1designfilt와 대화형 방식의 필터 디자이너 앱을 사용하여 FIR 필터를 설계하고 구현하는 방법을 보여줍니다.

이 예제에 사용할 신호를 생성합니다. 이 신호는 가산성 N(0,1/4) 백색 가우스 잡음(AWGN)이 있는 100Hz의 사인파입니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 디폴트 상태로 설정합니다.

rng default

Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));

필터 설계는 차수가 20이고 차단 주파수가 150Hz인 FIR 저역통과 필터를 사용합니다. 필터 차수보다 길이가 한 샘플 더 크고 β=3인 카이저 윈도우를 사용합니다. 카이저 윈도우에 대한 자세한 내용은 kaiser를 참조하십시오.

fir1을 사용하여 필터를 설계합니다. fir1은 구간 [0,1]에서 정규화 주파수를 필요로 합니다. 여기서 1은 π rad/sample에 해당합니다. fir1을 사용하려면 모든 주파수 사양을 정규화 주파수로 변환해야 합니다.

필터를 설계하고 크기 응답을 확인합니다.

fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,'low',kaiser(21,3));

[h,f] = freqz(b,1,[],Fs);
plot(f,mag2db(abs(h)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Magnitude (dB) contains an object of type line.

필터를 신호에 적용하고 100Hz 정현파의 처음 10개 주기에 대한 결과를 플로팅합니다.

y = filter(b,1,x);

plot(t,x,t,y)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 2 objects of type line. These objects represent Original Signal, Filtered Data.

designfilt를 사용하여 동일한 필터를 설계합니다. 필터 응답을 'lowpassfir'로 설정하고 사양을 Name,Value 쌍으로 입력합니다. designfilt를 사용하여 필터 설계의 단위를 Hz로 지정할 수 있습니다.

Fs = 1000;
Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150, ...
       'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

데이터를 필터링하고 결과를 플로팅합니다.

y1 = filter(Hd,x);

plot(t,x,t,y1)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 2 objects of type line. These objects represent Original Signal, Filtered Data.

필터 디자이너를 사용하여 저역통과 FIR 필터 설계

이 예제에서는 대화형 방식의 필터 디자이너 앱과 함께 윈도우 방법을 사용하여 저역통과 FIR 필터를 설계하고 구현하는 방법을 보여줍니다.

  • 명령줄에 filterDesigner를 입력하여 앱을 시작합니다.

  • 응답 유형저역통과로 설정합니다.

  • 설계 방법FIR로 설정하고 윈도우 방법을 선택합니다.

  • 필터 차수에서 차수 지정을 선택합니다. 차수를 20으로 설정합니다.

  • 주파수 사양에서 단위Hz로, Fs를 1000으로, Fc를 150으로 설정합니다.

  • 필터 설계를 클릭합니다.

  • 파일 > 내보내기...를 선택하여 FIR 필터를 MATLAB® 작업 공간에 계수나 필터 객체로 내보냅니다. 이 예제에서는 필터를 객체로 내보냅니다. 변수 이름을 Hd로 지정합니다.

  • 내보내기를 클릭합니다.

  • 명령 창에서 내보낸 필터 객체를 사용하여 입력 신호를 필터링합니다. 100Hz 정현파의 처음 10개 주기에 대한 결과를 플로팅합니다.

y2 = filter(Hd,x);

plot(t,x,t,y2)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 2 objects of type line. These objects represent Original Signal, Filtered Data.

  • 파일 > MATLAB 코드 생성 > 필터 설계 함수를 선택하여 여기서 사용한 사양으로 필터 객체를 만드는 MATLAB 함수를 생성합니다.

또한 대화형 방식의 툴 filterBuilder를 사용하여 필터를 설계할 수도 있습니다.

대역통과 필터 – 최소 차수 FIR 및 IIR 시스템

이 예제에서는 대역통과 필터를 설계하고 최소 차수 FIR 등리플 필터와 IIR 버터워스 필터를 사용하여 데이터를 필터링하는 방법을 보여줍니다. 많은 실제 신호를 진동 성분과 저주파수 추세, 가산성 잡음이 중첩된 형태로 모델링할 수 있습니다. 예를 들어, 경제 데이터에는 천천히 변화하는 상승 또는 하강 추세에 중첩된 주기를 나타내는 진동이 포함되는 경우가 많습니다. 또한, 데이터 수집 중에 발생하는 측정 오차 및 태생적으로 존재하는 임의 변동이 결합된 가산성 잡음 성분이 있을 수도 있습니다.

이 예제에서는 1년 동안 매일 어떠한 과정을 샘플링한다고 가정하겠습니다. 이 과정은 약 1주 및 1개월 범위에서 진동이 나타난다고 가정하겠습니다. 또한, 데이터와 가산성 N(0,1/4) 백색 가우스 잡음(AWGN)에 저주파수 상승 추세가 있습니다.

주파수가 각각 1/7 주기/일 및 1/30 주기/일인 두 사인파가 중첩된 형태를 띤 신호를 생성합니다. 저주파수가 증가하는 추세 항과 N(0,1/4) 백색 가우스 잡음을 더합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 데이터는 1샘플/일로 샘플링합니다. 결과로 생성되는 신호와 파워 스펙트럼 밀도(PSD) 추정값을 플로팅합니다.

rng default

Fs = 1;
n = 1:365;

x = cos(2*pi*(1/7)*n)+cos(2*pi*(1/30)*n-pi/4);
trend = 3*sin(2*pi*(1/1480)*n);

y = x+trend+0.5*randn(size(n));

[pxx,f] = periodogram(y,[],[],Fs);

subplot(2,1,1)
plot(n,y)
xlim([1 365])
xlabel('Days')
grid

subplot(2,1,2)
plot(f,10*log10(pxx))
xlabel('Cycles/day')
ylabel('dB')
grid

Figure contains 2 axes objects. Axes object 1 with xlabel Days contains an object of type line. Axes object 2 with xlabel Cycles/day, ylabel dB contains an object of type line.

저주파수 추세가 파워 스펙트럼 밀도 추정값에 증가된 저주파수 전력으로 나타납니다. 저주파수 전력은 1/30 주기/일의 진동보다 대략 10dB 높게 나타납니다. 필터 저지대역에 대한 사양에 이 정보를 사용합니다.

통과대역이 [1/40,1/4]주기/일부터이고 저지대역이 각각 [0,1/60]주기/일과 [1/4,1/2]주기/일부터인 사양을 사용하여 최소 차수 FIR 등리플 필터와 IIR 버터워스 필터를 설계합니다. 두 저지대역 감쇠량을 모두 10dB로 설정하고 통과대역 리플 허용오차를 1dB로 설정합니다.

Hd1 = designfilt('bandpassfir', ...
    'StopbandFrequency1',1/60,'PassbandFrequency1',1/40, ...
    'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 , ...
    'StopbandAttenuation1',10,'PassbandRipple',1, ...
    'StopbandAttenuation2',10,'DesignMethod','equiripple','SampleRate',Fs);
Hd2 = designfilt('bandpassiir', ...
    'StopbandFrequency1',1/60,'PassbandFrequency1',1/40, ...
    'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 , ...
    'StopbandAttenuation1',10,'PassbandRipple',1, ...
    'StopbandAttenuation2',10,'DesignMethod','butter','SampleRate',Fs);

FIR 필터 및 IIR 필터의 차수와 펼쳐진 위상 응답을 비교합니다.

fprintf('The order of the FIR filter is %d\n',filtord(Hd1))
The order of the FIR filter is 78
fprintf('The order of the IIR filter is %d\n',filtord(Hd2))
The order of the IIR filter is 8
[phifir,w] = phasez(Hd1,[],1);
[phiiir,w] = phasez(Hd2,[],1);

figure
plot(w,unwrap(phifir))
hold on
plot(w,unwrap(phiiir))
hold off

xlabel('Cycles/Day')
ylabel('Radians')
legend('FIR Equiripple Filter','IIR Butterworth Filter')
grid

Figure contains an axes object. The axes object with xlabel Cycles/Day, ylabel Radians contains 2 objects of type line. These objects represent FIR Equiripple Filter, IIR Butterworth Filter.

IIR 필터는 FIR 필터보다 훨씬 더 낮은 차수를 가집니다. 그러나, FIR 필터는 통과대역에서 선형 위상 응답을 갖는 반면, IIR 필터는 그렇지 않습니다. FIR 필터는 필터 통과대역의 모든 주파수를 균등하게 지연시키는 반면, IIR 필터는 그렇지 않습니다.

또한, 주파수 단위당 위상 변화율이 IIR 필터보다 FIR 필터에서 더 큽니다.

비교를 위해 저역통과 FIR 등리플 필터를 설계합니다. 저역통과 필터 사양의 경우, 통과대역은 [0,1/4]주기/일, 저지대역 감쇠량은 10dB, 통과대역 리플 허용오차는 1dB로 설정합니다.

Hdlow = designfilt('lowpassfir', ...
    'PassbandFrequency',1/4,'StopbandFrequency',1/2, ...
    'PassbandRipple',1,'StopbandAttenuation',10, ...
    'DesignMethod','equiripple','SampleRate',1);

대역통과 필터와 저역통과 필터를 사용하여 데이터를 필터링합니다.

yfir = filter(Hd1,y);
yiir = filter(Hd2,y);
ylow = filter(Hdlow,y);

대역통과 IIR 필터 출력값에 대한 PSD 추정값을 플로팅합니다. 다음 코드에서 yiiryfir로 대체하면 FIR 대역통과 필터 출력값에 대한 PSD 추정값을 볼 수 있습니다.

[pxx,f] = periodogram(yiir,[],[],Fs);

plot(f,10*log10(pxx))

xlabel('Cycles/day')
ylabel('dB')
grid

Figure contains an axes object. The axes object with xlabel Cycles/day, ylabel dB contains an object of type line.

PSD 추정값을 보면 대역통과 필터가 저주파수 추세와 고주파수 잡음을 감쇠한다는 것을 알 수 있습니다.

처음 120일 기간에 대한 FIR 및 IIR 필터 출력값을 플로팅합니다.

plot(n,yfir,n,yiir)

axis([1 120 -2.8 2.8])
xlabel('Days')
legend('FIR bandpass filter output','IIR bandpass filter output', ...
    'Location','SouthEast')

Figure contains an axes object. The axes object with xlabel Days contains 2 objects of type line. These objects represent FIR bandpass filter output, IIR bandpass filter output.

FIR 필터에서 증가된 위상 지연이 필터 출력값에 명확히 나타납니다.

비교를 위해 7일 주기와 30일 주기가 중첩된 형태 위에 겹쳐서 저역통과 FIR 필터 출력값을 플로팅합니다.

plot(n,x,n,ylow)

xlim([1 365])
xlabel('Days')
legend('7-day and 30-day cycles','FIR lowpass filter output', ...
    'Location','NorthWest')

Figure contains an axes object. The axes object with xlabel Days contains 2 objects of type line. These objects represent 7-day and 30-day cycles, FIR lowpass filter output.

앞의 플롯에 나와 있는 저역통과 필터 출력값에서 저주파수 추세가 명확히 표시되는 것을 확인할 수 있습니다. 저역통과 필터는 7일 주기와 30일 주기를 유지하지만, 대역통과 필터가 저주파수 추세도 제거하기 때문에 이 예제에서는 대역통과 필터의 성능이 더 좋습니다.

영위상 필터링

이 예제에서는 영위상 필터링을 수행하는 방법을 보여줍니다.

fir1designfilt를 사용하여 신호 생성 및 저역통과 필터 설계를 반복합니다. 작업 공간에 아래 변수가 이미 있는 경우에는 다음 코드를 실행할 필요가 없습니다.

rng default

Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));

% Using fir1
fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,'low',kaiser(21,3));

% Using designfilt
Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150, ...
       'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

filter를 사용하여 데이터를 필터링합니다. 진폭과 초기 위상이 입력 신호와 같은 처음 100개 점의 필터 출력값을 정현파와 함께 겹쳐 플로팅합니다.

yout = filter(Hd,x);
xin = cos(2*pi*100*t);

plot(t,xin,t,yout)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Input Sine Wave','Filtered Data')
grid

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 2 objects of type line. These objects represent Input Sine Wave, Filtered Data.

초기 0.01초 동안의 필터링된 데이터를 살펴보면 출력값이 입력값에 비해 지연되었음을 확인할 수 있습니다. 지연은 대략 0.01초로 보이며, 이는 샘플 (10×0.001)에서 FIR 필터 길이의 거의 절반에 해당합니다.

이 지연은 필터의 위상 응답으로 인해 발생합니다. 이 예제에 나와 있는 FIR 필터는 유형 I 선형 위상 필터입니다. 필터의 군지연은 10개 샘플입니다.

군지연을 플로팅합니다.

[gd,f] = grpdelay(Hd,[],Fs);

plot(f,gd)
xlabel('Frequency (Hz)')
ylabel('Group Delay (samples)')
grid

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel Group Delay (samples) contains an object of type line.

많은 응용 사례에서 위상 왜곡이 허용됩니다. 이는 위상 응답이 선형인 경우 특히 그렇습니다. 다른 응용 사례에서는 영위상 응답을 가진 필터를 사용하는 것이 바람직합니다. 영위상 응답은 기술적으로 비인과적 필터에서 가능하지 않습니다. filtfilt로 인과적 필터를 사용하는 영위상 필터링을 구현할 수 있습니다.

filtfilt를 사용하여 입력 신호를 필터링합니다. 응답을 플로팅하여 filterfiltfilt로 구한 필터 출력값을 비교합니다.

yzp = filtfilt(Hd,x);

plot(t,xin,t,yout,t,yzp)

xlim([0 0.1])
xlabel('Time (s)')
ylabel('Amplitude')
legend('100-Hz Sine Wave','Filtered Signal','Zero-phase Filtering',...
    'Location','NorthEast')

Figure contains an axes object. The axes object with xlabel Time (s), ylabel Amplitude contains 3 objects of type line. These objects represent 100-Hz Sine Wave, Filtered Signal, Zero-phase Filtering.

위 그림에서 filtfilt의 출력값이 FIR 필터의 위상 응답으로 인한 지연을 나타내지 않는 것을 확인할 수 있습니다.