Main Content

디지털 필터링에 대한 실용적 입문

이 예제에서는 디지털 필터를 설계 및 분석하여 데이터에 적용하는 방법을 보여줍니다. 이는 다음과 같은 질문에 답하는 데 도움이 됩니다. 필터 사용으로 인해 발생하는 지연을 보정하려면 어떻게 해야 합니까?, 신호가 왜곡되는 것을 방지하려면 어떻게 해야 합니까?, 신호에서 원치 않는 성분을 제거하려면 어떻게 해야 합니까?, 신호를 미분하려면 어떻게 해야 합니까?, 신호를 적분하려면 어떻게 해야 합니까?

필터는 신호 스펙트럼의 형태를 원하는 방식으로 만들거나 미분 및 적분과 같은 수학 연산을 수행하는 데 사용할 수 있습니다. 다음에서는 필터를 필요에 따라 쉽게 사용할 수 있도록 몇 가지 실용적인 개념을 설명합니다.

이 예제에서는 디지털 필터의 설계보다 디지털 필터의 응용에 대해 중점적으로 다룹니다. 디지털 필터를 설계하는 방법에 대한 자세한 내용을 보려면 디지털 필터 설계에 대한 실용적 입문 예제를 참조하십시오.

필터링으로 인해 발생하는 지연 보정하기

디지털 필터는 신호에 지연을 일으킵니다. 필터 특성에 따라 지연은 모든 주파수에 걸쳐 일정할 수도 있고, 주파수에 따라 다를 수도 있습니다. 지연의 유형에 따라 이를 보정하는 데 수행해야 하는 작업이 결정됩니다. grpdelay 함수를 사용하면 필터 지연을 주파수 함수로 간주할 수 있습니다. 이 함수의 출력값을 보면 필터의 지연이 일정한지 아니면 주파수에 따라 달라지는지(즉, 주파수 종속 지연인지) 식별할 수 있습니다.

모든 주파수에 걸쳐 일정한 필터 지연은 시간상에서 신호를 이동하는 방법으로 쉽게 보정할 수 있습니다. FIR 필터는 일반적으로 일정한 지연을 가집니다. 반면에 ,주파수에 따라 달라지는 지연은 위상 왜곡을 초래하며 신호 파형을 크게 변형시킬 수 있습니다. 주파수 종속 지연에 대한 보정은 일정한 지연에 대한 경우만큼 간단하지 않습니다. IIR 필터는 주파수 종속 지연을 일으킵니다.

일정한 필터 지연에 대한 보정

앞에서 설명한 것처럼 필터의 군지연을 측정하여 지연이 주파수의 상수 함수인지 확인할 수 있습니다. grpdelay 함수를 사용하면 필터 지연 D를 측정하고 입력 신호에 D개의 0을 추가하고 시간상에서 D개 샘플만큼 출력 신호를 이동하여 이 지연을 보정할 수 있습니다.

잡음이 포함된 심전도 신호에서 75Hz보다 높은 고주파 잡음을 제거하도록 필터링하는 경우를 가정하겠습니다. 잡음이 포함된 신호와 필터링된 신호가 올바르게 정렬되고 비교를 위해 서로 겹쳐칠 수 있도록 FIR 저역통과 필터를 적용하고 필터 지연을 보정합니다.

Fs = 500;                    % Sample rate in Hz
N = 500;                     % Number of signal samples
rng default;
x = ecg(N)'+0.25*randn(N,1); % Noisy waveform
t = (0:N-1)/Fs;              % Time vector

75Hz의 차단 주파수를 갖는 70차 저역통과 FIR 필터를 설계합니다.

Fnorm = 75/(Fs/2);           % Normalized frequency
df3 = designfilt("lowpassfir",FilterOrder=70,CutoffFrequency=Fnorm);

필터의 군지연을 플로팅하여 지연이 모든 주파수에 걸쳐 일정한지 확인합니다. 이 경우 필터는 선형 위상입니다. 군지연을 사용하여 필터의 지연을 측정합니다.

grpdelay(df3,2048,Fs)         % Plot group delay

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

D = mean(grpdelay(df3)) % Filter delay in samples
D = 35

필터링하기 전에 입력 데이터 벡터 x의 끝에 D개의 0을 추가합니다. 그러면 모든 유용한 샘플이 필터에서 제거되고 입력 신호와 지연 보정된 출력 신호의 길이가 같게 됩니다. 데이터를 필터링하고 출력 신호를 D개 샘플만큼 이동하여 지연을 보정합니다. 이 마지막 단계는 필터의 과도 특성을 효과적으로 제거합니다.

