Initial commit

Codes used in papers published since 2016 (CBEB'16, CBEB'18, WCMPBE'18, and EMBC'19).
NIATS-UFU · Aug 30, 2019 · 9ea87fd · 9ea87fd
1 parent 2ad08ba
commit 9ea87fd
Show file tree

Hide file tree

Showing 415 changed files with 89,305 additions and 0 deletions.
diff --git a/AudioOnset.m b/AudioOnset.m
@@ -0,0  1,34 @@
 % Function name....: AudioOnset
 % Date.............: May 04, 2016
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                    --
 % Parameters.......:
 %                    inputSignal ->
 %                    sampleFrequency ->
 %                    threshold ->
 %                    timeSkip -> in seconds
 % Return...........:
 %                    onsets -> indices of calculated onsets
 %
 % Remarks..........:
 
 function [onsets] = AudioOnset(inputSignal, sampleFrequency, threshold, timeSkip)
 onsets = [];
 
 if nargin < 4
     error('Insufficient arguments');
     return;
 end
 
 pointsSkip = timeSkip * sampleFrequency;
 onsets(1) = find(inputSignal > threshold, 1);
 nextIndex = onsets(1) pointsSkip;
 inputSignalLength = length(inputSignal);
 while nextIndex < inputSignalLength
     onsets(end 1) = (nextIndex   find(inputSignal(nextIndex:end) > threshold, 1));
     nextIndex = onsets(end)   pointsSkip;
 end
 
 end
 
diff --git a/Envelope.m b/Envelope.m
@@ -0,0  1,43 @@
 % Function name....: Envelope
 % Date.............: May 03, 2016
 % Mod date.........: Jan  12, 2019 (v2)
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                    Envelope estimates upper and lower envelopes of time-series.
 % Parameters.......:
 %                    inputSignal -> input time-series
 %                    (optional param) minPeakDistance -> minimum peak (and valley) separation
 %                    (optional param) windowing -> boolean flag to enable
 %                    finding peaks using the carrier frequency information
 % Return...........:
 %                    upperEnv -> upper envelope
 %                    upperEnvIndex -> upper envelope indexes
 %                    lowerEnv -> lower envelope
 %                    lowerEnvIndex -> lower envelope indexes
 % Remarks..........:
 %                    Envelope uses the PeakFinder (non matlab function) function to estimate the
 %                    local maxima and minima of the input signal. This
 %                    function do not use interpolation methods.
 function [upperEnv,upperEnvIndex,lowerEnv,lowerEnvIndex] = Envelope(inputSignal,minPeakDistance,windowing)
 upperEnv = [];
 upperEnvIndex = [];
 lowerEnv = [];
 lowerEnvIndex = [];
 
 if nargin == 1
     [upperEnv,upperEnvIndex] = PeakFinder(inputSignal);
     [lowerEnv,lowerEnvIndex] = PeakFinder(-inputSignal);
 elseif nargin == 2
     [upperEnv,upperEnvIndex] = PeakFinder(inputSignal, minPeakDistance);
     [lowerEnv,lowerEnvIndex] = PeakFinder(-inputSignal, minPeakDistance);
 elseif nargin == 3
     [upperEnv,upperEnvIndex] = PeakFinder(inputSignal, minPeakDistance, windowing);
     [lowerEnv,lowerEnvIndex] = PeakFinder(-inputSignal, minPeakDistance, windowing);
 else
     error('Insufficient arguments');
     return;
 end
 % Adjust lower envelope
 lowerEnv = -lowerEnv;
 end
 
diff --git a/GetActivityWn.m b/GetActivityWn.m
@@ -0,0  1,68 @@
 % Function name....: GetActivityWn
 % Date.............: May 19, 2019
 % Mod date.........: 
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                   Get activity windows based on y (trigger signal)
 % Parameters.......:
 %                   x -> time
 %                   y -> signal (trigger signal)
 %                   ref -> reference time-signal(adjust wnIndex for
 %                   different sample rate)
 % Return...........:
 %                   wnTime  -> Initial time and end time of each window
 %                   wnIndex -> Initial and end index of each window
 % Remarks..........:
 %                    
 function [wnTime, wnIndex] = GetActivityWn(x, y, refTime, plotFlag, humanReactionTime)
 % x = time;
 % y = digInA;
 % refTime = psEnvTime;
 % Invert signal if the initial state starts with 1 (high level)
 if y(1) == 1
     y = y * -1;
 end
 initialState = y(1);
 % figure;findpeaks(y);
 % figure;findpeaks(y * -1);
 [~,locs1] = findpeaks(y);
 [~,locs2] = findpeaks(y * -1);
 locs = sort(cat(1,locs1,locs2));
 locs(end 1) = find(y > initialState, 1, 'last');
 
 % wnIndex = locs;
 wnIndex = [];
 % humanReationTime = 0.25; % seconds
 hrtIndex = round(humanReactionTime * (length(y) / x(end)),0);
 for i = 1:length(locs)
     % Init of window
     if mod(i,2) ~= 0
         wnIndex(end 1) = locs(i) - hrtIndex;
     % End of the window
     else
         % Test
         wnIndex(end 1) = locs(i) - hrtIndex;
     end
 end
 wnIndex = unique(wnIndex);
 wnTime = x(wnIndex);
 
 if plotFlag
     figure;
     plot(x,y,'k-',x(wnIndex),y(wnIndex),'ro');
 end
 
 if length(refTime) ~= 0
 %     wnTime = wnTime;
     wnIndex2 = [];
     for i = 1:length(wnTime)
         idx = find(refTime > wnTime(i)-0.01 & ...
             refTime < wnTime(i) 0.01);
         wnIndex2(end 1) = idx(1);
     end
     wnIndex = unique(wnIndex2);
 %     wnTime = x(wnIndex);
 end
 
 end
 
