filtfilt
영위상(Zero-Phase) 디지털 필터링
구문
설명
y = filtfilt(b,a,x)x를 순방향과 역방향 모두로 처리하여 영위상 디지털 필터링을 수행합니다. 이 함수는 데이터를 순방향으로 필터링한 후, 시작 부분과 끝부분의 과도를 최소화하기 위해 초기 조건과 일치하도록 만들고, 필터링된 시퀀스를 역방향으로 뒤집은 다음, 반전된 시퀀스를 다시 필터에 통과시킵니다. 결과는 다음과 같은 특성을 가집니다.
영위상 왜곡
원래 필터 전달 함수 크기의 제곱과 같은 필터 전달 함수
b와a로 지정된 필터 차수의 2배에 해당하는 필터 차수
filtfilt 함수는 Gustafsson[1]이 제안한 알고리즘을 구현합니다.
미분기 FIR 필터와 힐베르트 FIR 필터는 위상 응답에 따라 동작이 크게 달라지므로 이러한 필터와 함께 filtfilt 함수를 사용하지 마십시오.
y = filtfilt(d,x)d를 사용하여 입력 데이터 x를 영위상 필터링합니다. designfilt를 사용하여 주파수-응답 사양을 기반으로 d를 생성합니다.
예제
영위상 필터링을 사용하면 필터링되지 않은 신호에서 발생하는 특징들을, 필터링한 시간 파형의 정확히 동일한 위치에 유지할 수 있습니다.
합성 심전도(ECG) 파형을 영위상 필터링합니다. 파형을 생성하는 이 함수가 예제의 맨 끝에 있습니다. QRS 복합파는 심전도의 중요한 특징입니다. 여기서는 시간 지점 160 부근에서 QRS 복합파가 시작됩니다.
wform = ecg(500); plot(wform) axis([0 500 -1.25 1.25]) text(155,-0.4,"Q") text(180,1.1,"R") text(205,-1,"S")
가산성 잡음으로 심전도를 손상시킵니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 저역통과 FIR 등리플 필터를 생성하고 영위상 필터링과 일반적인 필터링을 둘 다 사용하여 잡음이 있는 파형을 필터링합니다.
rng("default") x = wform' + 0.25*randn(500,1); d = designfilt("lowpassfir", ... PassbandFrequency=0.15,StopbandFrequency=0.2, ... PassbandRipple=1,StopbandAttenuation=60, ... DesignMethod="equiripple"); y = filtfilt(d,x); y1 = filter(d,x); tiledlayout("flow") nexttile plot([y y1]) title("Filtered Waveforms") legend(["Zero-phase Filtering" "Conventional Filtering"]) nexttile plot(wform) title("Original Waveform")
영위상 필터링은 신호에서 잡음을 줄이고, QRS 복합파가 원래 신호에서 발생하는 시점과 동일한 시간으로 QRS 복합파를 유지합니다. 일반적인 필터링은 신호에서 잡음을 줄이지만, QRS 복합파를 지연시킵니다.
버터워스 2차섹션형(SOS) 필터를 사용하여 필터링을 반복합니다.
d1 = designfilt("lowpassiir",FilterOrder=12, ... HalfPowerFrequency=0.15,DesignMethod="butter"); y = filtfilt(d1,x); figure plot(x) hold on plot(y,LineWidth=3) hold off legend(["Noisy ECG" "Zero-Phase Filtering"])
다음 함수는 심전도 파형을 생성합니다.
function x = ecg(L) %ECG Electrocardiogram (ECG) signal generator. % ECG(L) generates a piecewise linear ECG signal of length L. % % EXAMPLE: % x = ecg(500).'; % y = sgolayfilt(x,0,3); % Typical values are: d=0 and F=3,5,9, etc. % y5 = sgolayfilt(x,0,5); % y15 = sgolayfilt(x,0,15); % plot(1:length(x),[x y y5 y15]); % Copyright 1988-2002 The MathWorks, Inc. a0 = [0,1,40,1,0,-34,118,-99,0,2,21,2,0,0,0]; % Template d0 = [0,27,59,91,131,141,163,185,195,275,307,339,357,390,440]; a = a0/max(a0); d = round(d0*L/d0(15)); % Scale them to fit in length L d(15)=L; for i=1:14 m = d(i):d(i+1)-1; slope = (a(i+1)-a(i))/(d(i+1)-d(i)); x(m+1) = a(i)+slope*(m-d(i)); end end
R2024b 이후
종속 연결 전달 함수를 사용하여 영위상 필터링을 수행합니다.
통과대역 리플과 저지대역 감쇠량이 각각 0.1dB와 40dB인 타원 필터를 설계합니다. 샘플 레이트를 2000Hz로 지정합니다. 필터의 주파수 응답을 플로팅합니다.
Fs = 2000; [B,A] = ellip(20,0.1,40,[0.3 0.7],"ctf"); freqz(B,A,2048,Fs,"ctf")
나이퀴스트 주파수가 1초에 발생하는 선형 스윕 처프 신호를 필터링합니다. 입력 신호와 출력 신호 사이의 스펙트럼을 비교합니다.
t = 0:1/Fs:1;
x = chirp(t,0,t(end),Fs/2)';
y = filtfilt(B,A,x,"ctf");
pspectrum([x y],Fs,Leakage=1,FrequencyResolution=1)
R2024b 이후
전달 함수 영위상 필터링을 사용하여 잡음이 있는 정현파 아티팩트를 다시 만듭니다. CTF와 스케일 값을 사용하여 진동 신호를 필터링합니다.
정규 분포 잡음과 세 개의 정현파 파형으로 구성된 신호 u를 생성합니다. 샘플 레이트는 1kHz입니다.
rng("default")
Fs = 1e3;
ts = (0:1/Fs:2)';
a0 = [3 2 1];
f0 = [0.1 0.5 0.9]*Fs/2;
p0 = [0 pi/4 pi/2];
u = 0.1*randn(size(ts)) + 0.1*sin(2*pi*f0.*ts+p0)*a0';3차 버터워스 대역저지 디지털 필터로 n0을 필터링하여 잡음이 있는 정현파 아티팩트를 다시 만들고 신호 v를 생성합니다.
[b,a] = butter(3,[0.15 0.85],"stop");
v = filtfilt(b,a,u);u와 v를 비교합니다. 두 신호의 위상이 일치하는 것을 확인합니다.
tiledlayout("flow") nexttile strips([u(ts<0.1) v(ts<0.1)],0.1,Fs) legend(["u" "v"],Location="southeast") xlabel("Time (seconds)") nexttile pspectrum([u v],Fs) legend(["u" "v"],Location="southeast")
전압 제어 진동 신호 x를 만듭니다. 신호 v로 표현되는 잡음이 있는 정현파 아티팩트를 추가합니다.
vo = exp(-2*abs(ts-1)).*sin(8*pi*ts); x = vco(vo,[0.25 0.75]*Fs/2,Fs) + v;
24차 체비쇼프 유형 II 필터로 신호 x를 필터링합니다. CTF 형식과 스케일 값(B,A,g)을 사용합니다.
[B,A,g] = cheby2(24,50,[0.2 0.8],"ctf");
y = filtfilt({B,A,g},x);단시간 푸리에 변환의 크기 제곱을 비교합니다. 저지대역에서 크기가 급격히 감소하는 것을 확인합니다.
tiledlayout("flow") nexttile stft(x,Fs,Window=bohmanwin(128),OverlapLength=120, ... FFTLength=512,FrequencyRange="onesided") title("Input x") nexttile stft(y,Fs,Window=bohmanwin(128),OverlapLength=120, ... FFTLength=512,FrequencyRange="onesided") title("Output y")
입력 인수
출력 인수
세부 정보
팁
참고 문헌
[1] Gustafsson, F. “Determining the initial states in forward-backward filtering.” IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992. https://doi.org/10.1109/78.492552.
[2] Lyons, Richard G. Understanding Digital Signal Processing. Upper Saddle River, NJ: Prentice Hall, 2004.
[3] Mitra, Sanjit K. Digital Signal Processing. 2nd Ed. New York: McGraw-Hill, 2001.
[4] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.
