折腾来折腾去

pikipity的blog

尝试两种利用 EMD 来计算瞬时频率的方法

这两天利用 Matlab 尝试了两种利用 EMD 来计算瞬时频率的方法。EMD 的介绍请看这里。我利用的是 G. Rilling 的程序计算的 IMF,程序在这里下载。尽管这个程序存在着这样那样的问题,但是他利用自己的假设找到了方法来解决,我曾经尝试着利用我在这篇文章中提到的看法来解决那些问题,但是发现,速度大大降低,结果改善却不明显,所以放弃而在这里使用他的程序。

这两种方法中,一种方法先使用 Hilbert Transform 求得 IMF 的瞬时相位,然后利用瞬时相位求得瞬时频率,另一种则不适用 Hilbert Transform,直接根据 IMF 类似弦波这一特性,利用局部极值点与零点来计算瞬时频率。

方法一:瞬时相位求瞬时频率

关于 Hilber Transform,我在这篇文章中已经介绍过了,它可以将信号变换为其对应的虚数形式:

$$z(t)=s(t)+j\hat{s}(t)$$

其中 $$s(t)$$ 是待分析的 IMF 分量,$$\hat{s}(t)$$ 是 Hilbert Transform 的结果,利用 Matlab 函数 z=hilbert(s) 可以直接得到 $$z(t)$$(注意,得到的是 $$z(t)$$,不是 $$\hat{s}(t)$$)。由虚数形式我们就可以很轻松地得到瞬时相位 $$\theta(t)$$ 了。然后利用下式就可以求得瞬时频率了。

$$2\pi f(t)=\omega(t)=\frac{d\theta(t)}{dt}$$

但是要注意的是,我们只能求得取样点的相位,也就是说,我们求得的相位是离散的,无法求导,这里我使用减法来求近似导数,例如我们要求第 $$n$$ 个点的瞬时频率,那么就利用下式求得

$$2\pi f(n)=\frac{\theta(n+1)-\theta(n-1)}{2\frac{1}{f_s}}$$

这里的 $$f_s$$ 是采样频率。对于数据的第一个点,就将 $$n-1$$ 改为 $$n$$,下面除以一个 $$\frac{1}{f_s}$$;对于数据的最后一个点,就将 $$n+1$$ 换为 $$n$$,下面除以一个 $$\frac{1}{f_s}$$。Matlab 函数如下,可以在这里下载。

function [ XPhase,f,Imag,t,f_OGZ,maxf_emdvalue ]=IF_hilbert(x,fs)
%x is the signal that will be analysised
%fs is sample frequency
%XPhase is instanepus phase
%f is instaneous phase
%t is the time string
%f_OGZ is the mean frequency
%maxf_emdvalue is the frequency whose property is maximum.

if length(x)<2
    error('The length of x is too short')
end
Ts=1/fs;
x=hilbert(x);
XPhase=phase(x);
for n=1:length(x)
     t(n)=n*Ts;
    Imag(n)=abs(x(n));
    if n==1
        f(n)=1/2/pi*(XPhase(2)-XPhase(1))/Ts;
    elseif n==length(x)
        f(n)=1/2/pi*(XPhase(n)-XPhase(n-1))/Ts;
    else
        f(n)=1/2/pi*(XPhase(n+1)-XPhase(n-1))/2/Ts;
    end
end
f_OGZ=sum(f)/length(f);
f_range=0:1:100;
f_pro=hist(f,f_range);
maxf_emd_pro=max(f_pro);
maxf_emd_pro_index=find(f_pro==maxf_emd_pro);
maxf_emdvalue=f_range(maxf_emd_pro_index);
end

为了方便以后的对比,我又求了瞬时频率的平均频率f_OGZ和瞬时频率出现概率最高的频率值maxf_emdvalue

方法二:利用极值点与零点来估算瞬时频率

根据 IMF 的特性,IMF都可以看做多个弦波的叠加,那么对于弦波,一个极值点与零点之间的时间就可以看做四分之一个波长,我们可以根据这个来估算频率,具体推导和理论请看这篇论文。重点就是利用下面这个公式

$$f=\frac{1}{12}\left{\frac{1}{T_4}+\left(\frac{1}{T_21}+\frac{1}{T_22}\right)+\left(\frac{1}{T_11}+\frac{1}{T_12}+\frac{1}{T_13}+\frac{1}{T_14}\right)\right}$$

或者

$$f=\frac{1}{7}\left{\frac{1}{4T_4}+\left(\frac{1}{2T_21}+\frac{1}{2T_22}\right)+\left(\frac{1}{T_11}+\frac{1}{T_12}+\frac{1}{T_13}+\frac{1}{T_14}\right)\right}$$

符号代表的意思见下图

上式符号代表的意思

