hht

힐베르트-황 변환(Hilbert-Huang Transform)

구문

hs = hht(imf)

hs = hht(imf,fs)

[hs,f,t] = hht(___)

[hs,f,t,imfinsf,imfinse] = hht(___)

[___] = hht(___,Name=Value)

hht(___)

hht(___,freqlocation)

설명

hs = hht(imf)는 내재 모드 함수 imf에 의해 지정된 신호의 힐베르트 스펙트럼 hs를 반환합니다. hs는 시간에 따라 스펙트럼 성분이 변하는 혼합 신호로 구성된 신호를 분석하는 데 유용합니다. hht를 사용하여 신호에 대한 힐베르트 스펙트럼 분석을 수행함으로써 국소화된 특징을 확인할 수 있습니다.

hs = hht(imf,fs)는 fs 레이트로 샘플링된 신호의 힐베르트 스펙트럼 hs를 반환합니다.

[hs,f,t] = hht(___)는 hs 외에 주파수 벡터 f와 시간 벡터 t도 반환합니다. 이 출력 인수는 위에 열거된 입력 구문과 함께 사용할 수 있습니다.

[hs,f,t,imfinsf,imfinse] = hht(___)는 신호 진단에 사용할, 내재 모드 함수의 순시 주파수 imfinsf와 순시 에너지 imfinse도 반환합니다.

[___] = hht(___,Name=Value)는 하나 이상의 이름-값 인수로 지정된 추가 옵션을 사용하여 힐베르트 스펙트럼 파라미터를 추정합니다.

