bilinear

아날로그에서 디지털로의 필터 변환을 위한 쌍선형 변환 방법

구문

[zd,pd,kd] = bilinear(z,p,k,fs)

[numd,dend] = bilinear(num,den,fs)

[Ad,Bd,Cd,Dd] = bilinear(A,B,C,D,fs)

[___] = bilinear(___,fp)

설명

[zd,pd,kd] = bilinear(z,p,k,fs)는 z, p, k와 샘플 레이트 fs로 지정된 극점-영점 형식의 s 영역 전달 함수를 이에 상응하는 이산 함수로 변환합니다.

[numd,dend] = bilinear(num,den,fs)는 분자 num과 분모 den으로 지정된 s 영역 전달 함수를 이에 상응하는 이산 함수로 변환합니다.

[Ad,Bd,Cd,Dd] = bilinear(A,B,C,D,fs)는 A, B, C, D 행렬로 표현된 연속시간 상태공간 시스템을 이산시간 시스템으로 변환합니다.

[___] = bilinear(___,fp)는 파라미터 fp를 "일치" 주파수로 사용하여 사전 워핑을 지정합니다.

예제

체비쇼프 유형 I 아날로그 필터를 사용한 대역통과 IIR 필터 설계

6dB의 통과대역 리플을 갖는 10차 체비쇼프 유형 I 필터의 프로토타입을 설계합니다. 프로토타입을 상태공간 형식으로 변환합니다.

[z,p,k] = cheb1ap(10,6);
[A,B,C,D] = zp2ss(z,p,k);  

프로토타입을 대역통과 필터로 변환하여, ${\mathit{f}}_{\mathit{s}}=2\mathrm{kHz}$ 레이트로 샘플링될 때 상응하는 디지털 필터의 통과대역 경계가 100Hz와 500Hz가 되도록 합니다. 이 변환을 위해, 사전 워핑이 적용된 대역 경계 ${\mathit{u}}_1$과 ${\mathit{u}}_2$(단위: rad/s), 중심 주파수 ${\mathit{W}}_{\mathit{o}}=\sqrt{{\mathit{u}}_1{\mathit{u}}_2}$, 대역폭 ${\mathit{B}}_{\mathit{w}}={\mathit{u}}_2-{\mathit{u}}_1$을 지정합니다.

fs = 2e3;

f1 = 100; u1 = 2*fs*tan(f1*(2*pi/fs)/2);
f2 = 500; u2 = 2*fs*tan(f2*(2*pi/fs)/2);

[At,Bt,Ct,Dt] = lp2bp(A,B,C,D,sqrt(u1*u2),u2-u1);

freqs를 사용하여 아날로그 필터의 주파수 응답을 계산합니다. 크기 응답과 사전 워핑이 적용된 주파수 대역 경계를 플로팅합니다.

[b,a] = ss2tf(At,Bt,Ct,Dt);
[h,w] = freqs(b,a,2048);

plot(w,mag2db(abs(h)))
xline([u1 u2],"-",["Lower" "Upper"]+" passband edge", ...
    LabelVerticalAlignment="middle")

ylim([-165 5])
xlabel("Angular frequency (rad/s)")
ylabel("Magnitude (dB)")
grid

bilinear 함수를 사용하여 샘플 레이트가 ${\mathit{f}}_{\mathit{s}}$인 디지털 대역통과 필터를 만듭니다.

[Ad,Bd,Cd,Dd] = bilinear(At,Bt,Ct,Dt,fs);

freqz를 사용하여 디지털 필터를 상태공간 형식에서 2차섹션형(SOS)으로 변환하고 주파수 응답을 계산합니다. 크기 응답과 통과대역 경계를 플로팅합니다.

[hd,fd] = freqz(ss2sos(Ad,Bd,Cd,Dd),2048,fs);

plot(fd,mag2db(abs(hd)))
xline([f1 f2],"-",["Lower" "Upper"]+" passband edge", ...
    LabelVerticalAlignment="middle")

ylim([-165 5])
xlabel("Frequency (Hz)")
ylabel("Magnitude (dB)")
grid

타원 필터의 이산시간 표현

5dB의 통과대역 리플, 90dB의 저지대역 감쇠량 및 차단 주파수 ${\mathit{f}}_{\mathit{c}}=20\mathrm{Hz}$를 갖는 6차 타원 아날로그 저역통과 필터를 설계합니다.

fc = 20;

[z,p,k] = ellip(6,5,90,2*pi*fc,"s");

아날로그 타원 필터의 크기 응답을 시각화합니다. 차단 주파수를 표시합니다.

[num,den] = zp2tf(z,p,k);
[h,w] = freqs(num,den,1024);

plot(w/(2*pi),mag2db(abs(h)))
xline(fc,Color=[0.8500 0.3250 0.0980])

axis([0 100 -125 5])
grid
legend(["Magnitude response" "Cutoff frequency"])
xlabel("Frequency (Hz)")
ylabel("Magnitude (dB)")

