用希尔伯特黄变换(HHT)求时频谱和边际谱

1、Hilbert边际谱我觉得既然已经做出EMD了，也就是得到了IMF。

这个时候就是做hilbert幅值谱，然后对它积分就可以了。

程序不是很难搞到吧！我是用hspec画谱图的，自己又在后面添加了求边际谱的代码for k=1:size(E)bjp(k)=sum(E(k,:))*1/fs; %fs为采样频率;endfigureplot(bjp);xlabel('频率/ Hz');ylabel('幅值');比如我用两个正弦信号作仿真fs=1000;t=1/fs:1/fs:1;y1=2*sin(40*pi*t);y2=5*sin(80*pi*t);y=[y1,y2]; % 信号画出来的图很粗糙，更不用说对实际信号分析了，所以大家看看如何来修正？黄文章中边际谱对实际信号分析是很好的一条曲线我用hhspectrum算了一下谱图，同时求了一下边际谱，边际谱程序基本想法同form。

结果也不太好，20HZ处还行，40HZ就有些问题了，见附图你自己再用这个试试我没有用rilling的hhspectrumnspab:function h1= nspab(data,nyy,minw,maxw,dt)% The function NSPAB generates a smoothed HHT spectrum of data(n,k)% in time-frequency space, where% n specifies the length of time series, and% k is the number of IMF components.% The frequency-axis range is prefixed.% Negative frequency sign is reversed.%% MA TLAB Library function HILBERT is used to calculate the Hilbert transform. %% Example, [h,xs,w] = nspab(lod78_p',200,0,0.12,1,3224).%% Functions CONTOUR or IMG can be used to view the spectrum,% for example contour(xs,w,h) or img(xs,w,h).%% Calling sequence-% [h,xs,w] = nspab(data,nyy,minw,maxw,t0,t1)%% data - 2-D matrix data(n,k) of IMF components% nyy - the frequency resolution% minw - the minimum frequency% maxw - the maximum frequency% t0 - the start time% t1 - the end time% Output-% h - 2-D matrix of the HHT spectrum, where% the 1st dimension specifies the number of frequencies, % the 2nd dimension specifies the number of time values % xs - vector that specifies the time-axis values% w - vector that specifies the frequency-axis values% Z. Shen (JHU) July 2, 1995 Initial%----- Get dimensions (number of time points and components)[npt,knb] = size(data);%----- Get time interval%----- Apply Hilbert Transformdata=hilbert(data);a=abs(data);omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);%----- Smooth amplitude and frequencyfiltr=fir1(8,.1);for i=1:knba(:,i)=filtfilt(filtr,1,a(:,i));omg(:,i)=filtfilt(filtr,1,omg(:,i));end%----- Limit frequency and amplitudefor i=1:knbfor i1=1:npt-1if omg(i1,i) >=maxw,omg(i1,i)=maxw;a(i1,i)=0;elseif omg(i1,i)<=minw,omg(i1,i)=minw;a(i1,i)=0;endendendclear filtr data%va=var(omg(200:1200))%----- Get local frequencydw=maxw - minw;wmx=maxw;wmn=minw;%----- Construct the ploting matrixclear p;h1=zeros(npt-1,nyy+1);p=round(nyy*(omg-wmn)/dw)+1;for j1=1:npt-1for i1=1:knbii1=p(j1,i1);h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);endend%----- Do 3-point to 1-point averaging[nx,ny]=size(h1);%n1=fix(nx/3);%h=zeros(n1,ny);%for i1=1:n1%h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:)); %end%clear h1;%----- Do 3-points smoothing in x-directionfltr=1./3*ones(3,1);for j1=1:nyh1(:,j1)=filtfilt(fltr,1,h1(:,j1));endclear fltr;%----- Define the results%w=linspace(wmn,wmx,ny-1)';%xs=linspace(t0,t1,nx)';h1=flipud(rot90(h1));h1=h1(1:ny-1,:);form求边际谱时所用程序是没有问题的，用的是矩形积分公式。

hht边际谱计算
在光学中，Hartmann–Hartmann–Tscherning（HHT）边际谱计算是一种用于分析眼睛的光学性能的方法。

它可以用来评估眼睛的像差、波前畸变和其他光学特性。

HHT边际谱计算的过程通常包括以下步骤：
1. 采集数据：首先，需要采集眼睛通过适当的仪器（如自动视力检查仪或自适应光学系统）所产生的数据。

这可能包括从点光源发出的光线通过眼睛后在屏幕上形成的光斑的信息。

2. 分析数据：接下来，通过计算机程序对采集到的数据进行分析。

这可能涉及将数据转换成波前形状或像差的信息。

通常使用Zernike多项式来表示波前形状。

3. 边际谱计算：在分析数据的过程中，可以利用边际谱计算方法来评估眼睛的光学性能。

这个过程会衡量不同波长下的像差水平以及光学系统的性能。

总的来说，HHT边际谱计算是一个用于评估眼睛光学性能的重要工具，它可以帮助我们理解眼睛的视觉质量以及潜
在的视觉问题。

这对于眼科医生和光学工程师来说都非常有价值。

Python希尔伯特黄变换（Python Hilbert-Huang Transform，简称HHT）是一种复杂非线性信号分析方法，结合了希尔伯特变换和黄变换的优势，能够有效地对非线性和非平稳信号进行时频谱分析。

本文将从HHT的原理、基本步骤和Python实现方法三个方面进行介绍。

一、HHT的原理1.希尔伯特变换希尔伯特变换是一种将实数信号转换为解析信号的数学方法，通过对原信号进行傅立叶变换得到频谱信息，再对频谱信息进行一定的处理得到解析频谱，从而实现信号的解析表示。

希尔伯特变换的核心是求出原信号的解析函数，即原信号的复数形式，其中实部是原信号本身，虚部是原信号的希尔伯特变换。

希尔伯特变换在信号处理领域有着广泛的应用，能够提取信号的瞬时特征，对非平稳信号进行时频分析具有很高的效果。

2.黄变换黄变换是一种局部线性和非线性信号分解方法，可以将非线性和非平稳信号分解成若干个固有模态函数（Intrinsic Mode Function，简称IMF）的线性组合。

黄变换首先对原信号进行极值点的提取，然后通过极值点之间的插值得到包络线，再将原信号减去包络线得到一维信号，并对得到的一维信号进行数据挑选和插值，最终得到IMF。

多次重复以上步骤，直到原信号能够被分解为若干个IMF，再通过IMF的线性组合得到原信号的近似表示。

3.HHT的结合HHT将希尔伯特变换和黄变换结合在一起，利用希尔伯特变换提取信号的瞬时特征，再通过黄变换将信号分解成若干个IMF，从而能够更准确地描述信号的时频特性。

HHT的优势在于能够适用于非线性和非平稳信号，对信号的局部特征具有很好的描述能力，因此在振动信号分析、生物医学信号处理等领域有着广泛的应用。

二、HHT的基本步骤1.信号分解HHT首先对原信号进行希尔伯特变换，得到信号的瞬时频率特征，然后通过黄变换将信号分解成若干个IMF。

2.IMF的提取针对得到的IMF，需要对每个IMF进行较为严格的判别，确定其是否符合IMF的特征：极值点交替出现、包络线对称、局部频率单调。

希尔伯特-黄变换（Hilbert-Huang Transform，HHT）0 前言传统的数据分析方法都是基于线性和平稳信号的假设，然而对实际系统，无论是自然的还是人为建立的，数据最有可能是非线性、非平稳的。

希尔伯特-黄变换（Hilbert-Huang Transform，HHT）是一种经验数据分析方法，其扩展是自适应性的，所以它可以描述非线性、非平稳过程数据的物理意义。

1 HHT简介[贺礼平.希尔伯特-黄变换在电力谐波分析中的应用研究[D].湖南:中南大学,2009]HHT的发展。

1995年，Norden E.Huang为研究水表面波构思出一种所谓“EMD--HSA”的时间序列分析法，通过这种方法他发现水波的演化不是连续的，而是突变、离散、局部的。

1998年，Norden E.Huang等人提出了经验模态分解方法，并引入了Hilbert谱的概念和Hilbert谱分析的方法，美国国家航空和宇航局（NASA）将这一方法命名为Hilbert-Huang Transform，简称HHT，即希尔伯特-黄变换。

HHT是一种新的分析非线性非平稳信号的时频分析方法，由两部分组成：第一部分为经验模态分解（Empirical Mode Decomposition，EMD）（the sifting process，筛选过程），它是由Huang提出的,基于一个假设：任何复杂信号都可以分解为有限数目且具有一定物理定义的固有模态函数（Intrinsic Mode Function，IMF；也称作本征模态函数）；EMD方法能根据信号的特点，自适应地将信号分解成从高到低不同频率的一系列IMF；该方法直接从信号本身获取基函数，因此具有自适应性，同时也存在计算量大和模态混叠的缺点。

第二部分为Hilbert谱分析（Hilbert Spectrum Analysis，HSA），利用Hilbert变换求解每一阶IMF 的瞬时频率，从而得到信号的时频表示，即Hilbert谱。

频谱管理中基于希尔伯特 - 黄变换的干扰诊断张长根;王玉文;孟凡计【摘要】为了给频谱管理的干扰协调机制提供科学有效的依据,在分析系统间EMC 干扰的特点基础上,采用希尔伯特 - 黄变换分析接收机端的干扰信号样本的时 - 频信息、希尔伯特谱和边际谱,得到干扰信号的时频特征.采用多级频率搬移的方法,解决信号的经验模态分解中出现的模态混叠.文章以互调干扰信号的诊断分析为例,验证此法能够将互调干扰从信号样本中分离.【期刊名称】《软件》【年(卷),期】2011(032)006【总页数】4页(P70-72,77)【关键词】干扰诊断;多级频率搬移;希尔伯特;-;黄变换;模态混叠【作者】张长根;王玉文;孟凡计【作者单位】电子科技大学空天科学技术研究院,成都,610054;电子科技大学空天科学技术研究院,成都,610054;电子科技大学空天科学技术研究院,成都,610054【正文语种】中文【中图分类】TP391由于信息技术的不断发展，电磁环境越来越复杂，这使得电磁兼容性（Electromagnetic com patibility，即EMC）问题越来越突出。

为了保障正常的通信，有效地利用有限的频谱资源，建立有效的干扰协调机制显得尤为重要[1]。

为了给干扰协调提供科学、准确的技术依据，就需掌握干扰信号特征，包括基本参数和频谱特性参数。

信号的干扰是非常复杂的，可能是平稳信号或非平稳信号，频域上也可能非常相近。

传统的方法主要是从信号的频率、频率误差、射频电平、发射带宽和调制度等方面来测量信号的特征[1]。

这些方法大都只适用于平稳信号，并且对平稳信号的分析能力还是有限的，特别是频率相近的干扰[2]。

因此，本文从分析系统间电磁兼容性干扰场景出发，根据受干扰台站提供的干扰信号样本，采用希尔伯特-黄变换[3]（H ilbert-Huang Transform，HHT）对干扰信号进行细化分析，获得详细的干扰信号的时频特征和边际谱。

用希尔伯特黄变换（HHT）求时频谱和边际谱1.什么是HHT？HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。

2.EMD分解的步骤。

EMD分解的流程图如下：3.实例演示。

给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t)+5*sin(2*pi*35t)(1)为了对比，先用fft对求上述信号的幅频和相频曲线。

