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

emd

경험적 모드 분해

설명

예제

[imf,residual] = emd(X)X의 경험적 모드 분해에 해당하는, 내재 모드 함수 imf와 잔차 신호 residual을 반환합니다. emd를 사용하여 복잡한 신호를 분해한 다음, 힐베르트 스펙트럼 분석을 수행하는 데 필요한 유한개의 내재 모드 함수(intrinsic mode function)로 단순화합니다.

예제

[imf,residual,info] = emd(X)는 진단을 목적으로 IMF와 잔차 신호에 대한 추가 정보 info를 반환합니다.

예제

[___] = emd(___,Name,Value)는 하나 이상의 Name,Value 쌍의 인수로 지정된 추가 옵션을 사용하여 emd를 추정합니다.

예제

emd(___)는 동일한 Figure에 원래 신호, IMF, 잔차 신호를 플로팅합니다.

예제

모두 축소

뚜렷한 주파수 변화를 보이는 정현파로 구성된 비정상 연속 신호를 불러오고 시각화합니다. 비정상 연속 신호의 예에는 착암기의 진동이나 폭죽 소리가 있습니다. 이 신호는 fs의 레이트로 샘플링됩니다.

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

혼합된 신호에는 서로 다른 진폭과 주파수 값을 갖는 정현파가 포함되어 있습니다.

힐베르트 스펙트럼 플롯을 생성하려면 신호의 내재 모드 함수(IMF)가 필요합니다. 경험적 모드 분해를 수행하여 신호의 IMF와 잔차를 계산합니다. 신호가 매끄럽지 않기 때문에 보간 방법으로 'pchip'을 지정합니다.

[imf,residual,info] = emd(X,'Interpolation','pchip');
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

명령 창에 생성된 테이블에는 생성된 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준이 나와 있습니다. 이 정보는 info에도 포함되어 있습니다. 'Display',0 이름 값 쌍을 추가하여 테이블을 숨길 수 있습니다.

경험적 모드 분해를 사용하여 얻은 imf 성분을 사용하여 힐베르트 스펙트럼 플롯을 생성합니다.

hht(imf,fs)

시간에 대한 주파수의 플롯은 IMF의 각 지점에서의 순시 에너지를 나타내는 세로 컬러바가 있는 희소 플롯입니다. 이 플롯은 혼합된 원래 신호에서 분해된 각 성분의 순시 주파수 스펙트럼을 나타냅니다. 플롯의 1초 시점에 주파수 변화가 뚜렷한 3개의 IMF가 표시됩니다.

다음 삼각 함수 항등식은 동일한 물리적 신호에 대한 서로 다른 2가지 표현을 보여줍니다.

52cos2πf1t+14(cos2π(f1+f2)t+cos2π(f1-f2)t)=(2+cos2πf2t)cos2πf1t.

2개의 정현파 sz를 생성합니다. s는 세 사인파의 합이고, z는 변조된 진폭의 단일 사인파입니다. 두 신호 간 차분의 무한대 노름을 계산하여 두 신호가 같음을 확인합니다.

t = 0:1e-3:10;
omega1 = 2*pi*100;
omega2 = 2*pi*20;
s = 0.25*cos((omega1-omega2)*t) + 2.5*cos(omega1*t) + 0.25*cos((omega1+omega2)*t);
z = (2+cos(omega2/2*t).^2).*cos(omega1*t);

norm(s-z,Inf) 
ans = 3.2729e-13

정현파를 플로팅하고 2초 지점에서 시작하는 1초 구간을 선택합니다.

plot(t,[s' z'])
xlim([2 3])
xlabel('Time (s)')
ylabel('Signal')

이 신호의 스펙트로그램을 계산합니다. 스펙트로그램은 서로 다른 3가지 정현파 성분을 보여줍니다. 푸리에 분석에서는 신호를 사인파가 중첩된 것으로 인식합니다.

pspectrum(s,1000,'spectrogram','TimeResolution',4)

emd를 사용하여 신호의 내재 모드 함수(IMF) 및 추가 진단 정보를 계산합니다. 이 함수는 기본적으로 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준을 보여주는 표를 출력합니다. 경험적 모드 분해에서는 이 신호를 z로 인식합니다.

