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

tfestimate

전달 함수 추정값

설명

txy = tfestimate(x,y)는 입력 신호 x와 출력 신호 y가 주어진 경우 전달 함수 추정값 txy를 구합니다.

  • xy가 모두 벡터인 경우 길이가 같아야 합니다.

  • 두 신호 중 하나는 행렬이고 다른 하나는 벡터인 경우, 벡터의 길이는 행렬에 포함된 행 개수와 동일해야 합니다. 이 함수는 벡터를 확장하고 열마다 단위 전달 함수 추정값을 가진 행렬을 반환합니다.

  • xy가 행의 개수가 같지만 열의 개수가 다른 행렬인 경우 txy는 모든 입력 신호와 출력 신호를 결합하는 MIMO(다중 입력/다중 출력) 전달 함수입니다. txy는 3차원 배열입니다. x는 m개의 열을 가지고 y는 n개의 열을 가지는 경우, txy는 n개의 열과 m개의 페이지를 가집니다. 자세한 내용은 전달 함수 항목을 참조하십시오.

  • xy가 같은 크기의 행렬인 경우, tfestimate는 각 열에 대해 연산을 수행합니다. 즉, txy(:,n) = tfestimate(x(:,n),y(:,n))입니다. MIMO 추정값을 구하려면 인수 목록에 'mimo'를 추가하십시오.

예제

txy = tfestimate(x,y,window)window를 사용하여 xy를 세그먼트로 나누고 윈도우를 적용합니다.

txy = tfestimate(x,y,window,noverlap)은 인접 세그먼트 간에 noverlap개의 샘플이 중첩되도록 적용합니다.

txy = tfestimate(x,y,window,noverlap,nfft)nfft개 샘플링 점을 사용하여 이산 푸리에 변환을 계산합니다.

txy = tfestimate(___,'mimo')는 행렬 입력값에 대해 MIMO 전달 함수를 계산합니다. 이 구문은 위에 열거된 구문에 사용할 수 있습니다.

[txy,w] = tfestimate(___)는 전달 함수가 추정된 정규화 주파수로 구성된 벡터 w를 반환합니다.

예제

[txy,f] = tfestimate(___,fs)는 전달 함수가 추정된 주파수(샘플 레이트 fs의 값으로 표현됨)로 구성된 벡터 f를 반환합니다. fstfestimate에 대한 6번째 숫자형 입력값이어야 합니다. 샘플 레이트를 입력하고 위에 열거된 옵션 인수의 디폴트 값을 그대로 사용하려면 이 인수를 빈 값 []로 지정하십시오.

[txy,w] = tfestimate(x,y,window,noverlap,w)w에 지정된 정규화 주파수에서 전달 함수 추정값을 반환합니다.

[txy,f] = tfestimate(x,y,window,noverlap,f,fs)f에 지정된 주파수에서 전달 함수 추정값을 반환합니다.

[___] = tfestimate(x,y,___,freqrange)freqrange로 지정된 주파수 범위에 대한 전달 함수 추정값을 반환합니다. freqrange에 유효한 옵션은 'onesided', 'twosided', 'centered'입니다.

예제

[___] = tfestimate(___,'Estimator',est)는 추정량 est를 사용하여 전달 함수를 추정합니다. est에 유효한 옵션은 'H1''H2'입니다.