bilinear 함수를 사용하여 아날로그 필터를 이산시간 IIR 필터로 변환합니다. 샘플 레이트 ${\mathit{f}}_{\mathit{s}}=200\mathrm{Hz}$와 사전 워핑 일치 주파수 ${\mathit{f}}_{\mathit{p}}=20\mathrm{Hz}$를 지정합니다.

fs = 200;
fp = 20;

[zd,pd,kd] = bilinear(z,p,k,fs,fp);

이산시간 필터의 크기 응답을 시각화합니다. 차단 주파수를 표시합니다.

[hd,fd] = freqz(zp2sos(zd,pd,kd),[],fs);

plot(fd,mag2db(abs(hd)))
xline(fc,Color=[0.8500 0.3250 0.0980])

axis([0 100 -125 5])
grid
legend(["Magnitude response" "Cutoff frequency"])
xlabel("Frequency (Hz)")
ylabel("Magnitude (dB)")

입력 인수

z, p, k — 아날로그 영역 영점, 극점, 이득

열 벡터, 스칼라

s 영역 전달 함수의 영점, 극점, 이득으로, 두 개의 열 벡터와 하나의 스칼라로 지정됩니다.

fs — 샘플 레이트

양의 스칼라

샘플 레이트로, 양의 스칼라로 지정됩니다.

num, den — 아날로그 영역 분자 및 분모 계수

행 벡터

아날로그 전달 함수의 분자 및 분모 계수로, 행 벡터로 지정됩니다.

A, B, C, D — 아날로그 영역 상태공간 행렬

행렬

s 영역의 상태공간 표현으로, 행렬로 지정됩니다. 시스템이 p개 입력값과 q개 출력값을 갖고 n개 상태 변수로 설명된다면 A는 n×n, B는 n×p, C는 q×n, D는 q×p입니다.

데이터형: single | double

fp — 일치 주파수

양의 스칼라

일치 주파수로, 양의 스칼라로 지정됩니다.

출력 인수

zd, pd, kd — 디지털 영역 영점, 극점, 이득

열 벡터, 스칼라

z 영역 전달 함수의 영점, 극점, 이득으로, 열 벡터와 스칼라로 반환됩니다.

numd, dend — 디지털 영역 분자 및 분모 계수

행 벡터

디지털 전달 함수의 분자 및 분모 계수로, 행 벡터로 반환됩니다.

Ad, Bd, Cd, Dd — 디지털 영역 상태공간 행렬

행렬

z 영역의 상태공간 표현으로, 행렬로 반환됩니다. 시스템이 n개 상태 변수로 설명되고 q개 출력값을 갖는 경우 Ad는 n×n, Bd는 n×1, Cd는 q×n, Dd는 q×1입니다.

데이터형: single | double

알고리즘

쌍선형 변환은 변수에 대한 수학적 매핑입니다. 디지털 필터링에서 이는 아날로그 평면 s를 디지털 평면 z에 매핑하는 표준 방법입니다. 이는 고전적인 필터 설계 기법을 사용하여 설계된 아날로그 필터를 이에 상응하는 이산 필터로 변환합니다.

쌍선형 변환은 다음과 같이 s 평면을 z 평면에 매핑합니다.

$${H(z)=H(s)|}_{s=2f_s\times \frac{z-1}{z+1}}.$$

이 변환은 다음과 같이 jΩ 축(범위: Ω = –∞ ~ +∞)을 단위원(e^jω, 범위: ω = –π ~ π) 주위에 반복해서 매핑합니다.

$$\omega =2{\tan }^{-1}\left(\frac{\Omega }{2f_s}\right).$$

bilinear는 사전 워핑을 지정하는 선택적 파라미터 f_p를 받을 수 있습니다. f_p(단위: 헤르츠)는 매핑 전의 주파수 응답과 매핑 후의 주파수 응답이 정확히 일치하는 일치 주파수를 나타냅니다. 사전 워핑이 적용된 모드에서 쌍선형 변환은 다음을 사용하여 s 평면을 z 평면에 매핑합니다.

$${H(z)=H(s)|}_{s=\frac{2\pi f_p}{\tan \left(\pi \frac{f_p}{f_s}\right)}\times \frac{z-1}{z+1}}.$$

사전 워핑 옵션을 사용하면 bilinear가 다음과 같이 jΩ 축(범위: Ω = –∞ ~ +∞)을 단위원(e^jω, 범위: ω = –π ~ π) 주위에 반복해서 매핑합니다.

$$\omega =2{\tan }^{-1}\left(\frac{\Omega \tan \left(\pi \frac{f_p}{f_s}\right)}{2\pi f_p}\right).$$

사전 워핑이 적용된 모드에서 bilinear는 s 평면의 주파수 2πf_p(단위: 초당 라디안)를 z 평면의 정규화 주파수 2πf_p/f_s(단위 초당 라디안)와 일치시킵니다.