[imf,~,info] = emd(s);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |   1.8025e-06   |  SiftMaxRelativeTolerance
The decomposition stopped because the current energy ratio is greater than 'MaxEnergyRatio'.

영점교차의 개수와 국소 극값은 최대 1만큼 차이가 납니다. 이는 신호가 IMF가 되는 데 필요한 조건을 충족합니다.

info.NumZerocrossing - info.NumExtrema
ans = 1

IMF를 플로팅하고 2초 지점에서 시작하는 0.5초 구간을 선택합니다. emd가 신호를 변조된 진폭으로 인식하기 때문에 IMF는 AM 신호입니다.

plot(t,imf)
xlim([2 2.5])
xlabel('Time (s)')
ylabel('IMF')

손상된 베어링의 진동 신호를 시뮬레이션합니다. 신호의 IMF를 시각화하고 결함이 있는지 찾아보기 위해 경험적 모드 분해를 수행합니다.

베어링의 피치 직경은 12cm이며 8개의 구름 요소가 있습니다. 각 구름 요소의 지름은 2cm입니다. 내륜이 초당 25 사이클로 구동되는 동안 외륜은 부동 상태로 유지됩니다. 가속도계가 10kHz로 베어링 진동을 샘플링합니다.

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

정상적인 베어링의 진동 신호에는 여러 차수의 구동 주파수가 포함되어 있습니다.

t = 0:1/fs:10-1/fs;
yHealthy = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

측정 과정의 중간에 베어링 진동에 공명이 가해집니다.

yHealthy = (1+1./(1+linspace(-10,10,length(yHealthy)).^4)).*yHealthy;

공명으로 인해 베어링의 외륜에 결함이 생겨 점점 더 마모됩니다. 이 결함 때문에 베어링의 외륜 통과 주파수(ball pass frequency outer race, BPFO)에서 일련의 충격이 발생합니다.

BPFO=12nf0[1-dpcosθ],

f0은 구동 속도, n은 구름 요소의 수, d는 구름 요소의 직경, p는 베어링의 피치 직경, θ는 베어링 접촉각입니다. 15°의 접촉각을 가정하고 BPFO를 계산합니다.

ca = 15;
bpfo = n*f0/2*(1-d/p*cosd(ca));

pulstran 함수를 사용하여 충격을 5밀리초 정현파로 구성된 주기적 열로 모델링합니다. 각 3kHz 정현파에는 플랫 탑 윈도우가 적용됩니다. 베어링 진동 신호에 점진적 마모를 적용하기 위해 멱법칙을 사용합니다.

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
wImpact = flattopwin(length(tImpact))'/10;
xImpact = sin(2*pi*fImpact*tImpact).*wImpact;

tx = 0:1/bpfo:t(end);
tx = [tx; 1.3.^tx-2];

