希尔伯特变换求瞬时频率的matlab动态实现

通过希尔伯特变换可以求得一个信号的包络曲线，瞬时频率，瞬时相位，瞬时谱等等。

计算方法在此不赘述，本文只演示希尔伯特变换求包络线和瞬时频率的效果。

代码：

clc;close all;clear;
x= 0:pi/500:5*pi;
y2=sin(x);
a=30;t=0;
figure(1);
y1=1+0.8*sin(0.3*2*pi*x); subplot(3,1,1);
plot(x,y1,'r--','Linewidth',1.2);
axis([0,15,-2,2]);
hold on;
P1=plot(x,y1,'k')
legend('模拟信号','载波信号')
title('模拟信号和载波信号');subplot(3,1,2);
P2=plot(x,y1,'r','Linewidth',1.2);
axis([0,15,-2,2]);
legend('包络线');
title('希尔伯特变换解包络');subplot(3,1,3);
P3=stem(a,'k')
axis([0,40,-1,40]);
legend(['为了便于观察，f2比f1提前一个单位']);
title('实际的瞬时频率f1与希尔伯特变换得到的瞬时频率f2');while a<40a=a+0.01;y2=sin(a*x);y3=y2.*y1;z = hilbert(y3');                               %运用希尔伯特变换生成解析信号inst_amplitude = abs(z)';                       %提取包络信号（瞬时幅值）inst_phase = unwrap(angle(z));                  %计算瞬时相位inst_frequency = diff(inst_phase)/pi*500;    %计算瞬时频率subplot(3,1,1);set(P1,'XData',x,'YData',y3);drawnow;subplot(3,1,2);set(P2,'XData',x,'YData',inst_amplitude);drawnow;%pause(0.01);t=t+0.01;subplot(3,1,3);set(P3,'XData',[t,t+1],'YData',[a,inst_frequency(1250)],'color','red');drawnow;
end

演示效果：

一开始载波频率比较低时，可以看出希尔伯特变换求得的包络线并不准确，求得的瞬时频率也并不完全与实际的频率相符合。

随着频率逐渐增大：

可以看出此时包络曲线与实际的调幅信号已经基本符合，瞬时频率也与实际的瞬时频率基本相等。

当频率较大时：

可以看出此时包络线已经较完美的与调幅信号相匹配，瞬时频率也与实际的瞬时频率相等。

但包络线两端还是有失真的现象。

代码中的instfreq函数需要自己编写：

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;
end ;if  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 ;

实际的代码是动态演示的，看的非常清楚，我不知道怎么发视频，就不发了，感兴趣的话可以自己跑一下代码。

希尔伯特变换求瞬时频率的matlab动态实现相关推荐

离散信号的希尔伯特变换的计算公式_希尔伯特变换和瞬时频率问题--连载（二）...
写在开始的一段话: PS:OK,上一期关于希尔伯特变换的文章发出后,有知友在评论区说"看到最后--居然这--",哈哈,其实我也挺愧疚大家的,明明一篇知识分享的文章,却写到结尾都没进 ...
瞬时频率函数matlab,Hilbert 变换与瞬时频率
Hilbert 变换与瞬时频率 Hilbert 变换仅可估计单分量信号的瞬时频率.单分量信号在时频平面中用单一"脊"来描述.单分量信号包括单一正弦波信号和 chirp 等信号. 生 ...
VMD分解，matlab代码，包络线，包络谱，中心频率，峭度值，能量熵，近似熵，包络熵，频谱图，希尔伯特变换，包含所有程序MATLAB代码，-西储大学数据集为例
目录 1.选取数据 2.VMD函数-matlab代码 3.采用matlab脚本导入数据并做VMD分解 4. VMD分解图 5.计算中心频率 6.画包络线 7. 画包络谱 8. 计算峭度值 9.计算能量 ...
Hilbert变换求信号的包络线及MATLAB代码
一.Hilbert变换表达式 Hilbert变换是信号与的卷积,表达式如下: 由于本质是卷积,因此可以从"线性系统","调幅-调频"等角度思考.进一步可参考如下 ...
【20220207】【信号处理】希尔伯特变换定义及解调原理
一.解析信号 1. 定义解析信号是没有负频率分量的复值函数,解析信号的实部和虚部是由希尔伯特变换相关联的实值函数. (参考:解析信号) 2. 概念一个实值函数的 Hilbert 变换记作为 , ...
[HT/NHT/DQ]-三种基于EMD的瞬时频率计算方法的比较
文章目录 0 前言 1 瞬时频率和经验模态分解 1.1 瞬时频率的定义 1.2 经验模态分解 2 瞬时频率的计算方法 2.1 HT方法 2.2 NHT方法 2.3 DQ方法 3 三种瞬时频率计算方法的 ...
包络谱分析和希尔伯特变换（Hilbert transform）
参考: 希尔伯特变换将信号表示为复解析信号的物理意义是什么? 希尔伯特变换和瞬时频率问题–连载(二) 希尔伯特变换求包络原理题主本硕机械专业,自学转互联网算法岗成功,获得阿里.字节.美团.华为等 ...
初识希尔伯特变换(Hilbert Transform)
要想很好的理解,需要有信号处理的前导知识希尔伯特 David Hilbert(1862-1943),德国数学家.人类智慧丰碑的一个家伙. 希尔伯特变换的物理意义是一种积分变换,把信号的所有频率分量 ...
希尔伯特-黄变换（HHT）的前世今生——一个从瞬时频率讲起的故事
一.来自小X的疑问从前有一个国家,叫做实国(一维国度),里边有个叫小X的小人儿. 小X是一根线段,他每天最爱做的事情就是跳舞. 因为小X的舞姿十分稳定,同伴们都说他的头部跳出了频率为1hz的正弦曲线 ...

希尔伯特变换求瞬时频率的matlab动态实现

希尔伯特变换求瞬时频率的matlab动态实现相关推荐

最新文章

热门文章