tfestimate(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 전달 함수 추정값을 플로팅합니다.

예제

모두 축소

두 시퀀스 xy 간의 전달 함수 추정값을 계산하고 플로팅합니다. 시퀀스 x는 백색 가우스 잡음으로 구성되어 있습니다. y0.2π rad/sample의 정규화된 차단 주파수를 갖는 30차 저역통과 필터를 사용하여 x를 필터링한 결과입니다. 사각 윈도우를 사용하여 필터를 설계합니다. 전달 함수 추정값에 대해 샘플 레이트가 500Hz이고 길이가 1024인 해밍 윈도우를 지정합니다.

h = fir1(30,0.2,rectwin(31));
x = randn(16384,1);
y = filter(h,1,x);

fs = 500;
tfestimate(x,y,1024,[],[],fs)

fvtool을 사용하여 전달 함수가 필터의 주파수 응답에 근사하는지 확인합니다.

fvtool(h,1,'Fs',fs)

변수로 전달 함수 추정값을 반환하고 추정값의 절댓값을 데시벨 단위로 플로팅하여 동일한 결과를 구합니다.

[Txy,f] = tfestimate(x,y,1024,[],[],fs);

plot(f,mag2db(abs(Txy)))

간단한 SISO(단일 입력 단일 출력) 시스템에 대한 전달 함수를 추정하고 이를 정의와 비교합니다.

1차원의 이산시간 진동 시스템은 단위 탄성 상수를 갖는 용수철로 벽에 연결된 단위 질량 m으로 구성됩니다. 센서가 질량의 가속도 aFs=1Hz로 샘플링합니다. 감쇠기는 감쇠 상수 b=0.01을 사용하여 속도에 비례하는 힘을 가하여 질량의 운동을 저해합니다.

2000개의 시간 샘플을 생성합니다. 샘플링 간격을 Δt=1/Fs로 정의합니다.

Fs = 1;
dt = 1/Fs;
N = 2000;
t = dt*(0:N-1);
b = 0.01;

시스템은 다음 상태공간 모델로 설명할 수 있습니다.

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

여기서 x=[rv]T는 상태 벡터로, rv는 질량의 위치와 속도를 각각 나타내고, u는 구동력이며, y=a는 측정된 출력값입니다. 상태공간 행렬은 다음과 같습니다.

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[-1-b],D=1,

I2×2 단위 행렬이고, 연속시간 상태공간 행렬은 다음과 같습니다.

Ac=[01-1-b],Bc=[01].

Ac = [0 1;-1 -b];
A = expm(Ac*dt);

Bc = [0;1];
B = Ac\(A-eye(size(A)))*Bc;

C = [-1 -b];
D = 1;

질량은 측정 구간의 절반에 대한 임의의 입력값으로 구동됩니다. 상태공간 모델을 사용하여 모든 초기 상태값이 0인 시공간 시스템의 변화를 계산합니다. 질량의 가속도를 시간 함수로 플로팅합니다.

rng default

u = zeros(1,N);
u(1:N/2) = randn(1,N/2);

y = 0;
x = [0;0];
for k = 1:N
    y(k) = C*x + D*u(k);
    x = A*x + B*u(k);
end

plot(t,y)

시스템의 전달 함수를 주파수 함수로 추정합니다. 2048개의 DFT 점을 사용하고 형태 인자가 15인 카이저 윈도우를 지정합니다. 인접 세그먼트 간에 디폴트 중첩값을 적용합니다.

nfs = 2048;
wind = kaiser(N,15);

[txy,ft] = tfestimate(u,y,wind,[],nfs,Fs);

이산시간 시스템의 주파수-응답 함수는 단위원에서 계산된, 시스템의 시간 영역 전달 함수에 대한 z 변환으로 표현할 수 있습니다. tfestimate로 계산된 추정값이 이 정의와 일치하는지 확인합니다.

[b,a] = ss2tf(A,B,C,D);

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);
frf = polyval(b,z)./polyval(a,z);

plot(ft,20*log10(abs(txy)))
hold on
plot(fz,20*log10(abs(frf)))
hold off
grid
ylim([-60 40])

tfestimate의 내장 기능을 사용하여 추정값을 플로팅합니다.

tfestimate(u,y,wind,[],nfs,Fs)

간단한 다중 입력/다중 출력 시스템에 대해 전달 함수를 추정합니다.

이상적인 1차원 진동 시스템은 두 벽 사이에 있는 두 개의 질량 m1m2로 구성됩니다. 각각은 m1=1이고 m2=μ입니다. 각 질량은 탄성 상수 k를 갖는 용수철로 가장 가까운 벽에 연결됩니다. 동일한 용수철이 두 질량을 연결합니다. 세 개의 감쇠기가 감쇠 상수 b를 사용하여 속도에 비례하는 힘을 가하여 질량의 운동을 저해합니다. 센서가 질량의 가속도인 a1a2Fs=50Hz에서 샘플링합니다.

600초에 해당하는 30000개 시간 샘플을 생성합니다. 샘플링 간격을 Δt=1/Fs로 정의합니다.

Fs = 50;
dt = 1/Fs;
N = 30000;
t = dt*(0:N-1);

시스템은 다음 상태공간 모델로 설명할 수 있습니다.

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