hht(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 힐베르트 스펙트럼을 플로팅합니다. 이 구문은 위에 열거된 구문의 모든 입력 인수와 함께 사용할 수 있습니다.

hht(___,freqlocation)은 주파수 축의 위치를 지정하는 선택적 freqlocation 인수를 사용하여 힐베르트 스펙트럼을 플로팅합니다. 기본적으로 주파수는 y축에 표시됩니다.

예제

2차 처프의 힐베르트 스펙트럼

가우스 변조 2차 처프를 생성합니다. 2kHz의 샘플 레이트와 2초의 신호 지속 기간을 지정합니다.

fs = 2000;
t = 0:1/fs:2-1/fs;
q = chirp(t-2,4,1/2,6,'quadratic',100,'convex').*exp(-4*(t-1).^2);
plot(t,q)

emd를 사용하여 내재 모드 함수(IMF)와 잔차를 시각화합니다.

emd(q)

신호의 IMF를 계산합니다. 'Display' 이름-값 쌍을 사용하여 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준을 보여주는 테이블을 출력합니다.

imf = emd(q,'Display',1);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |    0.0063952   |  SiftMaxRelativeTolerance
      2      |        2     |       0.1007   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01189   |  SiftMaxRelativeTolerance
      4      |        2     |    0.0075124   |  SiftMaxRelativeTolerance
Decomposition stopped because the number of extrema in the residual signal is less than the 'MaxNumExtrema' value.

계산된 IMF를 사용하여 2차 처프의 힐베르트 스펙트럼을 플로팅합니다. 주파수 범위를 0Hz ~ 20Hz로 제한합니다.

hht(imf,fs,'FrequencyLimits',[0 20])

경험적 모드 분해 수행 및 신호의 힐베르트 스펙트럼 시각화하기

뚜렷한 주파수 변화를 보이는 정현파로 구성된 비정상 연속 신호를 불러오고 시각화합니다. 비정상 연속 신호의 예에는 착암기의 진동이나 폭죽 소리가 있습니다. 이 신호는 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");

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

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

hht(imf,fs)

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

흰긴수염고래 노래의 힐베르트 스펙트럼

태평양 흰긴수염고래의 오디오 데이터를 4kHz로 샘플링한 파일을 불러옵니다. 이 파일은 코넬대 생물 음향학 연구 프로그램(Cornell University Bioacoustics Research Program)에서 관리하는 동물 소리 라이브러리에서 생성된 것입니다. 피치를 올려서 울음소리가 더 잘 들리도록 데이터의 시간 스케일을 1/10로 압축했습니다. 신호를 MATLAB® 타임테이블로 변환하고 플로팅합니다. 이 신호의 잡음은 네 가지 특징이 두드러집니다. 첫 번째 소리는 짧게 반복되는 소리이고 나머지 세 번의 소리는 긴 울음소리입니다.

[w,fs] = audioread('bluewhale.wav');
whale = timetable(w,'SampleRate',fs);
stackedplot(whale);

emd를 사용하여 첫 3개의 내재 모드 함수(IMF)와 잔차를 시각화합니다.

emd(whale,'MaxNumIMF',3)

신호의 처음 3개 IMF를 계산합니다. 'Display' 이름-값 쌍을 사용하여 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준을 보여주는 테이블을 출력합니다.

imf = emd(whale,'MaxNumIMF',3,'Display',1);

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |      0.13523   |  SiftMaxRelativeTolerance
      2      |        2     |     0.030198   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01908   |  SiftMaxRelativeTolerance
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

계산된 IMF를 사용하여 신호의 힐베르트 스펙트럼을 플로팅합니다. 주파수 범위를 0Hz ~ 1400Hz로 제한합니다.

hht(imf,'FrequencyLimits',[0 1400])

동일한 주파수 범위에 대해 힐베르트 스펙트럼을 계산합니다. 짧게 반복되는 소리와 긴 울음소리의 힐베르트 스펙트럼을 메시 플롯으로 시각화합니다.

[hs,f,t] = hht(imf,'FrequencyLimits',[0 1400]);

mesh(seconds(t),f,hs,'EdgeColor','none','FaceColor','interp')
xlabel('Time (s)')
ylabel('Frequency (Hz)')
zlabel('Instantaneous Energy')

신호의 힐베르트 스펙트럼 파라미터 계산하기

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

load("sinusoidalSignalExampleData.mat","X","fs")
t = (0:length(X)-1)/fs;

plot(t,X)
xlabel("Time (s)")

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

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

[imf,residual,info] = emd(X,Interpolation="pchip");

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

힐베르트 스펙트럼 파라미터인 힐베르트 스펙트럼 hs, 주파수 벡터 f, 시간 벡터 t, 순시 주파수 imfinsf, 순시 에너지 imfinse를 계산합니다.

[hs,f,t,imfinsf,imfinse] = hht(imf,fs);

계산된 힐베르트 스펙트럼 파라미터를 시간-주파수 분석과 신호 진단에 사용합니다.

다중 성분 신호의 VMD

주파수 2Hz, 10Hz, 30Hz의 정현파 3개로 구성된 다중 성분 신호를 생성합니다. 정현파는 2초 동안 1kHz로 샘플링됩니다. 이 신호에 분산이 0.01²인 백색 가우스 잡음을 포함시킵니다.

fs = 1e3;
t = 1:1/fs:2-1/fs;
x = cos(2*pi*2*t) + ...
    2*cos(2*pi*10*t) + ...
    4*cos(2*pi*30*t) + ...
    0.01*randn(1,length(t));

잡음이 있는 신호의 IMF를 계산하고 3차원 플롯에 시각화합니다.

imf = vmd(x);
[p,q] = ndgrid(t,1:size(imf,2));
plot3(p,q,imf)
grid on
xlabel('Time Values')
ylabel('Mode Number')
zlabel('Mode Amplitude')

계산된 IMF를 사용하여 다중 성분 신호의 힐베르트 스펙트럼을 플로팅합니다. 주파수 범위를 [0, 40]Hz로 제한합니다.

hht(imf,fs,'FrequencyLimits',[0,40])

진동 신호의 힐베르트 스펙트럼 계산

손상된 베어링의 진동 신호를 시뮬레이션합니다. 이 신호의 힐베르트 스펙트럼을 계산하고 결함을 찾습니다.

베어링의 피치 직경은 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=\frac{1}{2}n{f}_0[1-\frac{d}{p}\cos \theta ],$$

$$f_0$$은 구동 속도, $$n$$은 구름 요소의 수, $$d$$는 구름 요소의 직경, $$p$$는 베어링의 피치 직경, $$\theta$$는 베어링 접촉각입니다. 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/150^2$$로 지정합니다.

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)를 계산합니다. 'Display' 이름-값 인수를 사용하여 각 IMF에 대한 선별 반복 횟수, 상대 허용오차, 선별 중지 기준을 보여주는 테이블을 출력합니다.

imfGood = emd(yGood,MaxNumIMF=5,Display=1);

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
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

출력 인수 없이 emd를 사용하여 첫 3개 IMF와 잔차를 시각화합니다.

