希尔伯特变换在MATLAB中的应用 您所在的位置:网站首页 希尔伯特谱分析的主要功能 希尔伯特变换在MATLAB中的应用

希尔伯特变换在MATLAB中的应用

2024-03-21 11:10| 来源: 网络整理| 查看: 265

原文出处:http://www.cnblogs.com/xingshansi/articles/6498913.html

一、基本理论

  A-Hilbert变换定义

对于一个实信号x(t)x(t),其希尔伯特变换为:

x~(t)=x(t)∗1πtx~(t)=x(t)∗1πt

式中*表示卷积运算。

Hilbert本质上也是转向器,对应频域变换为:

1πt⇔j⋅sign(ω)1πt⇔j⋅sign(ω)

即余弦信号的Hilbert变换时正弦信号,又有:

1πt∗1πt⇔j⋅sign(ω)⋅j⋅sign(ω)=−11πt∗1πt⇔j⋅sign(ω)⋅j⋅sign(ω)=−1

即信号两次Hilbert变换后是其自身相反数,因此正弦信号的Hilbert是负的余弦。

对应解析信号为:

z(t)=x(t)+jx~(t)z(t)=x(t)+jx~(t)

此操作实现了信号由双边谱到单边谱的转化。

  B-Hilbert解调原理

设有窄带信号:

x(t)=a(t)cos[2πfst+φ(t)]x(t)=a(t)cos⁡[2πfst+φ(t)]

其中fsfs是载波频率,a(t)a(t)是x(t)x(t)的包络,φ(t)φ(t)是x(t)x(t)的相位调制信号。由于x(t)x(t)是窄带信号,因此a(t)a(t)也是窄带信号,可设为:

a(t)=[1+∑m=1MXmcos(2πfmt+γm)]a(t)=[1+∑m=1MXmcos⁡(2πfmt+γm)]

式中,fmfm为调幅信号a(t)a(t)的频率分量,γmγm为fmfm的各初相角。

对x(t)x(t)进行Hilbert变换,并求解解析信号,得到:

z(t)=ej[2πfs+φ(t)][1+∑m=1MXmcos(2πfmt+γm)]z(t)=ej[2πfs+φ(t)][1+∑m=1MXmcos⁡(2πfmt+γm)]

A(t)=[1+∑m=1MXmcos(2πfmt+γm)]A(t)=[1+∑m=1MXmcos⁡(2πfmt+γm)]

Φ(t)=2πfst+φ(t)Φ(t)=2πfst+φ(t)

则解析信号可以重新表达为:

z(t)=A(t)ejΦ(t)z(t)=A(t)ejΦ(t)

对比x(t)x(t)表达式,容易发现:

a(t)=A(t)=x2(t)+x~2(t)−−−−−−−−−−√a(t)=A(t)=x2(t)+x~2(t)

φ(t)=Φ(t)−2πfst=arctanx(t)x~(t)−2πfstφ(t)=Φ(t)−2πfst=arctan⁡x(t)x~(t)−2πfst

由此可以得出:对于窄带信号x(t)x(t),利用Hilbert可以求解解析信号,从而得到信号的幅值解调a(t)a(t)和相位解调φ(t)φ(t),并可以利用相位解调求解频率解调f(t)f(t)。因为:

f(t)=12πdφ(t)dt=12πdΦ(t)dt−fsf(t)=12πdφ(t)dt=12πdΦ(t)dt−fs

  C-相关MATLAB指令

hilbert

功能:将实数信号x(n)进行Hilbert变换,并得到解析信号z(n).

调用格式:z = hilbert(x)

instfreq

功能:计算复信号的瞬时频率。

调用格式:[f, t] = insfreq(x,t)

示例:

1 2 z = hilbert(x); f = instfreq(z);

 

二、应用实例

 例1:给定一正弦信号,画出其Hilbert信号,直接给代码:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 clc clear  all close  all ts = 0.001; fs = 1/ts; N = 200; f = 50; k = 0:N-1; t = k*ts; % 信号变换 % 结论:sin信号Hilbert变换后为cos信号 y =  sin (2* pi *f*t); yh = hilbert(y);     % matlab函数得到信号是合成的复信号 yi =  imag (yh);       % 虚部为书上定义的Hilbert变换 figure subplot (211) plot (t, y) title ( '原始sin信号' ) subplot (212) plot (t, yi) title ( 'Hilbert变换信号' ) ylim ([-1,1])

  对应效果图:

例2:已知信号x(t)=(1+0.5cos(2π5t))cos(2π50t+0.5sin(2π10t))x(t)=(1+0.5cos⁡(2π5t))cos⁡(2π50t+0.5sin⁡(2π10t)),求解该信号的包络和瞬时频率。

分析:根据解包络原理知:

信号包络:(1+0.5cos(2π5t))(1+0.5cos⁡(2π5t))

瞬时频率:2π50t+0.5sin(2π10t)2π2π50t+0.5sin⁡(2π10t)2π

那么问题来了,实际情况是:我们只知道x(t)x(t)的结果,而不知道其具体表达形式,这个时候,上文的推导就起了作用:可以借助信号的Hilbert变换,从而求解信号的包络和瞬时频率。

对应代码:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 clear  all ;  clc ;  close  all ;   fs=400;                                  % 采样频率 N=400;                                   % 数据长度 n=0:1:N-1; dt=1/fs; t=n*dt;                                  % 时间序列 A=0.5;                                   % 相位调制幅值 x=(1+0.5* cos (2* pi *5*t)).* cos (2* pi *50*t+A* sin (2* pi *10*t));   % 信号序列 z=hilbert(x');                           % 希尔伯特变换 a= abs (z);                                % 包络线 fnor=instfreq(z);                        % 瞬时频率 fnor=[fnor(1); fnor; fnor( end )];         % 瞬时频率补齐 % 作图 pos =  get ( gcf , 'Position' ); set ( gcf , 'Position' ,[pos(1), pos(2)-100,pos(3),pos(4)]); subplot  211;  plot (t,x, 'k' );  hold  on; plot (t,a, 'r--' , 'linewidth' ,2); title ( '包络线' );  ylabel ( '幅值' );  xlabel ([ '时间/s'  10  '(a)' ]); ylim ([-2,2]); subplot  212;  plot (t,fnor*fs, 'k' );  ylim ([43 57]); title ( '瞬时频率' );  ylabel ( '频率/Hz' );   xlabel ([ '时间/s'  10  '(b)' ]);

  其中instfreq为时频工具包的代码,可能有的朋友没有该代码,这里给出其程序:

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 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 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 %   [email protected]   if  ( nargin  == 0),   error ( 'At least one parameter required' ); end ; [xrow,xcol] =  size (x); if  (xcol~=1),   error ( 'X must have only one column' ); end   if  ( 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 1e-6)&(iter


【本文地址】

公司简介

联系我们

今日新闻

    推荐新闻

    专题文章
      CopyRight 2018-2019 实验室设备网 版权所有