이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.

modalfit

주파수-응답 함수의 모드 파라미터

구문

fn = modalfit(frf,f,fs,mnum)

fn = modalfit(frf,f,fs,mnum,Name,Value)

[fn,dr,ms] = modalfit(___)

[fn,dr,ms,ofrf] = modalfit(___)

[___] = modalfit(sys,f,mnum,Name,Value)

설명

fn = modalfit(frf,f,fs,mnum)은 시스템의 mnum개 모드에 대한 고유 주파수를 추정합니다. 시스템은 샘플 레이트 fs 및 주파수 f에서 정의된 주파수-응답 함수 frf를 측정한 것입니다. modalfrf를 사용하여 측정 데이터에서 주파수-응답 함수의 행렬을 생성합니다. frf는 동적 유연성(리셉턴스) 형식이라고 가정합니다.

fn = modalfit(frf,f,fs,mnum,Name,Value)는 이름-값 인수를 사용하여 추가 옵션을 지정합니다.

[fn,dr,ms] = modalfit(___)은 위 구문의 입력값 조합을 사용하여 fn의 각 고유 주파수에 해당하는 모드 형상 벡터와 감쇠비도 반환합니다.

[fn,dr,ms,ofrf] = modalfit(___)은 추정된 모드 파라미터를 기반으로 재구성된 주파수-응답 함수 배열도 반환합니다.

[___] = modalfit(sys,f,mnum,Name,Value)는 식별된 모델 sys의 모드 파라미터를 추정합니다. ssest (System Identification Toolbox) 또는 tfest (System Identification Toolbox) 같은 추정 명령을 사용하여, 측정된 주파수-응답 함수 또는 시간 영역의 입력 신호와 출력 신호에서 sys를 생성합니다. 이 구문에서는 'DriveIndex', 'FreqRange', 'PhysFreq' 이름-값 인수를 사용할 수 있습니다. 일반적으로 이 방식은 비모수적 방법을 사용하는 구문보다 데이터를 더 적게 필요로 합니다. 이 구문을 사용하려면 System Identification Toolbox™ 라이선스가 있어야 합니다.

예제

SISO 시스템의 주파수-응답 함수

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

1차원의 이산시간 진동 시스템은 탄성 상수 $\mathit{k}=1$을 갖는 스프링으로 벽에 연결된 단위 질량 $\mathit{m}$으로 구성됩니다. 센서가 질량의 변위를 ${\mathit{F}}_{\mathrm{s}}=1$Hz로 샘플링합니다. 댐퍼는 감쇠 상수 $\mathit{b}=0.01$을 사용하여 속도에 비례하는 힘을 가해 질량의 운동을 저해합니다.

3000개의 시간 샘플을 생성합니다. 샘플링 간격을 $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$로 정의합니다.

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

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

$$\begin{array}{c}x(k+1)=Ax(k)+Bu(k),\\ y(k)=Cx(k)+Du(k),\end{array}$$

여기서 $\mathit{x}={\left[\begin{array}{cc}\mathit{r}& \mathit{v}\end{array}\right]}^{\mathit{T}}$는 상태 벡터로, $\mathit{r}$과 $\mathit{v}$는 질량의 변위와 속도를 각각 나타내고, $\mathit{u}$는 구동력이며, $\mathit{y}=\mathit{r}$은 측정된 출력값입니다. 상태공간 행렬은 다음과 같습니다.

$$A=\exp (A_c\Delta t),B=A_c^{-1}(A-I)B_c,C=\left[\begin{array}{cc}1& 0\end{array}\right],D=0,$$

$\mathit{I}$는 $2\times 2$ 단위 행렬이고, 연속시간 상태공간 행렬은 다음과 같습니다.

$$A_c=\left[\begin{array}{cc}0& 1\\ -1& -b\end{array}\right],B_c=\left[\begin{array}{c}0\\ 1\end{array}\right].$$

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

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

C = [1 0];
D = 0;

처음 2000초 동안은 임의의 입력값으로 질량을 구동시키고 그다음에는 자유롭게 움직이도록 둡니다. 상태공간 모델을 사용하여 모든 초기 상태값이 0인 시공간 시스템의 변화를 계산합니다. 질량의 변위를 시간 함수로 플로팅합니다.