여기서 x=[r1v1r2v2]T는 상태 벡터로, rivi는 각각 i번째 질량의 위치와 속도를 나타내고, u=[u1u2]T는 입력 구동력 벡터이며, y=[a1a2]T는 출력 벡터입니다. 상태공간 행렬은 다음과 같습니다.

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[-2k-2bkbk/μb/μ-2k/μ-2b/μ],D=[1001/μ],

I4×4 단위 행렬이고, 연속시간 상태공간 행렬은 다음과 같습니다.

Ac=[0100-2k-2bkb0001k/μb/μ-2k/μ-2b/μ],Bc=[00100001/μ].

k=400, b=0, μ=1/10로 설정합니다.

k = 400;
b = 0;
m = 1/10;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/m];
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];
D = [1 0;0 1/m];

질량은 전체 측정에서 임의의 입력값으로 구동됩니다. 상태공간 모델을 사용하여 모든 초기 상태값이 0인 시공간 시스템의 변화를 계산합니다.

rng default
u = randn(2,N);

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

입력 데이터와 출력 데이터를 사용하여 시스템의 전달 함수를 주파수 함수로 추정합니다. 'mimo' 옵션을 지정하여 네 개의 모든 전달 함수를 생성합니다. 5000개 샘플로 구성된 핸 윈도우를 사용하여 신호를 여러 세그먼트로 나눕니다. 인접 세그먼트 간에 2500개 샘플이 중첩되도록 지정하고 214개 DFT 점을 지정합니다. 추정값을 플로팅합니다.

wind = hann(5000);
nov = 2500;

[q,fq] = tfestimate(u',y',wind,nov,2^14,Fs,'mimo');

단위원에서 계산된, 시간 영역 전달 함수의 z 변환으로 이론상의 전달 함수를 계산합니다.

nfs = 2^14;

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