很容易发现,这个式子中很注重局部极值点与零点的位置,所以继续使用 G. Rilling 的近似会带来很大的误差,所以,我重新求取了极值点与零点。对于零点,如果利用 G. Rilling 求取的零点位置上信号值不是零,那么在其附近找到一个点,两者的信号值乘积为负数,那么零点必然就在两者中间,作直线,求直线与横坐标的交点作为零点。对于极值点,我在 G. Rilling 的函数找到的极值点的附近再取六个点,利用这七个点的坐标求取一个二次方程 $$y=ax2+bx+c$$ 的系数,这个二次方程的极值点作为信号的极值点。Matlab 函数如下,下载在这里。其中求新的极值点与零点的函数[ newindmin newindmax newindzer ] = refresh_index( x,indmin,indmax,indzer )这里。我同样求了一下瞬时频率的平均频率 f_GZC

function [f,t,f_GZC]=Instaneous_Frequency_Test(x,fs)
% [f,t]=Instaneous_Frequency_Test(x,fs)[f,t]=Instaneous_Frequency_Test(x,fs)
% x is the signal that will be analysised.
% fs is the sampling frequency
% f is the instaneous frquency
% t is the time range of each instaneous frequency
% f_GZC is the mean frequency

%the length of x must be larger than 8
if length(x)<8
    error('the length of x must be larger than 7')
end
%clear results
f=[]; 
t=[];
%calculate sampling periode
Ts=1/fs;
%calculate local extrema and zero
[indmin indmax indzer]=extr(x);
%refresh local extrema and zero
[ newindmin newindmax newindzer ] = refresh_index( x,indmin,indmax,indzer );
indzer=newindzer;
indmin=newindmin;
indmax=newindmax;
%get the index of all indmin, indmax and indzer
%[index,I] = rerange_index( indmin,indmax,indzer );
index=sort([indmin indmax indzer]);
if length(index)<5
    error('Too few extrema and zero')
end
%prepare Begin point index and finish point index
L_index=length(index);
%begin calculate instaneous frquency
Begin_point=4;
Finish_point=L_index-4;
for Number_index=1:L_index-1
    if Number_index==1
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/3*(1/t_q/4+1/t_h_2/2+1/t_f_4)];
    elseif Number_index==2
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/5*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_3+1/t_f_4)];
    elseif Number_index==3
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/6*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_2+1/t_f_3+1/t_f_4)];
    elseif Number_index==L_index-1
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        f=[f 1/3*(1/t_q/4+1/t_h_1/2+1/t_f_1)];
    elseif Number_index==L_index-2
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        f=[f 1/5*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2)];
    elseif Number_index==L_index-3
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        f=[f 1/6*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2+1/t_f_3)];
    else
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/7*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2+1/t_f_3+1/t_f_4)];
    end
end
t=t.*1/fs;
f_GZC=sum(f)/length(f);
end

两种方法连同 fft 的比较

下面来比较一下这两种方法,因为 现在使用最多的还是直接用 fft 来计算频率值,所以这里连同 fft 一起进行比较比较。分析的数据,我是用的是实验室曾经记录的 SSVEP 的信号,数据可以在这里下载,将信号保存到变量x中运行下面这段程序(直接加载入 Matlab 就已经保存在变量x中了),程序在这里下载

% signal -> x
%begin
fs=600;
Ts=1/fs;
t=1:length(x);
t=(t-1).*Ts;
%figure('name','original signal');plot(t,x)
%IMF=eemd_test(x); %try to use eemd to calculate IOF
IMF=emd(x);
% original signal
figure('name','original signal')
plot(t,x)
xlabel('t (s)')
ylabel('Amplitude')
%IMF components
figure('name','IMF')
for n=1:size(IMF,1)
    subplot(size(IMF,1),1,n)
    plot(t,IMF(n,:))
    xlabel('t (s)')
    ylabel('Amplitude')
    %hold on
    %plot(t,IMF(n,:),'rx')
end
%fft for each IMF components
figure('name','fft')
maxf_fft=[];
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    [frequency,fft_result]=fft_plot(IMF(n,:),fs);
    plot(frequency,abs(fft_result(1:length(frequency))))
    xlabel('f (Hz)')
    ylabel('Amplitude')
    axis([0 100 0 max(abs(fft_result(1:length(frequency))))])
    maxf_fft_index=find(abs(fft_result(1:length(frequency)))==max(abs(fft_result(1:length(frequency)))));
    maxf_fftvalue=frequency(maxf_fft_index);
    maxf_fft(n-1)=maxf_fftvalue(1);
