新旧版MATLAB中的希尔伯特-黄变换(HHT)及其边际谱的求取问题
2024-05-18 10:24:18
最近，由于科研工作需求，需要处理某一类波动信号，遂利用经验模态分解（EMD）、集合经验模态分解（EEMD）和希尔伯特-黄变换对上述信号进行处理，期望获取有用的本征模态函数（IMF）、频率信息和边际谱等。
在信号实施和处理时发现，网上大多代码均由较早时期的MATLAB完成，彼时还没有自带的HHT函数。而在MATLAB2018a及其以后的版本中，集成和自带了HHT函数。
请各位网友解答我的疑惑，不胜感激！！！
问题1：函数hhspectrum()返回的是否为归一化后的频率？为何？
函数instfreq中有如下注释
The result FNORMHAT lies between 0.0 and 0.5.
问题2：函数toimage()的原理是什么？为什么以各个IMF分量的瞬时幅值和瞬时频率为输入得到的im是N*N的矩阵（N表示信号长度）？
[im,tt,ff] = TOIMAGE(A,f,t,splx,sply) transforms a spectrum made of 1D functions (e.g., output of "hhspectrum") in an 2D image
问题3：老版MATLAB求取边际谱时为何乘以1fs\frac{1}{fs}fs1？
for k=1:size(E,1)bjp(k)=sum(E(k,:))*1/fs;
end
问题4：如何根据新版MATLAB求取边际谱？
[P,F,T,insf,inse]=hht(d(5:N_imf,:)',fs);
关于hht函数的解释，请参照MATLAB——Hilbert-Huang Transform （HHT）
问题5：如何根据新版MATLAB绘制三维时频谱图？
hht(d(5:N_imf,:)',fs);
上述代码可以直接绘制二维希尔伯特时频谱图，横轴表示时间、纵轴表示频率，颜色表示能量
那么三维的时频谱图该如何利用[P,F,T,insf,inse]=hht(d(5:N_imf,:)’,fs);的输出结果进行绘制？
下面首先展示老版MATLAB求取边际谱：
%% 获取源数据，signal1,signal2,signal3分别表示三种不同的信号类型
[date,time,signal1,signal2,signal3]=textread('BulgingData\aaa.txt','%s%s%f%f%f','headerlines',0);
timestring="05:25:00.00";
for i=1:length(date)if time(i)==timestringfor j=1:2048*4globalsignal(j)=signal1(i+j-1);endend
endglobalsignal=globalsignal';%% 定义采样频率、时间、采样点数
fs=25;
t = (0:length(globalsignal)-1)/fs;
N = length(globalsignal);%% 对globalsignal进行EEMD分解
d=eemd(globalsignal,3,100);
N_imf=fix(log2(length(globalsignal)));
d=d';%% 绘制分解结果
figure()
for i=1:length(d(:,1))subplot(length(d(:,1)),1,i)plot(t, d(i,:))
end%% 边际谱求取过程
c=d(2:N_imf,:);
[m,n]=size(c);
% 对IMF分量求取瞬时频率与振幅：
% A：每个IMF的振幅向量; fa:每个IMF对应的瞬时频率(归一化后); t：时间序列号
[A,fa,tt]=hhspectrum(c);
% 将每个IMF信号合成求取Hilbert谱
% E：对应的振幅值; Cenf：每个网格对应的中心频率
% 时频图(用颜色表示第三维值的大小); 三维图（三维坐标系：时间，中心频率，振幅）
[E,tt1]=toimage(A,fa,tt,length(tt)); % A,fa为固定参数, tt,length(tt)为可变参数figure()
% 边际谱是hilbert谱对时间的积分，从积分的角度来讲，
% 对任意一阶频率把所有的时间上的幅值都加起来，反映这阶频率在所有时间上的幅值积累
for k=1:size(E,1)bjp(k)=sum(E(k,:))*1/fs;
end
f=(1:N)/N*(fs/2);
plot(f(1:2000),bjp(1:2000));
xlabel('频率 / Hz');
ylabel('信号幅值');
ylim([0 80])
title('信号边际谱')%要求边际谱必须先对信号进行EMD分解
利用上述代码，可以求取信号的边际谱，代码中涉及的hhspectrum、toimage以及instfreq如下所示。
function [A,f,tt] = hhspectrum(x,t,l,aff)
%
% [A,f,tt] = HHSPECTRUM(x,t,l,aff) computes the Hilbert-Huang spectrum
%
% inputs:
% - x : matrix with one signal per row
% - t : time instants
% - l : estimation parameter for instfreq (integer >=1 (1:default))
% - aff : if 1, displays the computation evolution
%
% outputs:
% - A : instantaneous amplitudes
% - f : instantaneous frequencies
% - tt : truncated time instants
%
% calls:
% - hilbert : computes the analytic signal
% - instfreq : computes the instantaneous frequency
% - disprog : displays the computation evolution
%
%Examples:
%
%s = randn(1,512);
%imf = emd(s);
%[A,f,tt] = hhspectrum(imf(1:end-1,:));
%
%s = randn(10,512);
%[A,f,tt] = hhspectrum(s,1:512,2,1);
%
% rem: need the Time-Frequency Toolbox (http://tftb.nongnu.org)
%
% See also
% emd, toimage, disp_hhs
%
% G. Rilling, last modification 3.2007
% gabriel.rilling@ens-lyon.fr
error(nargchk(1,4,nargin));if nargin < 2t=1:size(x,2);endif nargin < 3l=1;endif nargin < 4aff = 0;endif min(size(x)) == 1if size(x,2) == 1x = x';if nargin < 2t = 1:size(x,2);endendNmodes = 1;
elseNmodes = size(x,1);
endlt=length(t);tt=t((l+1):(lt-l));for i=1:Nmodesan(i,:)=hilbert(x(i,:)')';f(i,:)=instfreq(an(i,:)',tt,l)';A=abs(an(:,l+1:end-l));if affdisprog(i,Nmodes,max(Nmodes,100))end
end
end
function [fnormhat,t]=instfreq(x,t,L,trace)
%INSTFREQ Instantaneous frequency estimation.
% [FNORMHAT,T]=INSTFREQ(X,T,L,TRACE) computes the instantaneous
% frequency of the analytic signal X at time instant(s) T, using the
% trapezoidal integration rule.
% The result FNORMHAT lies between 0.0 and 0.5.
%
% X : Analytic signal to be analyzed.
% T : Time instants (default : 2:length(X)-1).
% L : If L=1, computes the (normalized) instantaneous frequency
% of the signal X defined as angle(X(T+1)*conj(X(T-1)) ;
% if L>1, computes a Maximum Likelihood estimation of the
% instantaneous frequency of the deterministic part of the signal
% blurried in a white gaussian noise.
% L must be an integer (default : 1).
% TRACE : if nonzero, the progression of the algorithm is shown
% (default : 0).
% FNORMHAT : Output (normalized) instantaneous frequency.
% T : Time instants.
%
% Examples :
% x=fmsin(70,0.05,0.35,25); [instf,t]=instfreq(x); plot(t,instf)
% N=64; SNR=10.0; L=4; t=L+1:N-L; x=fmsin(N,0.05,0.35,40);
% sig=sigmerge(x,hilbert(randn(N,1)),SNR);
% plotifl(t,[instfreq(sig,t,L),instfreq(x,t)]); grid;
% title ('theoretical and estimated instantaneous frequencies');
%
% See also KAYTTH, SGRPDLAY.% F. Auger, March 1994, July 1995.
% Copyright (c) 1996 by CNRS (France).
%
% ------------------- CONFIDENTIAL PROGRAM --------------------
% This program can not be used without the authorization of its
% author(s). For any comment or bug report, please send e-mail to
% f.auger@ieee.orgif (nargin == 0)error('At least one parameter required');
end
[xrow xcol] = size(x)
if (xcol~=1)error('X must have only one column');
endif (nargin == 1)t=2:xrow-1; L=1; trace=0.0;
elseif (nargin == 2)L = 1; trace=0.0;
elseif (nargin == 3)trace=0.0;
endif L<1,error('L must be >=1');
end
[trow,tcol] = size(t);
if (trow~=1),error('T must have only one row');
end;if (L==1),if any(t==1)|any(t==xrow),error('T can not be equal to 1 neither to the last element of X');elsefnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi);end;
elseH=kaytth(L);if any(t<=L)|any(t+L>xrow),error('The relation L<T<=length(X)-L must be satisfied');elsefor icol=1:tcol,if trace, disprog(icol,tcol,10); end;ti = t(icol); tau = 0:L;R = x(ti+tau).*conj(x(ti-tau));M4 = R(2:L+1).*conj(R(1:L));diff=2e-6;tetapred = H * (unwrap(angle(-M4))+pi);while tetapred<0.0 , tetapred=tetapred+(2*pi); end;while tetapred>2*pi, tetapred=tetapred-(2*pi); end;iter = 1;while (diff > 1e-6)&(iter<50),M4bis=M4 .* exp(-j*2.0*tetapred);teta = H * (unwrap(angle(M4bis))+2.0*tetapred);while teta<0.0 , teta=(2*pi)+teta; end;while teta>2*pi, teta=teta-(2*pi); end;diff=abs(teta-tetapred);tetapred=teta; iter=iter+1;end;fnormhat(icol,1)=teta/(2*pi);end;end;
end
end
%TOIMAGE transforms a spectrum made of 1D functions in an 2D image
%
% [im,tt,ff] = TOIMAGE(A,f,t,splx,sply) transforms a spectrum made
% of 1D functions (e.g., output of "hhspectrum") in an 2D image
%
% inputs : - A : amplitudes of modes (1 mode per row of A)
% - f : instantaneous frequencies
% - t : time instants
% - splx : number of columns of the output im (time resolution).
% If different from length(t), works only for uniform
% sampling.
% - sply : number of rows of the output im (frequency resolution).
% outputs : - im : 2D image of the spectrum
% - tt : time instants in the image
% - ff : centers of the frequency bins
%
% Examples : [im,tt,ff] = toimage(A,f);[im,tt] = toimage(A,f,t);[im,tt,ff] = toimage(A,f,sply);
% [im,tt,ff] = toimage(A,f,splx,sply);[im,tt,ff] = toimage(A,f,t,splx,sply);
%
%
% See also
% emd, hhspectrum, disp_hhs
%
% G. Rilling, last modification 3.2007
% gabriel.rilling@ens-lyon.frfunction [im,tt,ff] = toimage(A,f,varargin)DEFSPL = 400;error(nargchk(2,5,nargin));switch nargincase 2t = 1:size(A,2);sply = DEFSPL;splx = length(t);case 3if isscalar(varargin{1})t = 1:size(A,2);splx = length(t);sply = varargin{1};elset = varargin{1};splx = length(t);sply = DEFSPL;endcase 4if isscalar(varargin{1})t = 1:size(A,2);sply = varargin{1};splx = varargin{2};elset = varargin{1};sply = varargin{2};splx = length(t);endcase 5t = varargin{1};splx = varargin{2};sply = varargin{3};
end
if isvector(A)A = A(:)';f = f(:)';
endif issparse(A) || ~isreal(A) || length(size(A)) > 2error('A argument must be a real matrix')
end
if issparse(f) || ~isreal(f) || length(size(f)) > 2error('f argument must be a real matrix')
end
if any(size(f)~=size(A))error('A and f matrices must have the same size')
end
if issparse(t) || ~isreal(t) || ~isvector(t) || length(t)~=size(A,2)error('t argument must be a vector and its length must be the number of columns in A and f inputs')
end
if ~isscalar(splx) || ~isreal(splx) || splx ~= floor(splx) || splx <= 0error('splx argument must be a positive integer')
end
if ~isscalar(sply) || ~isreal(sply) || sply ~= floor(sply) || sply <= 0error('splx argument must be a positive integer')
endif any(diff(diff(t))) && splx ~= length(t)warning('toimage:nonuniformtimeinsants','When splx differs from length(t), the function only works for equally spaced time instants. You may consider reformating your data (using e.g. interpolation) before using toimage.')
endf = min(f,0.5);
f = max(f,0);indf = round(2*f*(sply-1)+1);
indt = repmat(round(linspace(1,length(t),splx)),size(A,1),1);
im = accumarray([indf(:),indt(:)],A(:),[sply,splx]);indt = indt(1,:);
tt = t(indt);
ff = (0:sply-1)*0.5/sply+1/(4*sply);end
新旧版MATLAB中的希尔伯特-黄变换(HHT)及其边际谱的求取问题相关推荐
- MATLAB希尔伯特黄变换HHT
这两天在学习希尔伯特黄变换,也就是HHT,趁着学习的劲赶紧整理整理,用的是MATLAB进行编程,所用到的工具箱便是EMD工具箱,链接如下,请自行下载. 希尔伯特黄变换HHT_HHT-电信代码类资源-C ...
- matlab柯西主值积分,希尔伯特-黄变换基本概念
摘自<机械故障诊断理论与方法> 希尔伯特-黄变换(Hilbert-Huang Transform)是由N.E.Huang等人与1998年提出的一种非线性.非平稳信号的分析处理方法.这种方法 ...
- 信号处理：希尔伯特黄变换
目录: 目录: 前言 简介 基本原理 经验模态分解 希尔伯特变换 特点 (1)HHT能分析非线性非平稳信号. (2)HHT具有完全自适应性. (3)HHT不受Heisenberg测不准原理制约--适合 ...
- 希尔伯特黄变换(Hilbert-Huang)原理、HHT求时频谱、边际谱，及MATLAB（2018rb）实现
目录 1. 经验模态分解: 2. 希尔伯特变换: 3. 方法缺陷: 4. MATLAB(2018rb版本)实现和探讨 ##边际谱 [若觉文章质量良好且有用,请别忘了点赞收藏加关注,这将是我继续分享的动 ...
- 经验模式分解（EMD）及希尔伯特-黄变换（HHT）简介及matlab实现
本文介绍过程涉及到两个链接工具包,可以自己网上搜索下载,以下提供了网盘下载的地址,因为作者主要做语音方面工作,所以后面的说明主要以说话人识别为例.(链接:https://pan.baidu.com/s ...
- python新旧特性过渡_网站改版时的一种新旧版过渡方案
网站改版时,需要考虑一个周全的过渡方案,其中不容忽视的一点就是对旧版的处理问题.即使借助完美的数据迁移方案可以使新版从内容上完全取代旧版,但我们仍然不应该立即彻底废除掉旧版,因为: 1.网民有可能通过 ...
- 【信号处理】Matlab实现希尔伯特-黄变换
1 内容介绍 1998年,Norden E. Huang(黄锷:中国台湾海洋学家)等人提出了经验模态分解方法,并引入了Hilbert谱的概念和Hilbert谱分析的方法,美国国家航空和宇航局(NASA ...
- 天微TM1650数码管驱动IC新旧版 驱动和注意事项
天微TM1650数码管驱动IC新旧版 驱动和注意事项 项目场景: 项目需要一个控制板和显示,通过一条1米数据线连接主控制 TM1650市面上多,价格便宜,使用简单, 相对于用逻辑门或单片机做,开发简单 ...
- css消除全部css_消除旧版浏览器中CSS3头痛
css消除全部css 你已经看到了它的时间 和 时间 再次上Webdesigntuts +; 您可能会想尝试该CSS3教程,但是由于旧版浏览器缺乏支持,因此您无法进一步研究. 但是,有很多工具可以在这 ...
- 量化择时：基于经验模态分解的希尔伯特-黄变换（二）算法
量化择时:基于经验模态分解的希尔伯特-黄变换 part2部分是算法的介绍,抛开代码部分,其实就是所有人都能看得懂字面解释 Part2算法 在了解了基础的数理知识和学习了将实信号转换为复信号的处理方法之 ...
最新文章
- 标准W3C盒子模型和IE盒子模型CSS布局经典盒子模型(转)
- javaweb回顾第二篇tomcat和web程序部署
- C语言之详解#ifdef等宏
- 怎么查这个文件在linux下的哪个目录
- Yii1.1 CGridView 简单使用
- cmos存储器中存放了_天津大学姚建铨院士，张雅婷副教授JMCC：具有宽光谱调控特性的阻变存储器...
- Python | Tkinter中的文本区域和按钮
- SAP License：生产订单结算时候的几个差异
- codesys编程_明晚20:00，CODESYS教您制作可编程控制器
- java是解释执行么
- 由有向图的邻接矩阵生成其可达矩阵
- 一阶广义差分模型_计量经济学习题第5章 自相关性
- 商丘服务器维修,商丘联想服务器维修网点
- w8ndows 秒表,关闭 Windows Search，Win8 能变快？
- 山重水复疑无路，柳暗花明又一村。---找工作感想
- c语言实现7段数码管显示,FPGA入门--七段数码管显示
- layer.photos 查看本地图片，并实现缩放和旋转功能
- 我支持平板能代替笔记本电脑
- 《偏生要鲜花着景，应这急景流年》
- 2021漳州一中历年高考成绩查询,2021年漳州中考录取分数线,历年漳州各高中录取分数线排名...