bilinear 함수는 영점-극점-이득 형식, 전달 함수 형식, 상태공간 형식의 세 가지 선형 시스템 표현에 대해 사용할 수 있습니다.

bilinear는 사용자가 제공하는 입력 선형 시스템의 형식에 따라 두 알고리즘 중 하나를 사용합니다. 한 알고리즘은 영점-극점-이득 형식에서 작동하고 다른 알고리즘은 상태공간 형식에서 작동합니다. 전달 함수 표현인 경우에는 bilinear는 이를 상태공간 형식으로 전환하여 변환을 수행한 후, 결과로 생성된 상태공간 시스템을 다시 전달 함수 형식으로 전환합니다.

영점-극점-이득 알고리즘

영점-극점-이득 형식의 시스템에 대해 bilinear는 다음 4단계를 수행합니다.

1. fp가 존재하는 경우 다음과 같이 사전 워핑을 수행합니다.

   fp = 2*pi*fp;
   fs = fp/tan(fp/fs/2)
   

   그렇지 않으면 fs = 2*fs입니다.

2. 다음을 사용하여 ±∞에서 모든 영점을 제거합니다.

   z = z(finite(z));
   

3. 다음을 사용하여 영점, 극점, 이득을 변환합니다.

   pd = (1+p/fs)./(1-p/fs);    % Do bilinear transformation
   zd = (1+z/fs)./(1-z/fs);
   kd = real(k*prod(fs-z)./prod(fs-p));
   

4. -1에 영점을 추가하여 결과로 생성되는 시스템에서 분자와 분모의 차수가 동일해지도록 합니다.

상태공간 알고리즘

아날로그 시스템에 대한 다음 상태공간 방정식은

$$\begin{array}{l}\dot{x}=Ax+Bu\\ y=Cx+Du\end{array}$$

다음과 같은 이산 형식으로 변환할 수 있습니다.

$$\begin{array}{l}x[n+1]=A_{d}x[n]+B_{d}u[n],\\ y[n]=C_{d}x[n]+D_{d}u[n].\end{array}$$

변환을 위해 bilinear는 다음 두 단계를 수행합니다[1].

1. 일치 주파수 f_p가 지정된 경우 다음과 같이 설정합니다.

   $$\lambda =\frac{\pi f_p}{\tan (\pi f_p/f_s)}.$$

   f_p가 지정되지 않은 경우 다음과 같이 설정합니다. λ = f_s

2. 다음을 사용하여 A, B, C, D로 A_d, B_d, C_d, D_d를 계산합니다.

   $$\begin{array}{l}{A}_d={(}^I(I+A\frac{1}{2\lambda }),\\ B_d=\frac{1}{\sqrt{\lambda }}{(}^{I}B,\\ C_d=\frac{1}{\sqrt{\lambda }}C{(}^I,\\ D_d=\frac{1}{2\lambda }C{(}^{I}B+D.\end{array}$$

전달 함수

전달 함수 형식의 시스템에 대해 bilinear는 num과 den으로 주어진 s 영역 전달 함수를 이에 상응하는 이산 함수로 변환합니다. 행 벡터 num과 den은 각각 s의 거듭제곱 내림차순으로 분자와 분모의 계수를 지정합니다. B(s)를 분자 다항식으로 지정하고 A(s)를 분모 다항식으로 지정합니다. 전달 함수는 다음과 같습니다.

$$\frac{B(s)}{A(s)}=\frac{B(1)s^n+\cdots +B(n)s+B(n+1)}{A(1)s^m+\cdots +A(m)s+A(m+1)}$$

fs는 헤르츠 단위의 샘플 레이트입니다. bilinear는 상응하는 이산 함수를 z의 거듭제곱 내림차순(z^–1의 거듭제곱 오름차순)으로 정렬된 행 벡터 numd와 dend로 반환합니다. fp는 선택 사항이며 사전 워핑에 사용되는 헤르츠 단위의 일치 주파수입니다.

참고 문헌

[1] Al-Saggaf, Ubaid M., and Gene F. Franklin. “Model Reduction via Balanced Realizations: An Extension and Frequency Weighting Techniques.” IEEE Transactions on Automatic Control 33, no. 7 (July 1988): 687–92. https://doi.org/10.1109/9.1280.

[2] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice Hall, 1999.

[3] Parks, Thomas W., and C. Sidney Burrus. Digital Filter Design. New York: John Wiley & Sons, 1987.

[4] Tustin, Arnold. “A Method of Analysing the Behaviour of Linear Systems in Terms of Time Series.” Journal of the Institution of Electrical Engineers - Part IIA: Automatic Regulators and Servo Mechanisms 94, no. 1 (May 1947): 130–42. https://doi.org/10.1049/ji-2a.1947.0020.

확장 기능

C/C++ 코드 생성
MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

버전 내역

R2006a 이전에 개발됨

참고 항목

impinvar | lp2bp | lp2bs | lp2hp | lp2lp

도움말 항목

• Continuous-Discrete Conversion Methods (Control System Toolbox)