y = filter(df3,[x; zeros(D,1)]);  % Append D zeros to the input data
y = y(D+1:end);                  % Shift data to compensate for delay
plot(t,x,t,y,LineWidth=1.5)
title("Filtered Waveforms")
xlabel("Time (s)")
legend("Original Noisy Signal","Filtered Signal")
grid on
axis tight

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

주파수 종속 지연에 대한 보정

주파수 종속 지연은 신호에서 위상 왜곡을 일으킵니다. 이 유형의 지연에 대한 보정은 일정한 지연에 대한 보정만큼 간단하지 않습니다. 애플리케이션에서 오프라인 처리를 허용하는 경우 filtfilt 함수를 사용하여 영위상 필터링을 구현하는 방식으로 주파수 종속 지연을 제거할 수 있습니다. filtfilt는 순방향과 역방향 모두로 입력 데이터를 처리하여 영위상 필터링을 수행합니다. 주된 효과는 영위상 왜곡을 얻는 것입니다. 즉, 위상이 0인 상수 지연을 갖는 샘플 데이터를 필터로 필터링(제거)하는 것과 마찬가지입니다. 다른 효과는 원래 필터 전달 함수의 크기 제곱에 해당하는 필터 전달 함수와 원래 필터 차수의 2배에 해당하는 필터 차수를 얻는 것입니다.

이전 섹션에서 정의한 심전도 신호를 사용해 보겠습니다. 이 신호를 지연 보정을 적용한 상태와 적용하지 않은 상태 모두에서 필터링합니다. 75Hz의 차단 주파수를 갖는 7차 저역통과 IIR 타원 필터를 설계합니다.

Fnorm = 75/(Fs/2); % Normalized frequency
df4 = designfilt("lowpassiir", ...
    FilterOrder=7, ...
    PassbandFrequency=Fnorm, ...
    PassbandRipple=1, ...
    StopbandAttenuation=60);

필터의 군지연을 플로팅합니다. 군지연이 주파수에 따라 달라지는 것을 알 수 있습니다. 이는 필터 지연이 주파수에 종속적임을 나타냅니다.

grpdelay(df4,2048,Fs)

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

데이터를 필터링하고 각 필터 구현이 시간 신호에 미치는 영향을 검토합니다. 영위상 필터링(Zero-phase Filtering)이 효과적으로 필터 지연을 제거합니다.

y1 = filter(df4,x);    % Nonlinear phase filter - no delay compensation
y2 = filtfilt(df4,x);  % Zero-phase implementation - delay compensation
plot(t,x)
hold on
plot(t,y1,"r",LineWidth=1.5)
plot(t,y2,LineWidth=1.5)
title("Filtered Waveforms")
xlabel("Time (s)")
legend("Original Signal","Nonlinear Phase IIR Output", ...
  "Zero-Phase IIR Output")
xlim([0.25 0.55])
grid on

Figure contains an axes object. The axes object with title Filtered Waveforms, xlabel Time (s) contains 3 objects of type line. These objects represent Original Signal, Nonlinear Phase IIR Output, Zero-Phase IIR Output.

영위상 필터링은 애플리케이션이 비인과적 순방향/역방향 필터링 작업을 허용하는 경우와 원래 응답의 제곱으로 필터 응답을 변경하는 경우 사용할 수 있는 훌륭한 툴입니다.

상수 지연을 일으키는 필터는 선형 위상 필터입니다. 주파수 종속 지연을 일으키는 필터는 비선형 위상 필터입니다.

신호에서 원치 않는 스펙트럼 성분 제거하기

필터는 일반적으로 신호에서 원치 않는 스펙트럼 성분을 제거하는 데 사용됩니다. 다양한 필터를 선택하여 이러한 작업을 수행할 수 있습니다. 고주파 성분을 제거하려는 경우 저역통과 필터(Lowpass Filter)를 선택하고, 저주파 성분을 제거하려는 경우 고역통과 필터(Highpass Filter)를 선택합니다. 또한, 주파수의 중간 대역은 그대로 두고 저주파 성분과 고주파 성분만 제거하려면 대역통과 필터(Bandpass Filter)를 선택하면 됩니다. 지정된 대역에 있는 주파수를 제거하려면 대역저지 필터(Bandstop Filter)를 선택합니다.

전력 공급선 험(Hum)과 백색 잡음이 있는 오디오 신호가 있다고 가정하겠습니다. 전력 공급선의 험 잡음은 60Hz 톤으로 인해 발생합니다. 백색 잡음은 모든 오디오 대역폭에 걸쳐 존재하는 신호입니다.