rng default
u = randn(1,N)/2;
u(2001:end) = 0;

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)

시스템의 모드 주파수-응답 함수를 추정합니다. 측정된 신호 길이의 1/2인 핸 윈도우를 사용합니다. 출력값이 질량의 변위임을 지정합니다.

wind = hann(N/2);

[frf,f] = modalfrf(u',y',Fs,wind,'Sensor','dis');

이산시간 시스템의 주파수-응답 함수는 단위원에서 계산된, 시스템의 시간 영역 전달 함수에 대한 Z 변환으로 표현할 수 있습니다. modalfrf 추정값을 정의와 비교합니다.

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

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

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

진동 모드의 고유 주파수와 감쇠비를 추정합니다.

[fn,dr] = modalfit(frf,f,Fs,1,'FitMethod','PP')

fn = 0.1593

dr = 0.0043

고유 주파수를 비감쇠 시스템의 이론상 값인 $1/2\pi$와 비교합니다.

theo = 1/(2*pi)

theo = 0.1592

최소제곱(Least-squares) 유리 함수 방법을 사용한 모드 파라미터

우주 정거장 모듈의 주파수-응답 함수(FRF) 배열에서 시작하여 모드 파라미터를 계산합니다.

3개 입력값/3개 출력값으로 구성된 FRF 배열을 포함하는 구조체를 불러옵니다. 시스템은 320Hz로 샘플링되었습니다.

load modaldata SpaceStationFRF

frf = SpaceStationFRF.FRF; 
f = SpaceStationFRF.f; 
fs = SpaceStationFRF.Fs;

최소제곱 유리 함수 방법을 사용하여 가장 낮은 24개 모드의 모드 파라미터를 추출합니다.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,24,'FitMethod','lsrf');

복원한 FRF 배열을 측정된 FRF 배열과 비교합니다.

for ij = 1:3
 for ji = 1:3
    subplot(3,3,3*(ij-1)+ji)
    loglog(f,abs(frf(:,ji,ij)))
    hold on
    loglog(f,abs(ofrf(:,ji,ij)))
    hold off
    axis tight
    title(sprintf('In%d -> Out%d',ij,ji))
    if ij==3
        xlabel('Frequency (Hz)')
    end
 end
end

2체(Two-Body) 진동자의 모드 파라미터

간단한 MIMO(다중 입력/다중 출력) 시스템에 대해 주파수-응답 함수와 모드 파라미터를 추정합니다.

이상적인 1차원 진동 시스템은 두 벽 사이에 있는 두 개의 질량 ${\mathit{m}}_1$과 ${\mathit{m}}_2$로 구성됩니다. 각각은 ${\mathit{m}}_1=1$이고 ${\mathit{m}}_2=\mu$입니다. 각 질량은 탄성 상수 $\mathit{k}$를 갖는 스프링으로 가장 가까운 벽에 연결됩니다. 동일한 스프링이 두 질량을 연결합니다. 세 개의 댐퍼가 감쇠 상수 $\mathit{b}$를 사용하여 속도에 비례하는 힘을 가해 질량의 운동을 저해합니다. 센서가 질량의 변위인 ${\mathit{r}}_1$과 ${\mathit{r}}_2$를 ${\mathit{F}}_{\mathrm{s}}=50$Hz로 샘플링합니다.

600초에 해당하는 30,000개 시간 샘플을 생성합니다. 샘플링 간격을 $\Delta \mathit{t}=1/{\mathit{F}}_{\mathrm{s}}$로 정의합니다.

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

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

$$\begin{array}{c}x(k+1)=Ax(k)+Bu(k),\\ y(k)=Cx(k)+Du(k),\end{array}$$

$$x={\left[\begin{array}{cccc}{r}_1& v_1& r_2& v_2\end{array}\right]}^T$$는 상태 벡터, $$r_i$$ 및 $$v_i$$는 각각 $$i$$번째 질량의 위치와 속도, $$u={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$$는 입력 구동력의 벡터, $$y={\left[\begin{array}{cc}{r}_1& r_2\end{array}\right]}^T$$는 출력 벡터입니다. 상태공간 행렬은 다음과 같습니다.

$$A=\exp (A_c\Delta t),B=A_c^{-1}(A-I)B_c,C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],D=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],$$