nWear = 49000;
nSamples = 100000;
yImpact = pulstran(t,tx',xImpact,fs)/5;
yImpact = [zeros(1,nWear) yImpact(1,(nWear+1):nSamples)];

정상적인 신호에 충격을 더하여 BPFO 진동 신호를 생성합니다. 신호를 플로팅하고 5.0초 지점에서 시작하는 0.3초 구간을 선택합니다.

yBPFO = yImpact + yHealthy;

xLimLeft = 5.0;
xLimRight = 5.3;
yMin = -0.6;
yMax = 0.6;

plot(t,yBPFO)

hold on
[limLeft,limRight] = meshgrid([xLimLeft xLimRight],[yMin yMax]);
plot(limLeft,limRight,'--')
hold off

충격의 효과를 시각화하기 위해 선택된 구간을 확대합니다.

xlim([xLimLeft xLimRight])

백색 가우스 잡음을 신호에 추가합니다. 잡음 분산을 1/1502로 지정합니다.

rn = 150;
yGood = yHealthy + randn(size(yHealthy))/rn;
yBad = yBPFO + randn(size(yHealthy))/rn;

plot(t,yGood,t,yBad)
xlim([xLimLeft xLimRight])
legend('Healthy','Damaged')

emd를 사용하여 정상적인 베어링 신호의 경험적 모드 분해를 수행합니다. 처음 5개의 내재 모드 함수(IMF)를 계산합니다. 이 함수는 기본적으로 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준을 보여주는 표를 출력합니다.

imfGood = emd(yGood,'MaxNumIMF',5);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        3     |     0.017132   |  SiftMaxRelativeTolerance
      2      |        3     |      0.12694   |  SiftMaxRelativeTolerance
      3      |        6     |      0.14582   |  SiftMaxRelativeTolerance
      4      |        1     |     0.011082   |  SiftMaxRelativeTolerance
      5      |        2     |      0.03463   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

출력 인수 없이 emd를 사용하여 첫 3개 모드와 잔차를 시각화합니다. 'Display'0으로 설정하여 표를 숨깁니다.

emd(yGood,'MaxNumIMF',5,'Display',0)

결함 있는 베어링 신호의 IMF를 계산하고 시각화합니다. 첫 번째 경험적 모드에서 고주파수 영향이 나타납니다. 이 고주파수 모드는 마모가 진행되면서 에너지가 늘어납니다. 세 번째 모드는 진동 신호의 공명을 보여줍니다.

imfBad = emd(yBad,'MaxNumIMF',5);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.041274   |  SiftMaxRelativeTolerance
      2      |        3     |      0.16695   |  SiftMaxRelativeTolerance
      3      |        3     |      0.18428   |  SiftMaxRelativeTolerance
      4      |        1     |     0.037177   |  SiftMaxRelativeTolerance
      5      |        2     |     0.095861   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.
emd(yBad,'MaxNumIMF',5,'Display',0)

분석의 다음 단계는 추출된 IMF의 힐베르트 스펙트럼 계산입니다. 자세한 내용은 진동 신호의 힐베르트 스펙트럼 계산 예제를 참조하십시오.

뚜렷한 주파수 변화를 보이는 정현파로 구성된 비정상 연속 신호를 불러오고 시각화합니다. 비정상 연속 신호의 예에는 착암기의 진동이나 폭죽 소리가 있습니다. 이 신호는 fs의 레이트로 샘플링됩니다.

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

혼합된 신호에는 서로 다른 진폭과 주파수 값을 갖는 정현파가 포함되어 있습니다.

경험적 모드 분해를 수행하여 신호의 내재 모드 함수와 잔차를 플로팅합니다. 신호가 매끄럽지 않기 때문에 보간 방법으로 'pchip'을 지정합니다.

emd(X,'Interpolation','pchip')
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

원래 신호, 처음 3개의 IMF, 잔차가 포함된 대화형 플롯이 생성됩니다. 명령 창에 생성된 테이블에는 생성된 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준이 나와 있습니다. Display0으로 지정하여 테이블을 숨길 수 있습니다.

플롯에서 흰색 영역을 마우스 오른쪽 버튼으로 클릭하여, IMF 선택기 창을 엽니다. IMF 선택기를 사용하여 생성된 IMF, 원래 신호, 잔차를 선택적으로 볼 수 있습니다.

표시할 IMF를 목록에서 선택합니다. 플롯에 원래 신호와 잔차를 표시할지 여부를 선택합니다.

이제, 선택한 IMF가 플롯에 표시됩니다.

플롯을 사용하여, 원래 신호에서 분해된 각 성분을 잔차와 함께 시각화할 수 있습니다. 잔차는 총 IMF 개수에 대해 계산되며, IMF 선택기 창에서 선택한 IMF에 따라 변경되지 않습니다.

입력 인수

모두 축소

균일하게 샘플링된 시간 영역 신호로, 벡터 또는 데이터 열이 하나인 타임테이블로 지정됩니다.

이름-값 쌍의 인수

선택적으로 Name,Value 인수가 쉼표로 구분되어 지정됩니다. 여기서 Name은 인수 이름이고 Value는 대응값입니다. Name은 따옴표 안에 표시해야 합니다. Name1,Value1,...,NameN,ValueN과 같이 여러 개의 이름-값 쌍의 인수를 어떤 순서로든 지정할 수 있습니다.

예: 'MaxNumIMF',5

코시(Cauchy) 유형 수렴 기준으로, 'SiftRelativeTolerance'와 함께 양의 스칼라가 쉼표로 구분되어 지정됩니다. SiftRelativeTolerance는 선별 중지 기준 중 하나입니다. 즉, 현재 상대 허용오차가 SiftRelativeTolerance보다 작으면 선별이 중지됩니다.

선별 반복의 최대 횟수로, 'SiftMaxIterations'와 함께 양의 정수 스칼라가 쉼표로 구분되어 지정됩니다. SiftMaxIterations는 선별 중지 기준 중 하나입니다. 즉, 현재 반복 횟수가 SiftMaxIterations보다 크면 선별이 중지됩니다.

SiftMaxIterations는 양의 정수로만 지정할 수 있습니다.

추출된 최대 IMF 개수로, 'MaxNumIMF'와 함께 양의 정수 스칼라가 쉼표로 구분되어 지정됩니다. MaxNumIMF는 분해 중지 기준 중 하나입니다. 즉, 생성된 IMF 개수가 MaxNumIMF와 같으면 분해가 중지됩니다.

MaxNumIMF는 양의 정수로만 지정할 수 있습니다.

잔차 신호 극값의 최대 개수로, 'MaxNumExtrema'와 함께 양의 정수 스칼라가 쉼표로 구분되어 지정됩니다. MaxNumExtrema는 분해 중지 기준 중 하나입니다. 즉, 극값 개수가 MaxNumExtrema보다 적으면 분해가 중지됩니다.

MaxNumExtrema는 양의 정수로만 지정할 수 있습니다.

신호 대 잔차 에너지 비로, 'MaxEnergyRatio'와 함께 스칼라가 쉼표로 구분되어 지정됩니다. MaxEnergyRatio는 선별 시작 시의 신호 에너지 대 평균 포락선 에너지 비입니다. MaxEnergyRatio는 분해 중지 기준 중 하나로, 현재 에너지 비가 MaxEnergyRatio보다 크면 분해가 중지됩니다.

포락선 생성을 위한 보간 방법으로, 'Interpolation'과 함께 'spline' 또는 'pchip'이 쉼표로 구분되어 지정됩니다.

Interpolation은 다음과 같이 지정합니다.

  • 'spline' - X가 매끄러운 신호인 경우

  • 'pchip' - X가 매끄러운 신호가 아닌 경우

'spline' 보간 방법은 3차 스플라인을 사용하는 반면, 'pchip'은 조각별 3차 에르미트 보간 다항식 방법을 사용합니다.

명령 창에서의 정보 표시 전환으로, 'Display'와 함께 1 또는 0이 쉼표로 구분되어 지정됩니다. 명령 창에 생성된 테이블에는 생성된 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준이 나와 있습니다. 테이블을 표시하려면 Display를 1로 지정하고 테이블을 숨기려면 0으로 지정하십시오.

출력 인수

모두 축소

내재 모드 함수(intrinsic mode function)로, 행렬 또는 timetable형으로 반환됩니다. imf는 극값과 영점교차의 개수가 최대 1만큼 차이 나고 포락선이 0을 기준으로 대칭인 함수입니다. imf를 사용하면 힐베르트-황 변환(Hilbert-Huang Transform)을 적용하여 신호에 대한 스펙트럼 분석을 수행할 수 있습니다.

imf는 다음 형태로 반환됩니다.

  • X가 벡터인 경우 각 열이 imf인 행렬

  • X가 데이터 열이 하나인 타임테이블인 경우 타임테이블

신호의 잔차로, 열 벡터 또는 데이터 열이 하나인 타임테이블로 반환됩니다. residual은 원래 신호 X에서 emd로 분해되지 않은 부분을 나타냅니다.

residual은 다음 형태로 반환됩니다.

  • X가 벡터인 경우 열 벡터

  • X가 데이터 열이 하나인 타임테이블인 경우 단일 데이터 열의 타임테이블

진단에 대한 추가 정보로, 다음 필드가 있는 구조체로 반환됩니다.

  • NumIMF - 신호에서 추출된 IMF 개수

  • NumExtrema - 각 IMF의 극값 개수

  • NumZeroCrossing - 각 IMF의 영점교차 개수

  • NumSifting - 각 IMF에 대해 수행된 선별 횟수

  • MeanEnvelopeEnergy - 각 IMF 계산에서 획득한 상부 포락선과 하부 포락선의 평균이 갖는 에너지

  • RelativeTolerance - 각 IMF의 상대 허용오차

알고리즘

경험적 모드 분해

emd는 선별 과정을 사용하여 신호 X(t)k개의 내재 모드 함수(IMF)와 잔차 rk(t)로 분해합니다. [1][2]에 언급된 선별 과정에 대한 간략한 개요는 다음과 같습니다.

  1. 신호 X(t)에 대한 국소 최댓값과 국소 최솟값을 구하여 상부 포락선 s+(t) 및 하부 포락선 s(t)를 생성합니다.

  2. i번째 반복에 대한 평균 포락선 mk,i(t)를 계산합니다.

    mk,i(t)=12[s+(t)+s(t)]

  3. 첫 번째 반복에서 ck(t) = X(t)를 사용하여 잔차 신호에서 평균 포락선을 뺍니다.

    ck(t)=ck(t)mk,i(t)

    ck(t)가 IMF 기준과 일치하지 않을 경우 4단계와 5단계를 건너뜁니다. ck(t)의 새 값을 사용하여 절차가 1단계에서 다시 반복됩니다.

  4. ck(t)가 IMF 기준과 일치할 경우 새 잔차가 계산됩니다. 잔차 신호를 업데이트하려면 이전 잔차 신호에서 k번째 IMF를 빼십시오.

    rk(t)=rk1(t)ck(t)

  5. 그런 다음 새 신호 rk(t)로 획득한 잔차를 사용하여 1단계부터 시작한 후, ck(t)를 내재 모드 함수로 저장합니다.

N개의 내재 모드 함수가 있는 경우 원래 함수는 다음과 같이 표현됩니다.

X(t)=i=1NIMFi(t)+rN(t)

선별 과정에 대한 자세한 내용은 [1][2]를 참조하십시오.

SiftRelativeTolerance

SiftRelativeTolerance[4]에 제안된 코시 유형 중지 기준입니다. 현재 상대 허용오차가 SiftRelativeTolerance보다 작으면 선별 과정이 중지됩니다. 현재 상대 허용오차는 다음과 같이 정의됩니다.

Relative Tolerancerprev(t)rcur(t)22rprev(t)22.

코시 기준은 영점교차와 국소 극값의 개수를 직접 세지 않기 때문에, 분해에서 반환되는 IMF가 내재 모드 함수의 엄격한 정의를 충족하지 않을 가능성이 있습니다. 이 경우 SiftRelativeTolerance의 값을 디폴트 값에서 줄여 볼 수 있습니다. 중지 기준에 대한 자세한 내용은 [4] 항목을 참조하십시오. 참고 문헌에서는 경험적 모드 분해에서 엄격하게 정의된 IMF를 고집하는 것의 이점과 단점도 설명합니다.

MaxEnergyRatio

에너지 비는 선별 시작 시의 신호 에너지 대 평균 포락선 에너지의 비입니다[3]. 현재 에너지 비가 MaxEnergyRatio보다 크면 분해가 중지됩니다. k개의 IMF에 대한 EnergyRatio는 다음과 같이 정의됩니다.

Energy Ratio10log10(X(t)2ri(t)2).

참고 문헌

[1] Huang, Norden E., Zheng Shen, Steven R. Long, Manli C. Wu, Hsing H. Shih, Quanan Zheng, Nai-Chyuan Yen, Chi Chao Tung, and Henry H. Liu. "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis." Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences. Vol. 454, 1998, pp. 903–995.

[2] Rilling, G., Patrick Flandrin, and Paulo Gonçalves. "On empirical mode decomposition and its algorithms." IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing 2003. NSIP-03. Grado, Italy. 8–11.

[3] Rato, R.T., Manuel Ortigueira, and Arnaldo Batista. "On the HHT, its problems, and some solutions." Mechanical Systems and Signal Processing Vol. 22, 2008, pp. 1374–1394.

[4] Wang, Gang, Xian-Yao Chen, Fang-Li Qiao, Zhaohua Wu, and Norden Huang. "On Intrinsic Mode Function." Advances in Adaptive Data Analysis. Vol. 2, Number 3, 2010, pp. 277–293.

확장 기능

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

참고 항목

R2018a에 개발됨