오디오 신호를 불러옵니다. 샘플 레이트를 44.1kHz로 지정합니다.

Fs = 44100;
y = audioread("noisymusic.wav");

신호의 파워 스펙트럼을 플로팅합니다. 빨간색 삼각형 마커는 오디오 신호를 간섭하는 강력한 60Hz 톤을 표시합니다.

[P,F] = pwelch(y,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365], ...
  {"Original Signal Power Spectrum", "60 Hz Tone"})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Original Signal Power Spectrum, 60 Hz Tone.

먼저, 저역통과 필터를 사용하여 가능한 한 많이 백색 잡음 스펙트럼 성분을 제거할 수 있습니다. 필터의 통과대역은 잡음 감소와 고주파 성분 손실로 인한 음질 저하 간에 절충된 값으로 설정되어야 합니다. 대역 제한 신호를 다운샘플링할 수 있기 때문에 60Hz 험을 제거하기 전에 저역통과 필터를 적용하면 매우 간편합니다. 신호 레이트가 낮을수록 더 작은 필터 차수로 더 예리하고 더 좁은 60Hz 대역저지 필터를 설계할 수 있습니다.

통과대역 주파수가 1kHz이고 저지대역 주파수가 1.4kHz인 저역통과 필터를 설계합니다. 최소 차수 설계를 선택합니다.

Fp = 1e3;    % Passband frequency in Hz
Fst = 1.4e3; % Stopband frequency in Hz
Ap = 1;      % Passband ripple in dB
Ast = 95;    % Stopband attenuation in dB
df3 = designfilt("lowpassfir",PassbandFrequency=Fp, ...
    StopbandFrequency=Fst,PassbandRipple=Ap, ...
    StopbandAttenuation=Ast,SampleRate=Fs);

필터 응답을 분석합니다.

opts = filterAnalysisOptions;
opts.FrequencyRange = "freqvector";
opts.FrequencyVectorUnits = "Hz";
opts.FrequencyVector = F;
opts.FrequencyScale = "log";
filterAnalyzer(df3,AnalysisOptions=opts)

데이터를 필터링하고 지연을 보정합니다.

D = mean(grpdelay(df3)); % Filter delay
ylp = filter(df3,[y; zeros(D,1)]);
ylp = ylp(D+1:end);

저역통과 필터가 적용된 신호의 스펙트럼을 확인합니다. 1400Hz보다 높은 주파수 성분이 제거되었습니다.

[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Flp,Plp, ...
    {"Original Signal","Lowpass Filtered Signal"})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Original Signal, Lowpass Filtered Signal.

위에 나와 있는 파워 스펙트럼 플롯에서 저역통과 필터가 적용된 신호의 무시할 수 없는 최대 주파수 성분이 1400Hz에 있는 것을 알 수 있습니다. 샘플링 정리에 따라 샘플 주파수 2×1400=2800Hz는 신호를 올바르게 나타내는 데 충분합니다. 그러나 여기서는 44100Hz의 샘플 레이트를 사용합니다. 이는 필요한 샘플보다 더 많은 샘플을 처리해야 하므로 낭비입니다. 신호를 다운샘플링하여 처리해야 하는 샘플 수를 줄이는 방식으로 샘플 레이트를 줄이고 계산 부하를 줄일 수 있습니다. 또한, 샘플 레이트가 낮을수록 더 작은 필터 차수로, 60Hz 잡음을 제거하는 데 필요한 더 예리하고 더 좁은 대역저지 필터를 설계할 수 있습니다.

인자 10으로 저역통과 필터가 적용된 신호를 다운샘플링하여 Fs/10 = 4.41kHz의 샘플 레이트를 얻습니다. 다운샘플링을 수행하기 전과 후에 신호의 스펙트럼을 플로팅합니다.

Fs = Fs/10;
yds = downsample(ylp,10);
[Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Fds,Pds, ...
    {"Signal Sampled at 44100 Hz", "Downsampled Signal, Fs = 4410 Hz"})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Signal Sampled at 44100 Hz, Downsampled Signal, Fs = 4410 Hz.

이제, IIR 대역저지 필터를 사용하여 60Hz 톤을 제거합니다. 저지대역이 중심 60Hz, 폭 4Hz가 되도록 합니다. 예리한 주파수 노치, 작은 통과대역 리플, 상대적으로 낮은 차수를 실현하기 위해 IIR 필터를 선택합니다.