diff --git a/GetDir.m b/GetDir.m
@@ -0,0  1,25 @@
 % Function name....: GetDir
 % Date.............: Jun 16, 2019
 % Mod date.........:
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                   Get dir elements (folder names or file names)
 % Parameters.......:
 %                   path
 %                   folderOrFile -> 1 for folders or 2 for files
 % Return...........:
 %                   elementNames -> list of all element names
 % Remarks..........:
 %
 function [elementNames] = GetDir(path, folderOrFile)
 files = dir(path);
 if folderOrFile == 1
     flags = [files.isdir];
 else
     flags = ~[files.isdir];
 end
 elements = files(flags);
 elementNames = {elements(:).name};
 
 end
 
diff --git a/IMNFWrap.m b/IMNFWrap.m
@@ -0,0  1,73 @@
 % Function name....: IMNFWrap
 % Date.............: Jan 28, 2019
 % Mod date.........: -
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                    Estimates Hilbert Spectrum to track instantaneous mean
 %                    frequency (IMNF) and calculates boxplot stats variables
 %                    (descriptive statistics) of IMNF.
 % Parameters.......:
 %                    y -> input time-series
 %                    fs -> sample rate
 %                    rmHighestFreq -> flag to remove highest
 %                    frequency of input time-series
 %                    plotFlag -> flag to indicate to plot figures
 %                    label ->
 % Return...........:
 %                    imnf_out2 -> estimated instantaneous mean frequency
 %                    ds -> key values of the descriptive statistics
 % Remarks..........:
 %                    Ref.: A. O. ANDRADE, P. Kyberd, and S. J. Nasuto,
 %                    “The application of the Hilbert spectrum to the analysis of electromyographic signals,”
 %                    Inf. Sci. (Ny)., vol. 178, no. 9, pp. 2176–2193, May 2008.
 %                    
 %                    Used code for descriptive statistics: https://www.mathworks.com/matlabcentral/fileexchange/29305-descriptive-statistics
 
 function [imnf_out2, ds] = IMNFWrap(y, fs, rmHighestFreq, plotFlag, label)
 t = [0:1:length(y)-1]/fs;
 % Remove highest frequency component before proceed
 if rmHighestFreq == 1
     disp('Removing highest frequency component..');
     [IMFs,residue]= sig_to_imf(y,1e-5,2);
     y = y - IMFs(1,:);
 end
 
 % Estimating intrinsic mode functions
 [IMFs,residue]= sig_to_imf(y,1e-5,2);
 LFreq=0;
 UFreq=fs/2;
 n_bins= 400;
 [m_a_p,minFreq,maxFreq, hs_dt] = ...
     plotHS1(IMFs,t,fs,n_bins,[LFreq UFreq],0);
 %Generating auxiliary time and frequency vectors
 f = 1:1:n_bins;
 f = convScale(1,n_bins,f,minFreq,maxFreq);
 %Performing energy normalization (between 0 and 1)
 [m_a_p]= convScale(min(min(m_a_p)),max(max(m_a_p)),m_a_p,0,1);
 
 T=(0:1:n_bins-1)*hs_dt;
 B=m_a_p;
 F = f;
 [imnf_out2] = IMNF(y,B,F);
 
 ds = getDescriptiveStatistics(imnf_out2);
 
 if plotFlag == 1
     figure;
     subplot(2,1,1);
     plot(t,y); xlabel('time (s)'); ylabel('Unit');
     title(label,'FontSize',20);
 
     subplot(2,1,2);
     imagesc(T,f,m_a_p); axis xy; colormap('hot');drawnow;
     hold on; plot(T,imnf_out2,'g'); ylabel('HS (Hz)'); xlabel('time (s)');
 
     figure;
     %     subplot(3,1,3);
     boxplot(imnf_out2,'orientation','vertical');
     ylim([0 fs/2]); ylabel('Instantaneous mean frequency (Hz)');
     title(label,'FontSize',20);
     set(gca,'XTickLabel','');
 end
 end
 
