이 페이지의 최신 내용은 아직 번역되지 않았습니다. 최신 내용은 영문으로 볼 수 있습니다.

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

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

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

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

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

디지털 필터는 신호에 지연을 일으킵니다. 필터 특성에 따라 지연은 모든 주파수에 걸쳐 일정할 수도 있고, 주파수에 따라 다를 수도 있습니다. 지연의 유형에 따라 이를 보정하는 데 수행해야 하는 작업이 결정됩니다. 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

% Design a 70th order lowpass FIR filter with cutoff frequency of 75 Hz.

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

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

grpdelay(df,2048,Fs)   % plot group delay
D = mean(grpdelay(df)) % filter delay in samples
D =

    35

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

y = filter(df,[x; zeros(D,1)]); % Append D zeros to the input data
y = y(D+1:end);                  % Shift data to compensate for delay

figure
plot(t,x,t,y,'r','linewidth',1.5);
title('Filtered Waveforms');
xlabel('Time (s)')
legend('Original Noisy Signal','Filtered Signal');
grid on
axis tight

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

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

이전 섹션에서 정의한 심전도 신호를 사용해 보겠습니다. 이 신호를 지연 보정을 적용한 상태와 적용하지 않은 상태 모두에서 필터링합니다.

% Design a 7th order lowpass IIR elliptic filter with cutoff frequency
% of 75 Hz.

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

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

grpdelay(df,2048,'half',Fs)

데이터를 필터링하고 각 필터 구현이 시간 신호에 미치는 영향을 검토합니다.

y1 = filter(df,x);    % non-linear phase filter - no delay compensation
y2 = filtfilt(df,x);  % zero-phase implementation - delay compensation

figure
plot(t,x);
hold on
plot(t,y1,'r','linewidth',1.5);
plot(t,y2,'g','linewidth',1.5);
title('Filtered Waveforms');
xlabel('Time (s)')
legend('Original Signal','Non-linear phase IIR output',...
  'Zero-phase IIR output');
ax = axis;
axis([0.25 0.55 ax(3:4)])
grid on

영위상 필터링(Zero-phase Filtering)이 얼마나 효과적으로 필터 지연을 제거하는지 알 수 있습니다.

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

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

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

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

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

오디오 신호를 불러옵니다.

Fs = 44100; % Sample rate
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'})

먼저, 저역통과 필터를 사용하여 가능한 한 많이 백색 잡음 스펙트럼 성분을 제거할 수 있습니다. 필터의 통과대역은 잡음 감소와 고주파 성분 손실로 인한 음질 저하 간에 절충된 값으로 설정되어야 합니다. 대역 제한 신호를 다운샘플링할 수 있기 때문에 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

% Design the filter
df = designfilt('lowpassfir','PassbandFrequency',Fp,...
                'StopbandFrequency',Fst,'PassbandRipple',Ap,...
                'StopbandAttenuation',Ast,'SampleRate',Fs);

% Analyze the filter response
hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',F);

% Filter the data and compensate for delay
D = mean(grpdelay(df)); % filter delay
ylp = filter(df,[y; zeros(D,1)]);
ylp = ylp(D+1:end);

close(hfvt)

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

[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,Flp,Plp,...
  {'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'})

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

% Design the filter
df = designfilt('bandstopiir','PassbandFrequency1',55,...
               'StopbandFrequency1',58,'StopbandFrequency2',62,...
               'PassbandFrequency2',65,'PassbandRipple1',1,...
               'StopbandAttenuation',60,'PassbandRipple2',1,...
               'SampleRate',Fs,'DesignMethod','ellip');

% Analyze the magnitude response
hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',Fds(Fds>F(2)));

왜곡을 피하기 위해 영위상 필터링을 수행합니다.

ybs = filtfilt(df,yds);

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

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

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

[Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,'power');
close(hfvt)
helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal,...
  {'Original signal','Final filtered signal'})

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

% Play the original signal
hplayer = audioplayer(y, Fs);
play(hplayer);

% Play the noise-reduced signal
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로 설정합니다.

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

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

hfvt = fvtool(df,[1 -1],1,'MagnitudeDisplay','zero-phase','Fs',Fs);
legend(hfvt,'50th order FIR differentiator','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(df)); % filter delay
v2 = filter(df,[drift; zeros(D,1)]);
v2 = v2(D+1:end);
a2 = filter(df,[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)

신호 적분하기

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

인 누설 적분기를 사용합니다. 누설 적분기 필터의 크기 응답을 플로팅합니다. 필터가 저역통과 필터로 동작하여 고주파 잡음을 제거하는 것을 알 수 있습니다.

close(hfvt)
fvtool(1,[1 -.999],'Fs',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);

% Multiply by time differential
d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;

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

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

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

결론

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

추가 참고 자료

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

참고 문헌: J.G. Proakis and D. G. Manolakis, "Digital Signal Processing. Principles, Algorithms, and Applications", Prentice-Hall, 1996.

S. J. Orfanidis, "Introduction To Signal Processing", Prentice-Hall, 1996.

부록

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