df4 = designfilt("bandstopiir", ...
    PassbandFrequency1=55, ...
    StopbandFrequency1=58, ...
    StopbandFrequency2=62, ...
    PassbandFrequency2=65, ...
    PassbandRipple1=1, ...
    StopbandAttenuation=60, ...
    PassbandRipple2=1, ...
    DesignMethod="ellip", ...
    SampleRate=Fs);

크기 응답을 분석합니다.

opts = filterAnalysisOptions;
opts.FrequencyRange = "freqvector";
opts.FrequencyVectorUnits = "Hz";
opts.FrequencyVector = Fds(Fds>F(2));
opts.FrequencyScale = "log";
filterAnalyzer(df4,AnalysisOptions=opts)

위상 왜곡을 피하기 위해 영위상 필터링을 수행합니다. filtfilt 함수를 사용하여 데이터를 처리합니다.

ybs = filtfilt(df4,yds);

마지막으로, 신호를 업샘플링하여 오디오 사운드 카드와 호환되는 44.1kHz의 원래 오디오 샘플 레이트로 다시 되돌립니다.

yf = interp(ybs,10);
Fs = Fs*10;

원래 신호와 처리된 신호의 스펙트럼을 최종적으로 살펴봅니다. 고주파 바닥 잡음과 60Hz 톤이 필터에 의해 감쇠되었습니다.

[Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,"power");
helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal, ...
    {"Original Signal","Final Filtered Signal"})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Original Signal, Final Filtered Signal.

컴퓨터에 사운드 카드가 있다면 처리를 수행하기 전과 후의 신호를 들어볼 수 있습니다. 위에 설명된 대로, 최종적으로 오디오 파일에서 60Hz 험과 고주파 잡음이 효과적으로 감쇠되었습니다.

% To play the original signal, uncomment the next two lines
% hplayer = audioplayer(y, Fs);
% play(hplayer)

% To play the noise-reduced signal, uncomment the next two lines
% hplayer = audioplayer(yf, Fs);
% play(hplayer);

신호 미분하기

MATLAB diff 함수는 출력값에서 잠재적으로 잡음 수준을 늘릴 수 있는 문제를 안고 신호를 미분합니다. 더 좋은 방법은 관심 대역에서 미분기로 동작하고, 그 밖의 모든 주파수에서 감쇠기로 동작하는 미분기 필터를 사용하여 고주파 잡음을 효과적으로 제거하는 것입니다.

이에 대한 예로 지진이 일어나는 동안 건물 바닥의 변위 속도를 분석해 보겠습니다. 지진 조건하에서 3층짜리 테스트 구조물의 1층에서 변위 또는 변동 측정값이 기록되어 quakedrift.mat 파일에 저장되었습니다. 데이터 벡터의 길이는 10e3이고, 샘플 레이트는 1kHz이며, 측정값 단위는 cm입니다.

변위 데이터를 미분하여 지진이 일어나는 동안 건물 바닥의 속도와 가속도에 대한 추정값을 구합니다. diff를 사용한 결과와 FIR 미분기 필터를 사용한 결과를 비교합니다.

load quakedrift.mat 
Fs  = 1000;                 % Sample rate
dt = 1/Fs;                  % Time differential
t = (0:length(drift)-1)*dt; % Time vector

신호 에너지의 대부분이 검출되는 대역폭인 100Hz의 통과대역 주파수를 갖는 50차 미분기 필터를 설계합니다. 필터의 저지대역 주파수를 120Hz로 설정합니다.

df5 = designfilt("differentiatorfir", ...
    FilterOrder=50, ...
    PassbandFrequency=100, ...
    StopbandFrequency=120, ...
    SampleRate=Fs);

diff 함수는 응답 H(Z)=1-Z-1을 갖는 1차 FIR 필터로 볼 수 있습니다. 필터 분석기를 사용하여 50차 미분기 FIR 필터의 크기 응답과 diff 함수의 응답을 비교합니다. 확실히, 두 응답은 통과대역 영역(0Hz~100Hz)에서 동일합니다. 그러나, 저지대역 영역에서 50차 필터는 성분을 감쇠하지만, diff 응답은 성분을 증폭시킵니다. 이는 고주파 잡음 수준을 늘리게 됩니다.

fa = filterAnalyzer(df5, ...
    MagnitudeMode="zerophase", ...
    FilterNames="FIR_Differentiator_50th_Order");
addFilters(fa,[1 -1],1, ...
    SampleRate=Fs, ...
    FilterNames="Response_of_Diff_Function");

diff 함수를 사용하여 미분합니다. 0을 추가하여 diff 연산으로 인해 누락된 샘플을 보완합니다.