end
%IF not using hilbert transform
%IF property
figure('name','IF')
xlabel('frequency')
ylabel('probility (%)')
meanf=[];
maxf_emd=[];
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    [f,tf,f_GZC]=Instaneous_Frequency_Test(IMF(n,:),600);
    meanf(n-1)=f_GZC;
    f_range=0:1:100;%frequency range
    f_pro=hist(f,f_range);
    plot(f_range,f_pro./length(f)*100);
    xlabel('frequency')
    ylabel(strcat('IMF',num2str(n),'(%)'))
    maxf_emd_pro=max(f_pro);
    maxf_emd_pro_index=find(f_pro==maxf_emd_pro);
    maxf_emdvalue=f_range(maxf_emd_pro_index);
    maxf_emd(n-1)=maxf_emdvalue(1);
end
%IF in time domain
figure('name','IF value')
for h=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,h-1)
    xlabel('t (s)')
    ylabel('f (Hz)')
    [f,tf,f_GZC]=Instaneous_Frequency_Test(IMF(h,:),600);
    hold on
    for n=1:length(f)
        plot([tf(n,1),tf(n,2)],[f(n),f(n)],'k')
        if n~=length(f)
            plot([tf(n,2),tf(n,2)],[f(n+1),f(n)],'k')
        end
    end
end
%IF using hilbert transform
%IF property
figure('name','IF property hilbert')
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    xlabel('f (Hz)')
    ylabel('Number')
     [ XPhase,f,Imag,t,f_OGZ,maxf_emdvalue ]=IF_hilbert(IMF(n,:),fs);
     maxf_emd2(n-1)=maxf_emdvalue;
     f_range=0:1:100;%frequency range
    f_pro=hist(f,f_range);
    plot(f_range,f_pro./length(f)*100);
    xlabel('frequency')
    ylabel(strcat('IMF',num2str(n),'(%)'))
end
%IF in time domain
figure('name','IF value hilbert')
for h=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,h-1)
    xlabel('t (s)')
    ylabel('f (Hz)')
    [ XPhase,f,Imag,t,f_OGZ,maxf_emdvalue ]=IF_hilbert(IMF(h,:),fs);
    plot(t,f);
end
disp(strcat('frequency in fft from IMF2 to IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_fft)
disp(strcat('frequency in emd not uing heilbert from IMF2 to    IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_emd)
disp(strcat('frequency in emd using hilbert from IMF2 to    IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_emd2)

运行之后会生成多幅图片,对于 15Hz 的信号会得到如下的图片:

15Hz 原始信号

15Hz 各 IMF 分量时域图

15Hz 各 IMF fft 分析得到的频域图

15Hz 各 IMF 利用方法二得到的各瞬时频率在频域上的分布

15Hz 各 IMF 利用方法二得到的个瞬时频率时域分布

15Hz 各 IMF 利用方法一得到的各瞬时频率在频域上的分布

15Hz 各 IMF 利用方法一得到的个瞬时频率时域分布

对于 10Hz 的信号会得到如下的图片:

10Hz 原始信号

10Hz 各 IMF 分量时域图

10Hz 各 IMF fft 分析得到的频域图

10Hz 各 IMF 利用方法二得到的各瞬时频率在频域上的分布

10Hz 各 IMF 利用方法二得到的个瞬时频率时域分布

10Hz 各 IMF 利用方法一得到的各瞬时频率在频域上的分布

10Hz 各 IMF 利用方法一得到的个瞬时频率时域分布

在对于数据特征提取的时候,对于 fft,我们直接找频域图(就是第二幅图)上的最大值所对应的频率就是信号的频率,而对于方法一和方法二,我们对于所有瞬时频率进行统计(就是第三和第五幅图,我这里是从1到100,以1为单位进行统计,也就是统计瞬时频率为0.5~1.5的点数,然后统计瞬时频率为1.5~2.5的点数,以此类推),范围内点数最多的频率就是此信号的频率。所以每个频率下的六幅图中,我们重点关注第二、第三和第五幅图,会发现上面两种方法的结果的区别并不明显,并且显然 fft 的结果比上面两种方法的结果还明显。再看一下程序输出的结果,对于15Hz,fft显示频率为14.9414Hz,方法一的结果为13Hz,方法二的结果为13Hz。对于10Hz,fft显示的频率为9.9606Hz,方法一为10Hz,方法二为10Hz。很显然,方法一二在15Hz的时候变得极为不准确。再看图,会发现15Hz的时候方法一得到的瞬时频率时域图(第六幅图)竟然出现了负频,但是 IMF 应该是不存在负频的,而且边缘出现突变,说明算法并不完美,还需要解决这两个问题,而且应该尝试更多的数据并缩短时间(现在我使用的是全部的6秒钟数据)来看看结果如何。

最后总结一下,尽管尝试了两种方法,但是和论文上看到的结果相去甚远,下一步需要:

  1. 解决两个已发现问题:

    • 15Hz的时候方法一得到的瞬时频率时域图出现负频的问题
    • 瞬时频率时域图边缘突变
  2. 缩短时间看结果


Comments