如何得到波形包络线及瞬时频率

如何得到波形包络线及瞬时频率

引言

许多领域中都用到信号包络线的提取,又称为包络检测或者幅度解调。提取的算法也比较多,本文主要采用希尔伯特变化实现。

原理

详细原理请见维基百科

例子

  • hilbert变换matlab源码
    function x = hilbert(xr,n)
    %HILBERT  Discrete-time analytic signal via Hilbert transform.
    %   X = HILBERT(Xr) computes the so-called discrete-time analytic signal
    %   X = Xr + i*Xi such that Xi is the Hilbert transform of real vector Xr.
    %   If the input Xr is complex, then only the real part is used: Xr=real(Xr).
    %   If Xr is a matrix, then HILBERT operates along the columns of Xr.
    %
    %   HILBERT(Xr,N) computes the N-point Hilbert transform.  Xr is padded with 
    %   zeros if it has less than N points, and truncated if it has more.  
    %
    %   For a discrete-time analytic signal X, the last half of fft(X) is zero, 
    %   and the first (DC) and center (Nyquist) elements of fft(X) are purely real.
    %
    %   EXAMPLE:
    %          Xr = [1 2 3 4];
    %          X = hilbert(Xr)
    %          % produces X=[1+1i 2-1i 3-1i 4+1i] such that Xi=imag(X)=[1 -1 -1 1] is the
    %          % Hilbert transform of Xr, and Xr=real(X)=[1 2 3 4].  Note that the last half
    %          % of fft(X)=[10 -4+4i -2 0] is zero (in this example, the last half is just
    %          % the last element).  Also note that the DC and Nyquist elements of fft(X)
    %          % (10 and -2) are purely real.
    %
    %   See also FFT, IFFT.
    %   Copyright 1988-2008 The MathWorks, Inc.
    %   $Revision: 1.10.4.5 $  $Date: 2010/12/06 00:02:22 $
    %   References:
    %     [1] Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time
    %     Signal Processing, 2nd ed., Prentice-Hall, Upper Saddle River, 
    %     New Jersey, 1998.
    %
    %     [2] S. Lawrence Marple, Jr., Computing the discrete-time analytic 
    %     signal via FFT, IEEE Transactions on Signal Processing, Vol. 47, 
    %     No. 9, September 1999, pp.2600--2603.
    if nargin<2, n=[]; end
    if ~isreal(xr)
    warning(message('signal:hilbert:Ignore'))
    xr = real(xr);
    end
    % Work along the first nonsingleton dimension
    [xr,nshifts] = shiftdim(xr);
    if isempty(n)
    n = size(xr,1);
    end
    x = fft(xr,n,1); % n-point FFT over columns.
    h  = zeros(n,~isempty(x)); % nx1 for nonempty. 0x0 for empty.
    if n > 0 && 2*fix(n/2) == n
    % even and nonempty
    h([1 n/2+1]) = 1;
    h(2:n/2) = 2;
    elseif n>0
    % odd and nonempty
    h(1) = 1;
    h(2:(n+1)/2) = 2;
    end
    x = ifft(x.*h(:,ones(1,size(x,2))));
    % Convert back to the original shape.
    x = shiftdim(x,-nshifts);
  • demo1 模拟信号
    clc;clear all;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));  % 信号序列
    subplot 211
    plot(t,x,'k')
    hold on
    ht = hilbert(x');                          % 希尔伯特变换
    Rht = abs(ht);                               % 包络线
    plot(t,Rht,'r--','linewidth',2);
    xlabel('time/s')
    ylabel('Amp')
    ylim([-1.8 1.8])
    subplot 212
    fnor=instfreq(ht);                       % 瞬时频率
    fnor=[fnor(1); fnor; fnor(end)];        % 瞬时频率补齐
    plot(t,fnor*fs,'k'); 
    xlabel('time/s')
    ylabel('Hz')

    file
    *瞬时频率提取函数matlab源码

    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.org
    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<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');
    else
    fnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi);
    end;
    else
    H=kaytth(L);
    if any(t<=L)|any(t+L>xrow),
    error('The relation L<T<=length(X)-L must be satisfied');
    else
    for 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;