forked from xmameisu/MorphodynamicProfiling
-
Notifications
You must be signed in to change notification settings - Fork 0
/
hilbertHuangTransform.m
69 lines (57 loc) · 2.18 KB
/
hilbertHuangTransform.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
function [instAmp,instFreq,hht,imf] = hilbertHuangTransform(TS,Fs,specRes)
%This function calculates the Hilbert spectrum using empirical
%mode decomposition as preliminary step
%
%USAGE:
% [instAmp,instFreq,hht] = hilbertHuangTransform(TS,Fs,specRes)
%
%Input:
% TS - time series vector
% Fs - sampling frequency; real number in Hz
% specRes - interger that specifies the number of frequency points in
% the final spectrum [frequency axis goes from 0 to Fs/2 with specRes number of points]
%
%Output:
% instAmp - instantaneous Amplitute for each imf
% instFreq - instantaneous frequency for each imf
% hht - structure: hht.timeVector [in seconds]
% hht.freqVector [in Hz]
% hht.spectrum
%
%Ref:
% book: Hilbert Transform and Application in Mechanical Vibration; Michael Feldman
%
%Marco Vilela, 2013; Xiao Ma, 2017
%Decomposing signal into intrinsic modes
imf = emd(TS);
%Initialization
[nImf,nObs] = size(imf);
instAmp = nan(nImf,nObs);
instFreq = nan(nImf,nObs-1);
spectrum_revised = nan(specRes,nObs-1);
for i=1:nImf
[instAmp(i,:),instFreq(i,:)] = getInstantaneousFrequency(imf(i,:));
end
instAmp(:,1) = [];%Make it same size as instFreq
%Guarantee frequencies within the physical range
instFreq = max(instFreq,0);
instFreq = min(instFreq,0.5);
% 0.5 is the Nyquist frequency for 1 unit of frequency. [No physical units for time]
%Constructing the time-frequency spectrum
dFreq = 0.5/(specRes - 1);
freqIdx = round( instFreq./dFreq ) + 1;
%finalFreqIdx = sub2ind(size(spectrum),freqIdx,repmat(1:size(spectrum,2),size(freqIdx,1),1));
%spectrum(finalFreqIdx) = instAmp;
for jj=1:nObs-1
for ii=1:specRes
ampid=find(freqIdx(:,jj)==ii);
spectrum_revised(ii,jj)=sum(instAmp(ampid,jj));
[rowfind0,columnfind0]=find(spectrum_revised(:,jj)==0);
spectrum_revised(rowfind0,jj) = NaN;
end
end
%Adding unit to the instantaneous frequency [Hz]
instFreq = Fs*instFreq;
hht.timeVector = (1/Fs)*(0:(nObs - 1));
hht.freqVector = Fs*(0:dFreq:0.5);
hht.spectrum_revised = spectrum_revised;