$\mathit{I}$는 $4\times 4$ 단위 행렬이고, 연속시간 상태공간 행렬은 다음과 같습니다.

$$A_c=\left[\begin{array}{cccc}0& 1& 0& 0\\ -2k& -2b& k& b\\ 0& 0& 0& 1\\ k/\mu & b/\mu & -2k/\mu & -2b/\mu \end{array}\right],B_c=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 0\\ 0& 1/\mu \end{array}\right].$$

$\mathit{k}=400$, $\mathit{b}=0.1$, $\mu =1/10$로 설정합니다.

k = 400;
b = 0.1;
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 = [1 0 0 0;0 0 1 0];
D = zeros(2);

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

rng default
u = randn(2,N);

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

입력 데이터와 출력 데이터를 사용하여 시스템의 전달 함수를 주파수 함수로 추정합니다. 인접 세그먼트 간에 9000개 샘플이 중첩되는, 15000개 샘플로 구성된 핸 윈도우를 사용합니다. 측정된 출력값이 변위임을 지정합니다.

wind = hann(15000);
nove = 9000;
[FRF,f] = modalfrf(u',y',Fs,wind,nove,'Sensor','dis');

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

nfs = 2048;
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(f,20*log10(abs(FRF(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        axis([0 Fs/2 -100 0])
        title(sprintf('Input %d, Output %d',jk,kj))
    end
end

출력 인수 없이 modalfrf 구문을 사용하여 추정값을 플로팅합니다.

figure
modalfrf(u',y',Fs,wind,nove,'Sensor','dis')

시스템의 고유 주파수, 감쇠비, 모드 형상을 추정합니다. 피크 추출(peak-picking) 방법을 사용하여 계산합니다.

[fn,dr,ms] = modalfit(FRF,f,Fs,2,'FitMethod','pp');
fn

fn = 
fn(:,:,1) =

    3.8466    3.8466
    3.8495    3.8495


fn(:,:,2) =

    3.8492    3.8490
    3.8552   14.4684

고유 주파수를 비감쇠 시스템의 이론상 예측값과 비교합니다.

undamped = sqrt(eig([2*k -k;-k/m 2*k/m]))/2/pi

undamped = 2×1

    3.8470
   14.4259

MIMO 시스템의 모드 파라미터

몇 번의 랜덤 잡음 버스트가 가해진, 2입력/3출력 시스템에 대해 고유 주파수, 감쇠비, 모드 형상을 계산합니다. 각 버스트는 1초간 지속되고, 각 버스트가 끝나고 2초가 지난 다음 버스트가 시작되었습니다. 데이터는 4kHz로 샘플링되었습니다.

데이터 파일을 불러옵니다. 입력 신호와 출력 신호를 플로팅합니다.

load modaldata

subplot(2,1,1)
plot(Xburst)
title('Input Signals')
subplot(2,1,2)
plot(Yburst)
title('Output Signals')

주파수-응답 함수를 계산합니다. 길이가 버스트 기간과 동일하고 인접 세그먼트 간에 중첩이 없는 사각 윈도우를 지정합니다.

burstLen = 12000;
[frf,f] = modalfrf(Xburst,Yburst,fs,burstLen);

안정화 다이어그램을 시각화하고, 안정된 고유 주파수를 반환합니다. 최대 모델 차수를 30개 모드로 지정합니다.

figure
modalsd(frf,f,fs,'MaxModes',30);

플롯을 확대합니다. 평균 응답 함수가 시스템의 물리적 주파수에 해당하는 373Hz, 852Hz, 1371Hz에서 최댓값을 가집니다. 최댓값을 변수에 저장합니다.

phfr = [373 852 1371];

최소제곱 복소수 지수(LSCE) 알고리즘을 사용하여 모드 파라미터를 계산합니다. 모델 차수를 6개 모드로 지정하고, 안정화 다이어그램에서 확인된 3개 모드의 물리적 주파수를 지정합니다. 아래 함수는 각 입력값을 기준으로 고유 주파수와 감쇠비를 하나의 세트로 생성합니다.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,6,'PhysFreq',phfr);

복원한 주파수-응답 함수를 플로팅하고 원래 함수와 비교합니다.

for k = 1:2
    for m = 1:3
        subplot(2,3,m+3*(k-1))
        plot(f/1000,10*log10(abs(frf(:,m,k))))
        hold on
        plot(f/1000,10*log10(abs(ofrf(:,m,k))))
        hold off
        text(1,-50,[['Output ';' Input '] num2str([m k]')])
        ylim([-100 -40])
    end
end
subplot(2,3,2)
title('Frequency-Response Functions')

입력 인수

frf — 주파수-응답 함수

벡터 | 행렬 | 3차원 배열

주파수-응답 함수로, 벡터, 행렬 또는 3차원 배열로 지정됩니다. frf는 p×m×n의 크기를 갖습니다. 여기서 p는 주파수 Bin의 개수이고, m은 응답 신호의 개수이고, n은 전달 함수 추정에 사용되는 가진 신호의 개수입니다. frf는 동적 유연성(리셉턴스) 형식이라고 가정합니다.

modalfrf를 사용하여 측정 데이터에서 주파수-응답 함수의 행렬을 생성합니다.

예제: 비감쇠 조화 진동자

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

f — 주파수

벡터

주파수로, 벡터로 지정됩니다. f의 요소 개수는 frf의 행 개수와 같아야 합니다.

데이터형: single | double

fs — 측정 데이터의 샘플 레이트

양의 스칼라

측정 데이터의 샘플 레이트로, 헤르츠를 단위로 하는 양의 스칼라로 지정됩니다.

데이터형: single | double

mnum — 모드 개수

양의 정수

모드 개수로, 양의 정수로 지정됩니다.

데이터형: single | double

sys — 식별된 시스템

식별된 파라미터를 가진 모델

식별된 시스템으로, 식별된 파라미터를 가진 모델로 지정됩니다. ssest (System Identification Toolbox), n4sid (System Identification Toolbox) 또는 tfest (System Identification Toolbox) 같은 추정 명령을 사용하여, 측정된 주파수-응답 함수 또는 시간 영역의 입력 신호와 출력 신호에서 sys를 생성합니다. 이에 대한 예제는 Modal Analysis of Identified Models 항목을 참조하십시오. 이 입력 인수를 사용하려면 System Identification Toolbox 라이선스가 있어야 합니다.

예: idss([0.5418 0.8373;-0.8373 0.5334],[0.4852;0.8373],[1 0],0,[0;0],[0;0],1)은 식별된 상태공간 모델을 생성합니다. 생성된 모델은 단위 탄성 상수를 갖는 스프링과 상수 0.01을 갖는 댐퍼로 벽에 연결된 단위 질량에 해당합니다. 질량 변위는 1Hz로 샘플링됩니다.

예: idtf([0 0.4582 0.4566],[1 -1.0752 0.99],1)은 식별된 전달 함수 모델을 생성합니다. 생성된 모델은 단위 탄성 상수를 갖는 스프링과 상수 0.01을 갖는 댐퍼로 벽에 연결된 단위 질량에 해당합니다. 질량 변위는 1Hz로 샘플링됩니다.

이름-값 인수

선택적 인수 쌍을 Name1=Value1,...,NameN=ValueN으로 지정합니다. 여기서 Name은 인수 이름이고 Value는 대응값입니다. 이름-값 인수는 다른 인수 뒤에 와야 하지만, 인수 쌍의 순서는 상관없습니다.

R2021a 이전 릴리스에서는 쉼표를 사용하여 각 이름과 값을 구분하고 Name을 따옴표로 묶으십시오.

예: 'FitMethod','pp','FreqRange',[0 500]은 피크 추출(peak-picking) 방법을 사용하여 피팅을 수행하고, 주파수 범위를 0Hz와 500Hz 사이로 제한합니다.

Feedthrough — 추정된 전달 함수에 피드스루가 있는지 여부

false (디폴트 값) | true

추정된 전달 함수에 피드스루가 있는지 여부로, 논리값으로 지정됩니다. 이 인수는 'FitMethod'가 'lsrf'로 지정된 경우에만 사용할 수 있습니다.

데이터형: logical

FitMethod — 피팅 알고리즘

'lsce' (디폴트 값) | 'lsrf' | 'pp'

피팅 알고리즘으로, 'lsce', 'lsrf' 또는 'pp'로 지정됩니다.

• 'lsce' — 최소제곱 복소수 지수 방법. 'lsce'를 지정할 경우, fn은 frf 크기에 관계없이 mnum개 요소로 구성된 벡터입니다.

• 'lsrf' — 최소제곱(Least-squares) 유리 함수 추정 방법. 'lsrf'를 지정할 경우, fn은 frf 크기에 관계없이 mnum개 요소로 구성된 벡터입니다. 이 방법은 [3] 항목에 설명되어 있습니다. 자세한 내용은 Continuous-Time Transfer Function Estimation Using Continuous-Time Frequency-Domain Data (System Identification Toolbox) 항목을 참조하십시오. 일반적으로 이 알고리즘은 비모수적 방법을 사용하는 구문보다 데이터를 더 적게 필요로 하며, f가 균일하지 않을 때 사용할 수 있는 유일한 방법입니다.

• 'pp' — 피크 추출 방법. frf가 n개의 가진 신호와 m개의 응답 신호에서 계산된 경우 fn은 mnum×m×n 배열입니다. 이 배열은 각 frf의 fn 추정값 한 개와 dr 추정값 한 개로 구성됩니다.

FreqRange — 주파수 범위

증가하는 양의 값으로 구성된 요소를 2개 가진 벡터

주파수 범위로, f에 지정된 범위 내에 포함된 증가하는 양수 값 요소 2개로 구성된 벡터로 지정됩니다.

데이터형: single | double

PhysFreq — 물리적 모드의 고유 주파수

벡터

분석에 포함할 물리적 모드의 고유 주파수로, f에 의해 구간 지정된 범위 내의 주파수 값으로 구성된 벡터로 지정됩니다. 함수는 해당 모드와 함께, 벡터에 지정된 값에 가장 근접한 고유 주파수를 분석에 포함합니다. 벡터에 m개의 주파수 값이 포함되어 있으면, fn과 dr은 각각 m개의 행을 갖고 ms는 m개의 열을 갖습니다. 이 인수를 지정하지 않을 경우, f의 전체 주파수 범위가 함수에 사용됩니다.

데이터형: single | double

DriveIndex — 구동점 주파수-응답 함수의 인덱스

[1 1] (디폴트 값) | 양의 정수로 구성된 요소를 2개 가진 벡터

구동점 주파수-응답 함수의 인덱스로, 양의 정수로 구성된 요소를 2개 가진 벡터로 지정됩니다. 벡터의 첫 번째 요소는 시스템 응답 수보다 작거나 같아야 합니다. 벡터의 두 번째 요소는 시스템 가진 수보다 작거나 같아야 합니다. 모드 형상은 구동점을 기반으로 하여 단위 모드로 정규화됩니다.

예: 'DriveIndex',[2 3]은 구동점 주파수-응답 함수가 frf(:,2,3)임을 지정합니다.

데이터형: single | double

출력 인수

fn — 고유 주파수

행렬 | 3차원 배열

고유 주파수로, 행렬 또는 3차원 배열로 반환됩니다. fn의 크기는 다음과 같이 'FitMethod'에 지정된 피팅 알고리즘 선택에 따라 달라집니다.

• 'lsce' 또는 'lsrf'를 지정할 경우, fn은 frf 크기에 관계없이 mnum개 요소로 구성된 벡터입니다. 시스템이 mnum개가 넘는 진동 모드를 갖는 경우, 'lsrf' 방법은 최소 감쇠 모드 중 고유 주파수가 낮은 순으로 처음 mnum개를 반환합니다.

• 'pp'를 지정할 경우, fn은 mnum×m×n 크기의 배열입니다. 이 배열은 각 frf의 fn 추정값 한 개와 dr 추정값 한 개로 구성됩니다.

dr — 감쇠비

행렬 | 3차원 배열

fn의 고유 주파수에 대한 감쇠비로, fn과 크기가 같은 행렬 또는 3차원 배열로 반환됩니다.

ms — 모드 형상 벡터

행렬

모드 형상 벡터로, 행렬로 반환됩니다. ms는 mnum개의 열을 갖고 각 열에는 길이가 q인 모드 형상 벡터가 포함되어 있습니다. 여기서 q는 가진 채널 수와 응답 채널 수 중 더 큰 값입니다.

ofrf — 복원한 주파수-응답 함수

벡터 | 행렬 | 3차원 배열

복원한 주파수-응답 함수로, frf와 크기가 같은 벡터, 행렬 또는 3차원 배열로 반환됩니다.

알고리즘

최소제곱 복소수 지수 방법

최소제곱 복소수 지수(Least-squares Complex Exponential) 방법은 각 주파수-응답 함수에 대응하는 임펄스 응답을 계산하고 Prony 방법을 사용하여 일련의 복소수 감쇠 정현파를 그 응답에 피팅합니다.

샘플링된 감쇠 정현파는 다음 형식으로 표현될 수 있습니다.

$$\begin{array}{c}{s}_i(n)=A_i{e}^{-b_{i}n/f_{\text{s}}}\cos (2\pi f_{i}n/f_{\text{s}}+{\varphi }_i)\\ =\frac{1}{2}{A}_i{e}^{j{\varphi }_i}\exp \left(-(b_i/f_{\text{s}}-j2\pi f_i/f_{\text{s}})n\right)+\frac{1}{2}{A}_i{e}^{-j{\varphi }_i}\exp \left(-(b_i/f_{\text{s}}+j2\pi f_i/f_{\text{s}})n\right)\\ \equiv a_{i+}{x}_{i+}^n+a_{i-}{x}_{i-}^n,\end{array}$$

여기서,

• f_s는 샘플 레이트입니다.

• f_i는 정현파 주파수입니다.

• b_i는 감쇠 계수입니다.

• A_i와 ϕ_i는 정현파의 진폭과 위상입니다.

a_i는 진폭이라고 하고, x_i는 극점이라고 합니다. Prony 방법은 샘플링된 함수 h(n)을 N/2개 모드가 중첩된 형태(N개의 진폭과 극점)로 표현합니다.

$$\begin{array}{c}h(0)=a_1{x}_1^0+a_2{x}_2^0\cdots +a_N{x}_N^0\\ h(1)=a_1{x}_1^1+a_2{x}_2^1+\cdots +a_N{x}_N^1\\ ⋮\\ h(N-1)=a_1{x}_1^{N-1}+a_2{x}_2^{N-1}+\cdots +a_N{x}_N^{N-1}.\end{array}$$

극점은 계수 c₀, c₁, …, c_N–1을 갖는 다항식의 근입니다.

$$x_i^N+c_{N-1}{x}_i^{N-1}+\cdots +c_1{x}_i^1+c_0{x}_i^0=0.$$

계수는 h에 대한 자기회귀 모델을 사용하여 구합니다(L = 2N개 샘플).

$$\left[\begin{array}{cccc}h(0)& h(1)& \cdots & h(N-1)\\ h(1)& h(2)& \cdots & h(N)\\ ⋮& ⋮& \ddots & ⋮\\ h(L-N-1)& h(L-N)& \cdots & h(L-2)\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ c_{N-1}\end{array}\right]=-\left[\begin{array}{c}h(N)\\ h(N+1)\\ ⋮\\ h(L-1)\end{array}\right].$$

극점은 roots 함수를 사용하여 구합니다. 극점을 알게 되면, 극점 로그의 허수부와 실수부를 파악하여 주파수와 감쇠 인자를 확인할 수 있습니다. 마지막 단계는 다음 식을 사용하여 진폭을 구하고 임펄스 응답을 복원하는 것입니다.

$$\left[\begin{array}{c}h(0)\\ ⋮\\ h(N-1)\end{array}\right]=\left[\begin{array}{ccc}{x}_1^0& \cdots & x_N^0\\ ⋮& \ddots & ⋮\\ x_1^{N-1}& \cdots & x_N^{N-1}\end{array}\right]\left[\begin{array}{c}{a}_1\\ ⋮\\ a_N\end{array}\right].$$

다음의 간단한 MATLAB^® 구문에 이 절차가 요약되어 있습니다.

N = 4;
L = 2*N;
h = rand(L,1);
c = hankel(h(1:N),h(L-N:L-1))\-h(N+1:L);
x = roots([1;c(N:-1:1)]).';
p = log(x);
hrec = x.^((0:L-1)')*(x.^((0:L-1)')\h(1:L));
sum(h-hrec)

ans =

   3.2613e-15 - 1.9297e-16i

이 시스템은 또한 여러 주파수-응답 함수의 샘플을 포함하도록 구성할 수 있으며 최소제곱을 사용하여 풀 수 있습니다.

피크 추출 방법

피크 추출(peak-picking) 방법에서는 주파수-응답 함수에서 유의미한 각 피크가 정확히 하나의 고유 모드에 해당하는 것으로 가정합니다. 피크 가까이에서 시스템이 자유도 1의 감쇠 조화 진동자처럼 동작하는 것으로 가정됩니다.

$$H(f)=\frac{-1}{{(2\pi )}^2}\frac{1/m}{f^2-j2{\zeta }_{\text{r}}{f}_{\text{r}}f-f_{\text{r}}^2}\Rightarrow f_{\text{r}}^{2}H(f)+j2{\zeta }_{\text{r}}{f}_{\text{r}}fH(f)-\frac{1}{{(2\pi )}^{2}m}=f^{2}H(f),$$

여기서 H는 주파수-응답 함수이고, f_r은 비감쇠 공명 주파수, ζ_r = b/(4mk)^1/2은 상대적 감쇠, b는 감쇠 상수, k는 탄성 상수, m은 질량입니다.

f_p에 있는 피크가 주어진 경우, 이 절차에서는 해당 피크와 함께 그 양쪽으로 고정된 개수만큼의 점을 취하고, 질량 항을 가변수 d로 대체하고, 다음 연립방정식을 풀어 모드 파라미터를 계산합니다.

$$\left[\begin{array}{ccc}H(f_{p-k})& j2f_{p-k}H(f_{p-k})& -1\\ ⋮& ⋮& ⋮\\ H(f_p)& j2f_{p}H(f_p)& -1\\ ⋮& ⋮& ⋮\\ H(f_{p+k})& j2f_{p+k}H(f_{p+k})& -1\end{array}\right]\left[\begin{array}{c}{f}_{\text{r}}^2\\ {\zeta }_{\text{r}}{f}_{\text{r}}\\ d\end{array}\right]=\left[\begin{array}{c}{f}_{p-k}^{2}H(f_{p-k})\\ ⋮\\ f_p^{2}H(f_p)\\ ⋮\\ f_{p+k}^{2}H(f_{p+k})\end{array}\right].$$

참고 문헌

[1] Allemang, Randall J., and David L. Brown. “Experimental Modal Analysis and Dynamic Component Synthesis, Vol. III: Modal Parameter Estimation.” Technical Report AFWAL-TR-87-3069. Air Force Wright Aeronautical Laboratories, Wright-Patterson Air Force Base, OH, December 1987.

[2] Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.

[3] Ozdemir, Ahmet Arda, and Suat Gumussoy. "Transfer Function Estimation in System Identification Toolbox via Vector Fitting." Proceedings of the 20th World Congress of the International Federation of Automatic Control, Toulouse, France, July 2017.

버전 내역

R2017a에 개발됨

참고 항목

modalfrf | modalsd | n4sid (System Identification Toolbox) | tfest (System Identification Toolbox) | tfestimate

도움말 항목

• Modal Analysis of Identified Models
• System Identification Overview (System Identification Toolbox)
• System Identification Workflow (System Identification Toolbox)
• Supported Continuous- and Discrete-Time Models (System Identification Toolbox)