diff --git a/IMNFWrap1.m b/IMNFWrap1.m
@@ -0,0  1,90 @@
 % Function name....: IMNFWrap1
 % Date.............: Jan 28, 2019
 % Mod date.........: Aug 20, 2019
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                    Estimates Hilbert Spectrum to track instantaneous mean
 %                    frequency (IMNF) and calculates boxplot stats variables
 %                    (descriptive statistics) of IMNF.
 % Parameters.......:
 %                    y -> input time-series
 %                    fs -> sample rate
 %                    rmHighestFreq -> flag to remove highest
 %                    frequency of input time-series
 %                    plotFlag -> flag to indicate to plot figures
 %                    label ->
 % Return...........:
 %                    imnf_out2 -> estimated instantaneous mean frequency
 %                    ds -> key values of the descriptive statistics
 %                    hfRemoved -> estimated mean frequency of highest
 %                    frequency component, if rmHighestFreq = 1
 % Remarks..........:
 %                    Ref.: A. O. ANDRADE, P. Kyberd, and S. J. Nasuto,
 %                    “The application of the Hilbert spectrum to the analysis of electromyographic signals,”
 %                    Inf. Sci. (Ny)., vol. 178, no. 9, pp. 2176–2193, May 2008.
 %                    
 %                    Used code for descriptive statistics: https://www.mathworks.com/matlabcentral/fileexchange/29305-descriptive-statistics
 
 function [imnf_out2, ds, hfRemoved] = IMNFWrap1(y, fs, rmHighestFreq, plotFlag, label)
 t = [0:1:length(y)-1]/fs;
 % Remove highest frequency component before proceed
 if rmHighestFreq == 1
     disp('Removing highest frequency component..');
     [IMFs,residue]= sig_to_imf(y,1e-5,2);
     y = y - IMFs(1,:);
     % Estimate main frequency of highest frequency removed component
 %     [pxx,f] = pwelch(IMFs(1,:),60,10,60,fs);
 %     if plotFlag == 1
 %         figure; 
 %         plot(f,10*log10(pxx))
 %         xlabel('Frequency (Hz)')
 %         ylabel('PSD (dB/Hz)')
 %     end
 %     index = find(pxx == max(pxx));
 %     hfRemoved = f(index);
 
     hfRemoved = meanfreq(IMFs(1,:), fs);
 %     hfRemoved = medfreq(IMFs(1,:), fs);
 else
     hfRemoved = -1;
 end
 
 % Estimating intrinsic mode functions
 [IMFs,residue]= sig_to_imf(y,1e-5,2);
 LFreq=0;
 UFreq=fs/2;
 n_bins= 400;
 [m_a_p,minFreq,maxFreq, hs_dt] = ...
     plotHS1(IMFs,t,fs,n_bins,[LFreq UFreq],0);
 %Generating auxiliary time and frequency vectors
 f = 1:1:n_bins;
 f = convScale(1,n_bins,f,minFreq,maxFreq);
 %Performing energy normalization (between 0 and 1)
 [m_a_p]= convScale(min(min(m_a_p)),max(max(m_a_p)),m_a_p,0,1);
 
 T=(0:1:n_bins-1)*hs_dt;
 B=m_a_p;
 F = f;
 [imnf_out2] = IMNF(y,B,F);
 
 ds = getDescriptiveStatistics(imnf_out2);
 
 if plotFlag == 1
     figure;
     subplot(2,1,1);
     plot(t,y); xlabel('time (s)'); ylabel('Unit');
     title(label,'FontSize',20);
 
     subplot(2,1,2);
     imagesc(T,f,m_a_p); axis xy; colormap('hot');drawnow;
     hold on; plot(T,imnf_out2,'g'); ylabel('HS (Hz)'); xlabel('time (s)');
 
     figure;
     %     subplot(3,1,3);
     boxplot(imnf_out2,'orientation','vertical');
     ylim([0 fs/2]); ylabel('Instantaneous mean frequency (Hz)');
     title(label,'FontSize',20);
     set(gca,'XTickLabel','');
 end
 end
 
diff --git a/PeakFeatures.m b/PeakFeatures.m
@@ -0,0  1,38 @@
 % Function name....: PeakFeatures
 % Date.............: Feb 02, 2018
 % Author...........: Fabio Henrique, ([email protected])
 % Description......:
 %                    
 % Parameters.......:
 %                    pks -> 
 %                    locs -> in seconds
 %                    printFlag -> 
 % Return...........:
 %                    mpi -> median peak interval
 %                    mpa -> median peak amplitude
 %                    extremesDiff -> difference between last and first
 %                    peaks
 % Remarks..........:
 %                    
 function [nPeaks, mpi, mpa, extremesDiff] = PeakFeatures(pks, locs, printFlag)
 nPeaks = length(pks);
 
 % Calc median lag between peaks
 %tLocs = locs/fs;
 mpi = median(diff(locs));
 
 % Calc median amp peaks
 mpa = median(pks);
 
 % Calc time between last and first peak
 extremesDiff = locs(end) - locs(1);
 
 if printFlag
    fprintf('\nNumber of peaks = %d\n', nPeaks); 
    fprintf('Median peak interval = %.3f\n', mpi); 
    fprintf('Median peak amp = %.3f\n', mpa);
    fprintf('Last Peak - First Peak = %.3f\n', extremesDiff);
 end
 
 end