代码:function fftfenxiclear;clc;N=2048;%fft默认计算的信号是从0开始的t=linspace(1,2,N);deta=t(2)-t(1);1/detax=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);% N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;%x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*det a).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t);y = x;m=0:N-1;f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的%下面计算的Y就是x(t)的傅里叶变换数值%Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值Y=fft(y);z=sqrt(Y.*conj(Y));plot(f(1:100),z(1:100));title('幅频曲线')xiangwei=angle(Y);figure(2)plot(f,xiangwei)title('相频曲线')figure(3)plot(t,y,'r')%axis([-2,2,0,1.2])title('原始信号')（2）用Hilbert 变换直接求该信号的瞬时频率代码:clear;clc;clf;%假设待分析的函数是z=t^3N=2048;%fft 默认计算的信号是从0开始的t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);z=x;hx=hilbert(z);xr=real(hx);xi=imag(hx);%计算瞬时振幅sz=sqrt(xr.^2+xi.^2);%计算瞬时相位sx=angle(hx);%计算瞬时频率dt=diff(t);dx=diff(sx);sp=dx./dt;plot(t(1:N-1),sp)title('瞬时频率')小结：傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。

基于希尔伯特黄变换的刀具磨损特征提取
孙惠斌;牛伟龙;王俊阳
【期刊名称】《振动与冲击》
【年(卷),期】2015(000)004
【摘要】概述了希尔伯特黄变换（HHT）的基本理论和算法，对信号经过经验模态分解（EMD）后得到的固有模态函数（IMF）求取振幅均值，差值筛选出与刀具磨损相关的 IMF 分量，并对单分量固有模态函数求取边际谱，获取边际谱最大幅值点，建立他们与刀具磨损之间的映射关系，进行特征提取，将其作为神经网络的输入特征向量，结合希尔伯特三维时频谱进行刀具磨损状态的判断。

研究结果证明，该方法可以作为刀具磨损监测中信号特征提取的一种简单和可靠的方法。

【总页数】8页(P158-164,183)
【作者】孙惠斌;牛伟龙;王俊阳
【作者单位】西北工业大学机电学院，西安 710072;西北工业大学机电学院，西安 710072;西北工业大学机电学院，西安 710072
【正文语种】中文
【中图分类】TH165;TN911
【相关文献】
1.基于希尔伯特-黄变换和等距特征映射的刀具磨损状态监测 [J], 宋伟杰;关山;庞弘阳
2.钛合金车削加工过程中刀具磨损状态监测的小波包子带能量变换特征提取新方法
[J], 陈侃;傅攀;李威霖;曹伟青
3.基于Hough变换的刀具磨损监测加工表面纹理特征提取 [J], 郑建明;李鹏阳;李言;朱云飞
4.基于改进希尔伯特黄的故障特征提取方法研究 [J], 沈颉; 郭欣; 何嘉
5.基于改进希尔伯特黄的故障特征提取方法研究 [J], 沈颉;郭欣;何嘉
因版权原因，仅展示原文概要，查看原文内容请购买。