emd(yGood,MaxNumIMF=5)

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

imfBad = emd(yBad,MaxNumIMF=5,Display=1);

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
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

emd(yBad,MaxNumIMF=5)

결함 있는 베어링 신호의 첫 번째 경험적 모드의 힐베르트 스펙트럼을 플로팅합니다. 첫 번째 모드는 고주파수 충격의 효과를 캡처합니다. 베어링 마모가 진행되면 충격 에너지가 커집니다.

figure
hht(imfBad(:,1),fs)

세 번째 모드의 힐베르트 스펙트럼은 진동 신호의 공명을 보여줍니다. 주파수 범위를 0Hz ~ 100Hz로 제한합니다.

hht(imfBad(:,3),fs,FrequencyLimits=[0 100])

비교를 위해 정상인 베어링 신호의 첫 번째 및 세 번째 모드의 힐베르트 스펙트럼을 플로팅합니다.

subplot(2,1,1)
hht(imfGood(:,1),fs)
subplot(2,1,2)
hht(imfGood(:,3),fs,FrequencyLimits=[0 100])

입력 인수

imf — 내재 모드 함수

행렬 | timetable형

내재 모드 함수로, 행렬 또는 timetable형으로 지정됩니다. imf는 극값 및 영점교차의 개수와 최대 하나만큼 차이 나는 신호로, 포락선이 0을 기준으로 대칭입니다. emd는 복잡한 신호를 분해한 다음, 힐베르트 스펙트럼 분석을 수행하는 데 필요한 유한개의 내재 모드 함수로 단순화하는 데 사용됩니다.

hht는 imf의 각 열을 내재 모드 함수로 취급합니다. imf 계산 방법에 대한 자세한 내용은 emd 항목을 참조하십시오.

fs — 샘플 레이트

2π (디폴트 값) | 양의 스칼라

샘플 레이트로, 양의 스칼라로 지정됩니다. fs가 제공되지 않으면 힐베르트 스펙트럼을 계산할 때 정규화 주파수 2π가 사용됩니다. imf가 타임테이블로 지정된 경우 이 인수에서 샘플 레이트가 추론됩니다.

freqlocation — 플롯에서의 주파수 축 위치

"yaxis" (디폴트 값) | "xaxis"

플롯에서의 주파수 축 위치로, "yaxis" 또는 "xaxis"로 지정됩니다. 플롯의 y축이나 x축에 주파수 데이터를 표시하려면 freqlocation을 각각 "yaxis" 또는 "xaxis"로 지정하십시오.

이름-값 인수

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

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

예: 'FrequencyResolution',1

FrequencyLimits — 힐베르트 스펙트럼을 계산하기 위한 주파수 제한

[0,fs/2] (디폴트 값) | 1×2 정수 값 벡터

힐베르트 스펙트럼을 계산하기 위한 주파수 제한으로, 1×2 정수 값 벡터로 지정됩니다. FrequencyLimits는 Hz 단위로 지정됩니다.

FrequencyResolution — 주파수 범위를 이산화하기 위한 주파수 분해능

(f_high-f_low)/100 (디폴트 값) | 양의 스칼라

주파수 제한을 이산화하기 위한 주파수 분해능으로, 양의 스칼라로 지정됩니다. FrequencyResolution은 Hz 단위로 지정됩니다. FrequencyResolution이 지정되지 않은 경우 (f_high-f_low)/100의 값은 FrequencyLimits에서 추론됩니다. 여기서 f_high는 FrequencyLimits의 상한이고, f_low는 하한입니다.

MinThreshold — 힐베르트 스펙트럼의 최소 임계값

-inf (디폴트 값) | 스칼라

힐베르트 스펙트럼의 최소 임계값으로, 스칼라로 지정됩니다. MinThreshold는 $$10{\log }_{10}\left(hs\right)$$의 대응하는 요소가 MinThreshold보다 작은 경우 hs의 요소를 0으로 설정합니다.

출력 인수

hs — 신호의 힐베르트 스펙트럼

희소 행렬

신호의 힐베르트 스펙트럼으로, 희소 행렬로 반환됩니다. 시간-주파수 분석이나 신호의 국소화된 특징을 확인하는 데 hs를 사용합니다.

f — 주파수 값

벡터

