이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.
filtfilt
R2026a영위상(Zero-Phase) 디지털 필터링
구문
설명
y = filtfilt(b,a,x)x를 순방향과 역방향 모두로 처리하여 영위상 디지털 필터링을 수행합니다. 이 함수는 데이터를 순방향으로 필터링한 후, 시작 부분과 끝부분의 과도를 최소화하기 위해 초기 조건과 일치하도록 만들고, 필터링된 시퀀스를 역방향으로 뒤집은 다음, 반전된 시퀀스를 다시 필터에 통과시킵니다. 결과는 다음과 같은 특성을 가집니다.
영위상 왜곡
원래 필터 전달 함수 크기의 제곱과 같은 필터 전달 함수
b와a로 지정된 필터 차수의 2배에 해당하는 필터 차수
미분기 FIR 필터와 힐베르트 FIR 필터는 위상 응답에 따라 동작이 크게 달라지므로 이러한 필터와 함께 filtfilt 함수를 사용하지 마십시오.
y = filtfilt(d,x)d를 사용하여 입력 데이터 x를 영위상 필터링합니다. designfilt를 사용하여 주파수-응답 사양을 기반으로 d를 생성합니다.
y = filtfilt(___,Name=Value)
예제
영위상 필터링을 사용하면 필터링되지 않은 신호에서 발생하는 특징들을, 필터링한 시간 파형의 정확히 동일한 위치에 유지할 수 있습니다.
합성 심전도(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")
R2026a 이후
60Hz 간격으로 최대 600Hz까지 확장되는 여러 고조파를 포함하는 복소수 값 30Hz 정현파 톤에 영위상 필터링을 적용합니다.
샘플 레이트가 600Hz인 2초 길이의 복소수 값 신호를 만듭니다. 이 신호는 10V 30Hz 사인 톤으로 구성되며, 60Hz에서 600Hz 사이까지 균등하게 분포된 10개의 고조파가 포함되어 있고 각 고조파의 진폭은 1V입니다.
Fs = 600; t = (0:1/Fs:2)'; x = 10*exp(1i*2*pi*30*t) + sum(exp(1i*2*pi*60*t*(0:9)),2);
톤이 30Hz에서 진동하고 원치 않는 성분이 60Hz의 배수로 나타나 최대 10번째 배수까지 주파수 피크를 가지므로, 관심 신호인 톤을 복원하는 데에는 10차 IIR 콤(Comb) 노치 필터가 적합합니다.
b와 a를 품질 계수가 35인 10차 IIR 콤(Comb) 노치 필터의 분자 계수와 분모 계수를 나타내는 벡터로 정의합니다. 특정 차수와 품질 계수를 가진 IIR 콤(Comb) 필터를 설계하는 방법에 대한 자세한 내용은 iircomb (DSP System Toolbox) 항목을 참조하십시오.
b = [0.957 zeros(1,9) -0.957]; a = [1 zeros(1,9) -0.914];
0.9초에서 1.1초까지 신호의 실수부를 플로팅합니다. 0Hz에서 600Hz까지 필터 응답을 플로팅합니다.
tiledlayout("flow") nexttile plot(t,real(x)) xlim([0.9 1.1]) xlabel("Time (s)") title("real(x)") nexttile [h,f] = freqz(b,a,8192,"whole",Fs); plot(f,mag2db(abs(h))) xlabel("Frequency (Hz)") title("Magnitude Response (b,a)")
위상을 유지하면서 입력 신호를 필터링합니다. Gustafsson 기법을 사용하여 필터 상태의 초기 조건을 추정합니다. 과도 길이가 신호 길이와 동일하다고 가정하고 샘플을 대칭 복사하여 입력 신호를 채웁니다.
y = filtfilt(b,a,x, ... InitialStatesMethod="gustafsson", ... TransientLength=length(x), ... PaddingPattern="reflect");
시간 영역과 주파수 영역에서 입력 신호와 필터링된 신호를 비교합니다. 0.9초에서 1.1초까지 두 신호의 실수부와 0Hz에서 600Hz까지 Welch 파워 스펙트럼을 플로팅합니다. 필터링된 신호는 고조파가 제거된 30Hz 톤을 보여주며, 필터링된 신호 y는 입력 신호 x와 위상이 일치합니다.
figure; tiledlayout("vertical") nexttile plot(t,real(x),t,real(y),"--") legend("real(" + ["x" "y"] +")",Location="southeast") xlabel("Time (s)") xlim([0.9 1.1]) nexttile [p,f] = pwelch([x y],[],[],8192,Fs); pdb = pow2db(abs(p)); plot(f,pdb(:,1),"-",f,pdb(:,2),"--") xlabel("Frequency (Hz)") legend("PSD of " + ["x" "y"],Location="southeast")
입력 인수
이름-값 인수
출력 인수
세부 정보
팁
스케일링 이득을 포함한 CTF 형식의 필터를 구할 수 있습니다. butter, cheby1, cheby2, ellip와 같은 디지털 IIR 필터 설계 함수의 출력값을 사용하십시오. 이들 함수에서 "ctf" 필터 유형 인수를 지정하고 B, A, g를 반환하도록 지정하여 스케일 값을 구합니다. (R2024b 이후)
참고 문헌
[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.
