希尔伯特变换在MATLAB中的应用 | 您所在的位置:网站首页 › 希尔伯特谱分析的主要功能 › 希尔伯特变换在MATLAB中的应用 |
原文出处: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=arctanx(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 实验室设备网 版权所有 |