신호의 주파수 값으로, 벡터로 반환됩니다. hht는 주파수 벡터 f와 시간 벡터 t를 사용하여 힐베르트 스펙트럼 플롯을 생성합니다.

f는 수학적으로 f = f_low : f_res : f_high와 같이 나타냅니다. 여기서 f_res는 주파수 분해능입니다.

t — 시간 값

벡터 | duration형 배열

신호의 시간 값으로, 벡터 또는 duration형 배열로 반환됩니다. hht는 시간 벡터 t와 주파수 벡터 f를 사용하여 힐베르트 스펙트럼 플롯을 생성합니다.

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

• 배열 - imf가 배열로 지정된 경우.

• duration형 배열 - imf가 균일하게 샘플링된 timetable형으로 지정된 경우.

imfinsf — 각 IMF의 순시 주파수

벡터 | 행렬 | timetable형

각 IMF의 순시 주파수로, 벡터, 행렬 또는 timetable형으로 반환됩니다.

imfinsf는 imf와 같은 개수의 열을 가지며 다음 형태로 반환됩니다.

• 벡터 - imf가 벡터로 지정된 경우.

• 행렬 - imf가 행렬로 지정된 경우.

• timetable형 - imf가 균일하게 샘플링된 timetable형으로 지정된 경우.

imfinse — 각 IMF의 순시 에너지

벡터 | 행렬 | timetable형

각 IMF의 순시 에너지로, 벡터, 행렬 또는 timetable형으로 반환됩니다.

imfinse는 imf와 같은 개수의 열을 가지며 다음 형태로 반환됩니다.

• 벡터 - imf가 벡터로 지정된 경우.

• 행렬 - imf가 행렬로 지정된 경우.

• timetable형 - imf가 균일하게 샘플링된 timetable형으로 지정된 경우.

알고리즘

힐베르트-황 변환은 비정상(Nonstationary) 비선형 데이터의 시간-주파수 분석을 수행하는 데 유용합니다. 힐베르트-황 절차는 다음 단계로 구성됩니다.

1. emd 또는 vmd가 데이터 세트 x를 유한 개의 내재 모드 함수로 분해합니다.

2. 각 내재 모드 함수 x_i에 대해 함수 hht가 다음을 수행합니다.

   1. hilbert를 사용하여 해석적 신호 $$z_i(t)=x_i(t)+jH\{x_i(t)\}$$를 계산합니다. 여기서 H{x_i}는 x_i의 힐베르트 변환입니다.

   2. z_i를 $$z_i(t)=a_i(t)e^{j{\theta }_i(t)}$$으로 표현합니다. 여기서 a_i(t)는 순시 진폭이고 $${\theta }_i(t)$$는 순시 위상입니다.

   3. 순시 에너지 $$|a_i(t){|}^2$$과 순시 주파수 $${\omega }_i(t)\equiv d{\theta }_i(t)/dt$$를 계산합니다. 샘플 레이트가 주어진 경우 hht는 $${\omega }_i(t)$$를 Hz 단위의 주파수로 변환합니다.

   4. 순시 에너지를 imfinse에 출력하고 순시 주파수를 imfinsf에 출력합니다.

3. 출력 인수 없이 호출될 경우 hht는 진폭에 비례하는 색을 사용하여 신호의 에너지를 시간과 주파수의 함수로 플로팅합니다.

참고 문헌

[1] Huang, Norden E, and Samuel S P Shen. Hilbert–Huang Transform and Its Applications. 2nd ed. Vol. 16. Interdisciplinary Mathematical Sciences. WORLD SCIENTIFIC, 2014. https://doi.org/10.1142/8804.

[2] Huang, Norden E., Zhaohua Wu, Steven R. Long, Kenneth C. Arnold, Xianyao Chen, and Karin Blank. “ON INSTANTANEOUS FREQUENCY.” Advances in Adaptive Data Analysis 01, no. 02 (April 2009): 177–229. https://doi.org/10.1142/S1793536909000096.

확장 기능

C/C++ 코드 생성

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

사용법 관련 참고 및 제한 사항:

• 이름-값 인수는 컴파일타임 상수여야 합니다.

버전 내역

R2018a에 개발됨

참고 항목

emd | hilbert | vmd

도움말 항목

• Time-Frequency Gallery
• 해석적 신호와 힐베르트 변환
• 힐베르트 변환과 순시 주파수