v1 = diff(drift)/dt;
a1 = diff(v1)/dt;

v1 = [0; v1];
a1 = [0; 0; a1];

50차 FIR 필터를 사용하여 미분하고 지연에 대해 보정합니다.

D = mean(grpdelay(df5)); % Filter delay
v2 = filter(df5,[drift; zeros(D,1)]);
v2 = v2(D+1:end);
a2 = filter(df5,[v2; zeros(D,1)]);
a2 = a2(D+1:end);
v2 = v2/dt;
a2 = a2/dt^2;

바닥 변위에 대한 몇 개의 데이터 점을 플로팅합니다. 또한 diff로 계산된 경우와 50차 FIR 필터로 계산된 경우의 속도와 가속도에 대한 몇 개의 데이터 점도 플로팅합니다. 잡음이 diff로 구한 속도 추정값에서는 약간 증폭되고 가속도 추정값에서는 크게 증폭되는 것을 알 수 있습니다.

helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2)

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel Displacement (cm) contains an object of type line.

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel speed (cm/s) contains 2 objects of type line. These objects represent Estimated speed using diff, Estimated speed using FIR filter.

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel acceleration (cm/s Squared baseline ) contains 2 objects of type line. These objects represent Estimated acceleration using diff, Estimated acceleration using FIR filter.

신호 적분하기

누설 적분기 필터는 전달 함수 H(Z)=1/[1-cZ-1]을 사용하는 전극점 필터(All-pole Filter)이며, 여기서 c는 필터의 안정성을 보장하기 위해 1보다 작아야 하는 상수입니다. c가 1에 가까워질수록 누설 적분기가 diff 전달 함수의 역에 가까워진다는 것은 당연합니다. 이전 섹션에서 구한 가속도 추정값과 속도 추정값에 누설 적분기를 적용하여 각각 속도와 변동을 다시 얻습니다. diff 함수로 구한 추정값에 더 잡음이 많으므로 이 추정값을 사용합니다.

a=0.999인 누설 적분기를 사용합니다. 누설 적분기 필터의 크기 응답을 플로팅합니다. 필터가 저역통과 필터로 동작하여 고주파 잡음을 효과적으로 제거합니다.

filterAnalyzer(1,[1 -0.999],SampleRate=Fs)

누설 적분기로 속도와 가속도를 필터링합니다. 시간 미분소로 곱합니다.

v_original = v1;
a_original = a1;

d_leakyint = filter(1,[1 -0.999],v_original);
v_leakyint = filter(1,[1 -0.999],a_original);

d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;

변위 추정값과 속도 추정값을 플로팅하고 원래 신호와 비교합니다.

helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint)

Figure contains 2 axes objects. Axes object 1 with ylabel Displacement (cm) contains 2 objects of type line. These objects represent Leaky integrator estimate, Original displacement. Axes object 2 with xlabel Time in seconds, ylabel Speed (cm/s) contains 2 objects of type line. These objects represent Leaky integrator estimate, Original speed.

또한 cumsum 함수와 cumtrapz 함수를 사용하여 신호를 적분할 수도 있습니다. 결과는 누설 적분기로 구한 결과와 비슷합니다.

결론

이 예제에서는 선형 위상 필터와 비선형 위상 필터에 대해 배웠으며 각 필터 유형으로 인해 발생하는 위상 지연을 보정하는 방법에 대해 배웠습니다. 또한, 필터를 적용하여 신호에서 원치 않는 주파수 성분을 제거하는 방법과 저역통과 필터를 사용하여 신호의 대역폭을 제한한 후 신호를 다운샘플링하는 방법에 대해서도 배웠습니다. 마지막으로, 디지털 필터 설계를 사용하여 신호를 미분하고 적분하는 방법에 대해 배웠습니다. 또한 이 예제 전체에 걸쳐 분석 툴을 사용하여 필터의 응답과 군지연을 확인하는 방법에 대해서도 배웠습니다.

필터 응용에 대한 자세한 내용은 Signal Processing Toolbox™ 문서를 참조하십시오. 디지털 필터를 설계하는 방법에 대한 자세한 내용은 디지털 필터 설계에 대한 실용적 입문 예제를 참조하십시오.

참고 문헌

Proakis, J. G., and D. G. Manolakis. Digital Signal Processing: Principles, Algorithms, and Applications. Englewood Cliffs, NJ: Prentice-Hall, 1996.

Orfanidis, S. J. Introduction To Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1996.

부록

이 예제에서는 다음 헬퍼 함수가 사용되었습니다.

참고 항목

| | | | | | | |