이론상의 전달 함수와 이에 대응되는 추정값을 플로팅합니다.

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(fq,20*log10(abs(q(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        grid
        title(['Input ' int2str(kj) ', Output ' int2str(jk)])
        axis([0 Fs/2 -50 100])
    end
end

전달 함수는 예상 값 ω1,2/2π에서 최댓값을 가집니다. 여기서 ω는 모드 행렬의 고유값입니다.

sqrt(eig(k*[2 -1;-1/m 2/m]))/(2*pi)
ans = 2×1

    3.8470
   14.4259

b=0.1로 설정하여 시스템에 감쇠를 추가합니다. 동일한 구동력으로 감쇠된 시스템의 시간 변화를 계산합니다. 동일한 윈도우와 중첩을 사용하여 MIMO 전달 함수의 H2 추정값을 계산합니다. tfestimate 기능을 사용하여 추정값을 플로팅합니다.

b = 0.1;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

clf
tfestimate(u',y',wind,nov,[],Fs,'mimo','Estimator','H2')
legend('I1, O1','I1, O2','I2, O1','I2, O2')

yl = ylim;

추정값을 이론상의 예측값과 비교합니다.

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

plot(fz*Fs,20*log10(abs(reshape(permute(frf,[2 1 3]),[nfs/2 4]))))
legend('I1, O1','I1, O2','I2, O1','I2, O2')
ylim(yl)
grid

입력 인수

모두 축소

입력 신호로, 벡터나 행렬로 지정됩니다.

예: cos(pi/4*(0:159))+randn(1,160)은 백색 가우스 잡음을 포함하는 정현파를 지정합니다.

데이터형: single | double
복소수 지원 여부:

출력 신호로, 벡터나 행렬로 지정됩니다.

데이터형: single | double
복소수 지원 여부:

윈도우로, 정수로 지정되거나 행 벡터 또는 열 벡터로 지정됩니다. window를 사용하여 신호를 여러 세그먼트로 나눕니다.

  • window가 정수이면 tfestimatexy를 길이가 window인 세그먼트로 나누고 각 세그먼트에 해당 길이의 해밍 윈도우를 적용합니다.

  • window가 벡터이면 tfestimatexy를 벡터와 길이가 같은 세그먼트로 나누고 window를 사용하여 각 세그먼트에 윈도우를 적용합니다.

xy의 길이를 noverlap개의 샘플이 중첩된 정수 개수의 세그먼트로 정확히 나눌 수 없는 경우 신호가 이에 맞게 적절하게 잘립니다.

window를 빈 값으로 지정하면, tfestimate는 해밍 윈도우를 사용하여 xynoverlap개의 중첩 샘플을 가지는 8개의 세그먼트로 나눕니다.

사용 가능한 윈도우 목록은 윈도우 항목을 참조하십시오.

예: hann(N+1)(1-cos(2*pi*(0:N)'/N))/2는 모두 길이가 N + 1인 핸 윈도우를 지정합니다.

데이터형: single | double

중첩(Overlap) 샘플의 개수로, 양의 정수로 지정됩니다.

  • window가 스칼라이면 noverlapwindow보다 작아야 합니다.

  • window가 벡터이면 noverlapwindow의 길이보다 작아야 합니다.

noverlap을 빈 값으로 지정하면 tfestimate는 세그먼트 간에 50% 중첩을 만드는 값을 사용합니다. 세그먼트 길이가 지정되지 않은 경우 이 함수는 noverlap을 ⌊N/4.5⌋로 설정합니다. 여기서 N은 입력 신호와 출력 신호의 길이입니다.

데이터형: double | single

DFT를 적용할 지점의 개수로, 양의 정수로 지정됩니다. nfft를 빈 값으로 지정하면 tfestimate는 이 인수를 max(256,2p)로 설정합니다. 여기서, 길이가 N인 입력 신호의 경우 p = ⌈log2 N⌉입니다.

데이터형: single | double

샘플 레이트로, 양의 스칼라로 지정됩니다. 샘플 레이트는 단위 시간당 샘플 개수입니다. 시간 단위가 초이면 샘플 레이트의 단위는 Hz입니다.

정규화 주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 정규화 주파수는 rad/sample을 단위로 사용합니다.

예: w = [pi/4 pi/2]

데이터형: double

주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 주파수는 단위 시간당 주기를 단위로 사용합니다. 단위 시간은 샘플 레이트 fs로 지정됩니다. fs의 단위가 샘플/초이면 f는 Hz 단위입니다.

예: fs = 1000; f = [100 200]

데이터형: double

전달 함수 추정값의 주파수 범위로, 'onesided', 'twosided' 또는 'centered'로 지정됩니다. 디폴트 값은 실수 값을 갖는 신호의 경우 'onesided'이고, 복소수 값을 갖는 신호의 경우 'twosided'입니다.

  • 'onesided' — 실수 값을 갖는 두 개의 입력 신호 xy 간 전달 함수의 단측 추정값을 반환합니다. nfft가 짝수이면 txynfft/2 + 1개의 행을 가지고 구간 [0,π] rad/sample에 대해 계산됩니다. nfft가 홀수이면 txy가 (nfft + 1)/2개의 행을 가지고 구간은 [0,π) rad/sample이 됩니다. fs를 지정하면 이에 대응되는 구간은 nfft가 짝수인 경우 [0,fs/2] cycles/unit time이고 nfft가 홀수인 경우 [0,fs/2) cycles/unit time입니다.

  • 'twosided' — 실수 값 또는 복소수 값을 갖는 두 개의 입력 신호 xy 간 전달 함수의 양측 추정값을 반환합니다. 이 경우, txynfft개의 행을 가지고 구간 [0,2π) rad/sample에 대해 계산됩니다. fs를 지정하면 구간은 [0,fs) cycles/unit time이 됩니다.

  • 'centered' — 실수 값 또는 복소수 값을 갖는 두 개의 입력 신호 xy 간 전달 함수의, 중심이 맞춰진 양측 추정값을 반환합니다. 이 경우, txynfft개의 행을 가지고, nfft가 짝수이면 구간 (–π,π] rad/sample에 대해 계산되며, nfft가 홀수이면 구간 (–π,π) rad/sample에 대해 계산됩니다. fs를 지정하면 이에 대응되는 구간은 nfft가 짝수인 경우 (–fs/2, fs/2] cycles/unit time이고 nfft가 홀수인 경우 (–fs/2, fs/2) cycles/unit time입니다.

전달 함수 추정량으로, 'H1' 또는 'H2'로 지정됩니다.

  • 잡음이 입력 신호와 상관 관계가 없는 경우 'H1'을 사용합니다.

  • 잡음이 출력 신호와 상관 관계가 없는 경우 'H2'를 사용합니다. 이 경우, 입력 신호의 개수는 출력 신호의 개수와 같아야 합니다.

자세한 내용은 전달 함수 항목을 참조하십시오.

출력 인수

모두 축소

전달 함수 추정값으로, 벡터, 행렬 또는 3차원 배열로 반환됩니다.

정규화 주파수로, 실수 열 벡터로 반환됩니다.

주기적 주파수로, 실수 열 벡터로 반환됩니다.

세부 정보

모두 축소

전달 함수

입력값 x와 출력값 y 간의 관계는 선형 시불변 전달 함수 txy로 모델링됩니다. 주파수 영역에서 Y(f) = H(f)X(f)입니다.

  • SISO(단일 입력 단일 출력) 시스템의 경우, 전달 함수의 H1 추정값은 다음과 같이 지정됩니다.

    H1(f)=Pyx(f)Pxx(f),

    여기서 Pyxxy의 상호 전력 스펙트럼 밀도이고, Pxxx의 전력 스펙트럼 밀도입니다. 이 추정값은 잡음이 시스템 입력값과 상관관계가 없다고 가정합니다.

    MIMO(다중 입력/다중 출력) 시스템의 경우, H1 추정량은 m개 입력값과 n개 출력값에 대해 다음과 같이 됩니다.

    H1(f)=PYX(f)PXX1(f)=[Py1x1(f)Py1x2(f)Py1xm(f)Py2x1(f)Py2x2(f)Py2xm(f)Pynx1(f)Pynx2(f)Pynxm(f)][Px1x1(f)Px1x2(f)Px1xm(f)Px2x1(f)Px2x2(f)Px2xm(f)Pxmx1(f)Pxmx2(f)Pxmxm(f)]1

    여기서,

    • Pyixk는 k번째 입력값과 i번째 출력값의 상호 전력 스펙트럼 밀도입니다.

    • Pxixk는 k번째 입력값과 i번째 입력값의 상호 전력 스펙트럼 밀도입니다.

    입력값과 출력값이 각각 두 개인 경우, 추정량은 다음과 같은 행렬입니다.

    H1(f)=[Py1x1(f)Px2x2(f)Py1x2(f)Px2x1(f)Py1x2(f)Px1x1(f)Py1x1(f)Px1x2(f)Py2x1(f)Px2x2(f)Py2x2(f)Px2x1(f)Py2x2(f)Px1x1(f)Py2x1(f)Px1x2(f)]Px1x1(f)Px2x2(f)Px1x2(f)Px2x1(f).

  • SISO(단일 입력 단일 출력) 시스템의 경우, 전달 함수의 H2 추정값은 다음과 같이 지정됩니다.

    H2(f)=Pyy(f)Pxy(f),

    여기서 Pyyy의 전력 스펙트럼 밀도이고, Pxy = P*yxxy에 대한 상호 전력 스펙트럼 밀도의 켤레 복소수입니다. 이 추정값은 잡음이 시스템 출력값과 상관관계가 없다고 가정합니다.

    MIMO 시스템의 경우, H2 추정량은 입력값과 출력값의 개수가 같은 경우, 즉 n = m인 경우에만 제대로 정의됩니다. 이 추정량은 다음과 같이 됩니다.

    H2(f)=PYY(f)PXY1(f)=[Py1y1(f)Py1y2(f)Py1yn(f)Py2y1(f)Py2y2(f)Py2yn(f)Pyny1(f)Pyny2(f)Pynyn(f)][Px1y1(f)Px1y2(f)Px1yn(f)Px2y1(f)Px2y2(f)Px2yn(f)Pxny1(f)Pxny2(f)Pxnyn(f)]1,

    여기서,

    • Pyiyk는 k번째 출력값과 i번째 출력값의 상호 전력 스펙트럼 밀도입니다.

    • Pxiyk는 i번째 입력값과 k번째 출력값에 대한 상호 전력 스펙트럼 밀도의 켤레 복소수입니다.

알고리즘

tfestimate는 웰치의 평균화된 주기도 방법을 사용합니다. 자세한 내용은 pwelch 항목을 참조하십시오.

참고 문헌

[1] Vold, Håvard, John Crowley, and G. Thomas Rocklin. “New Ways of Estimating Frequency Response Functions.” Sound and Vibration. Vol. 18, November 1984, pp. 34–38.

참고 항목

| | |

R2006a 이전에 개발됨