为您推荐:
精华内容
最热下载
问答
  • 5星
    7.5MB xj1234chen 2019-03-28 14:22:50
  • 5.92MB xj1234chen 2019-03-28 14:20:22
  • 1KB weixin_42696271 2021-09-11 06:14:25
  • 11.58MB QIANG654001342 2021-01-12 05:01:28
  • 将数据先进行fast-ICA处理,再进行滤波之后得到的数据,进行HHT变换。 导入库: import numpy as np import pyhht #pyhht库,进行hht变换的主要库,可使用pip install pyhht安装 import matplotlib.pyplot as plt ...

    将数据先进行fast-ICA处理,再进行滤波之后得到的数据,进行HHT变换。
    导入库:

    import numpy as np
    import pyhht    #pyhht库,进行hht变换的主要库,可使用pip install pyhht安装
    import matplotlib.pyplot as plt
    from pyhht.visualization import plot_imfs    #同pyhht库,可视化操作
    import tftb.processing    #用来计算瞬时频率的包
    from scipy.signal import hilbert    #hilbert变换需要使用的包
    

    导入数据:

    LoadList = ["EEG-Logs/2-滤波后/EEG_1_0.txt", "EEG-Logs/2-滤波后/EEG_1_2.txt", "EEG-Logs/2-滤波后/EEG_1_3.txt"]   #由于可能需要导入多个数据,故使用List表示
    Data = np.loadtxt(LoadList[0], delimiter=',')   #data为一个(263,14)的数组
    

    进行hht处理:

    DataRaw = Data[:, 0]  # 原信号
    decomposer = pyhht.emd.EMD(DataRaw)
    imfs = decomposer.decompose()      #获取到imfs数据。
    

    画图:

    plot_imfs(DataRaw, imfs)    #得到7张图,一张原数据,一张IMF以及rn数据
    

    在这里插入图片描述
    再使用matplotlib库的画图可以知道imfs[5,:]就是rn

    for i in range(6):
        plt.subplot(6, 1, i + 1)
        plt.plot(imfs[i, :])
    plt.show()
    

    在这里插入图片描述
    根据hht的原理,可以知道,原始信号x(t)可以有如下表达:
    x ( t ) = ∑ i = 1 n c i ( t ) + r n ( t ) x(t)=\sum_{i=1}^{n}c_i(t)+r_n(t) x(t)=i=1nci(t)+rn(t)
    根据上图可知,ifms这个数组里面的前5行代表 c i ( t ) c_i(t) ci(t),即IMF分量,最后1行代表 r n ( t ) r_n(t) rn(t),即剩余信号。
    理论上如此,但是实际上把IMF和剩余信号相加,最后减去原信号,还是有一定的误差。误差如下:

    rn = DataRaw - ImfSum  # 误差
    plt.plot(rn)
    plt.show()
    

    在这里插入图片描述
    可以看到,误差为 1 − 15 1^{-15} 115级别。

    接下来要对数据IMF分量进行HT变换。

    ReImf0 = hilbert(imfs[0])
    ReImf1 = hilbert(imfs[1])
    ReImf2 = hilbert(imfs[2])
    ReImf3 = hilbert(imfs[3])
    ReImf4 = hilbert(imfs[4])
    ReImf = np.array([ReImf0, ReImf1, ReImf2, ReImf3, ReImf4])
    

    变换之后放进ReImf数组里面。

    amplitude = abs(ReImf0)     #计算幅值
    instf, timestamps = tftb.processing.inst_freq(ReImf0)    #计算瞬时频率
    plt.plot(timestamps, instf)     #画图看看效果
    plt.show()
    

    瞬时频率
    包络线以及IMF
    至此,HHT变换完毕,接下来要对获得的IMF0和IMF1构建AR参数,书上推荐使用Yule-Walker方程,但是我不会。。。网上找了最小二乘法的代码,放了进来

    def ar_least_square(sample, p):
        matrix_x = np.zeros((sample.size - p, p))
        matrix_x = np.matrix(matrix_x)
        array = sample.reshape(sample.size)
        j = 0
        for i in range(sample.size - p):
            matrix_x[i, 0:p] = array[j:j + p]
            j = j + 1
        matrix_y = np.array(array[p:sample.size])
        matrix_y = matrix_y.reshape(sample.size - p, 1)
        matrix_y = np.matrix(matrix_y)
        # cofe为AR系数
        cofe = np.dot(np.dot((np.dot(matrix_x.T, matrix_x)).I, matrix_x.T), matrix_y)
        return np.array(cofe)
    

    返回AR系数

    def cal(DataRaw):
        DataRaw = DataRaw  # 原信号
        decomposer = pyhht.emd.EMD(DataRaw)
        imfs = decomposer.decompose()
        ReImf0 = hilbert(imfs[0])
        ReImf1 = hilbert(imfs[1])
        ReImf2 = hilbert(imfs[2])
        ReImf3 = hilbert(imfs[3])
        ReImf4 = hilbert(imfs[4])
        instf0, timestamps0 = tftb.processing.inst_freq(ReImf0)
        instf1, timestamps1 = tftb.processing.inst_freq(ReImf1)
        instf2, timestamps2 = tftb.processing.inst_freq(ReImf2)
        instf3, timestamps3 = tftb.processing.inst_freq(ReImf3)
        instf4, timestamps4 = tftb.processing.inst_freq(ReImf4)
        IA0 = abs(ReImf0)
        IA1 = abs(ReImf1)
        instf = np.array([instf0, instf1, instf2, instf3, instf4])  # 频率均值
        AR0 = ar_least_square(IA0, 4)
        AR1 = ar_least_square(IA1, 4)
        f1 = instf.mean()
        return np.vstack((AR0, AR1, f1))    #纵向拼接
    

    然后把前面部分的代码写成函数,并且把AR0,AR1,f1,组成一个(9,1)的向量,其中f1本来应该使用
    I E ( t ) = ∫ ω 1 ω 2 H 2 ( ω , t ) d ω IE(t)=\int_{\omega 1}^{\omega 2}H^2(\omega ,t)d\omega IE(t)=ω1ω2H2(ω,t)dω
    公式计算的,但是找不到计算这个的库,于是把频率取了平均值组成了向量。。。

    val0 = np.zeros((126,1))
    val1 = np.zeros((126,1))
    val2 = np.zeros((126,1))
    Data = np.loadtxt(LoadList[0], delimiter=',')
    for i in range(14):
        val0[9 * i:9 * i + 9] = cal(Data[:, i])
    
    Data = np.loadtxt(LoadList[1], delimiter=',')
    for i in range(14):
        val1[9 * i:9 * i + 9] = cal(Data[:, i])
    
    Data = np.loadtxt(LoadList[2], delimiter=',')
    for i in range(14):
        val2[9 * i:9 * i + 9] = cal(Data[:, i])
    
    val=np.hstack((val0,val1,val2))    #横向拼接
    print(val)
    plt.plot(val[:,0])
    plt.plot(val[:,1])
    plt.plot(val[:,2])
    plt.show()
    

    三个数组对比
    感觉差了很多。然后把这个126*3的数组存一下

    np.savetxt("EEG-Logs/3-特征提取后/EEG1.txt", val, delimiter=',')
    
    展开全文
    qq_34769201 2019-03-11 21:15:08
  • 5KB weixin_38549721 2021-05-29 07:30:11
  • HHT就是先将信号进行经验模态分解(EMD分解),然后将分解后的每个IMF分量进行Hilbert变换,得到信号的时频属性的一种时频分析方法。 2.EMD分解的步骤。   EMD分解的流程图如下:   3.实例演示。 给定频率...
    1.什么是HHT?
    
    HHT就是先将信号进行经验模态分解(EMD分解),然后将分解后的每个IMF分量进行Hilbert变换,得到信号的时频属性的一种时频分析方法。


    2.EMD分解的步骤。
     

    EMD分解的流程图如下:
     



    3.实例演示。
    给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号,采样频率fs=2048Hz的信号,表达式如下:y=5sin(2*pi*10t)+5*sin(2*pi*35t)

    (1)为了对比,先用fft对求上述信号的幅频和相频曲线
    1. function fftfenxi
    2. clear;clc;
    3. N=2048;
    4. %fft默认计算的信号是从0开始的
    5. t=linspace(1,2,N);deta=t(2)-t(1);1/deta
    6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
    7. % N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;
    8. % x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t);
    9. y = x;
    10. m=0:N-1;
    11. f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的
    12. %下面计算的Y就是x(t)的傅里叶变换数值
    13. %Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值
    14. Y=fft(y);
    15. z=sqrt(Y.*conj(Y));
    16. plot(f(1:100),z(1:100));
    17. title('幅频曲线')
    18. xiangwei=angle(Y);
    19. figure(2)
    20. plot(f,xiangwei)
    21. title('相频曲线')
    22. figure(3)
    23. plot(t,y,'r')
    24. %axis([-2,2,0,1.2])
    25. title('原始信号')
    复制代码




    (2)用Hilbert变换直接求该信号的瞬时频率
    1. clear;clc;clf;
    2. %假设待分析的函数是z=t^3
    3. N=2048;
    4. %fft默认计算的信号是从0开始的
    5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;
    6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
    7. z=x;
    8. hx=hilbert(z);
    9. xr=real(hx);xi=imag(hx);
    10. %计算瞬时振幅
    11. sz=sqrt(xr.^2+xi.^2);
    12. %计算瞬时相位
    13. sx=angle(hx);
    14. %计算瞬时频率
    15. dt=diff(t);
    16. dx=diff(sx);
    17. sp=dx./dt;
    18. plot(t(1:N-1),sp)
    19. title('瞬时频率')
    复制代码

    小结:傅里叶变换不能得到瞬时频率,即不能得到某个时刻的频率值。Hilbert变换是求取瞬时频率的方法,但如果只用Hilbert变换求出来的瞬时频率也不准确。(出现负频,实际上负频没有意义!)
    (3)用HHT求取信号的时频谱与边际谱
    1. function HHT
    2. clear;clc;clf;
    3. N=2048;
    4. %fft默认计算的信号是从0开始的
    5. t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;
    6. x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
    7. z=x;
    8. c=emd(z);

    9. %计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性
    10. [m,n]=size(c);
    11. for i=1:m;
    12. a=corrcoef(c(i,:),z);
    13. xg(i)=a(1,2);
    14. end
    15. xg;

    16. for i=1:m-1
    17. %--------------------------------------------------------------------
    18. %计算各IMF的方差贡献率
    19. %定义:方差为平方的均值减去均值的平方
    20. %均值的平方
    21. %imfp2=mean(c(i,:),2).^2
    22. %平方的均值
    23. %imf2p=mean(c(i,:).^2,2)
    24. %各个IMF的方差
    25. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
    26. end;
    27. mmse=sum(mse);
    28. for i=1:m-1
    29. mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2; 
    30. %方差百分比,也就是方差贡献率
    31. mseb(i)=mse(i)/mmse*100;
    32. %显示各个IMF的方差和贡献率
    33. end;
    34. %画出每个IMF分量及最后一个剩余分量residual的图形
    35. figure(1)
    36. for i=1:m-1
    37. disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
    38. end;
    39. subplot(m+1,1,1)
    40. plot(t,z)
    41. set(gca,'fontname','times New Roman')
    42. set(gca,'fontsize',14.0)
    43. ylabel(['signal','Amplitude'])

    44. for i=1:m-1
    45. subplot(m+1,1,i+1);
    46. set(gcf,'color','w')
    47. plot(t,c(i,:),'k')
    48. set(gca,'fontname','times New Roman')
    49. set(gca,'fontsize',14.0)
    50. ylabel(['imf',int2str(i)])
    51. end
    52. subplot(m+1,1,m+1);
    53. set(gcf,'color','w')
    54. plot(t,c(m,:),'k')
    55. set(gca,'fontname','times New Roman')
    56. set(gca,'fontsize',14.0)
    57. ylabel(['r',int2str(m-1)])

    58. %画出每个IMF分量及剩余分量residual的幅频曲线
    59. figure(2)
    60. subplot(m+1,1,1)
    61. set(gcf,'color','w')
    62. [f,z]=fftfenxi(t,z);
    63. plot(f,z,'k')
    64. set(gca,'fontname','times New Roman')
    65. set(gca,'fontsize',14.0)
    66. ylabel(['initial signal',int2str(m-1),'Amplitude'])

    67. for i=1:m-1
    68. subplot(m+1,1,i+1);
    69. set(gcf,'color','w')
    70. [f,z]=fftfenxi(t,c(i,:));
    71. plot(f,z,'k')
    72. set(gca,'fontname','times New Roman')
    73. set(gca,'fontsize',14.0)
    74. ylabel(['imf',int2str(i),'Amplitude'])
    75. end
    76. subplot(m+1,1,m+1);
    77. set(gcf,'color','w')
    78. [f,z]=fftfenxi(t,c(m,:));
    79. plot(f,z,'k')
    80. set(gca,'fontname','times New Roman')
    81. set(gca,'fontsize',14.0)
    82. ylabel(['r',int2str(m-1),'Amplitude'])

    83. hx=hilbert(z);
    84. xr=real(hx);xi=imag(hx);
    85. %计算瞬时振幅
    86. sz=sqrt(xr.^2+xi.^2);
    87. %计算瞬时相位
    88. sx=angle(hx);
    89. %计算瞬时频率
    90. dt=diff(t);
    91. dx=diff(sx);
    92. sp=dx./dt;
    93. figure(6)
    94. plot(t(1:N-1),sp)
    95. title('瞬时频率')

    96. %计算HHT时频谱和边际谱
    97. [A,fa,tt]=hhspectrum(c);
    98. [E,tt1]=toimage(A,fa,tt,length(tt));
    99. figure(3)
    100. disp_hhs(E,tt1) %二维图显示HHT时频谱,E是求得的HHT谱
    101. pause
    102. figure(4)
    103. for i=1:size(c,1)
    104. faa=fa(i,:);
    105. [FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图
    106. surf(FA,TT1,E)
    107. title('HHT时频谱三维显示')
    108. hold on
    109. end
    110. hold off
    111. E=flipud(E);
    112. for k=1:size(E,1)
    113. bjp(k)=sum(E(k,:))*1/fs; 
    114. end
    115. f=(1:N-2)/N*(fs/2);
    116. figure(5)
    117. plot(f,bjp);
    118. xlabel('频率 / Hz');
    119. ylabel('信号幅值');
    120. title('信号边际谱')%要求边际谱必须先对信号进行EMD分解

    121. function [A,f,tt] = hhspectrum(x,t,l,aff)

    122. error(nargchk(1,4,nargin));

    123. if nargin < 2

    124. t=1:size(x,2);

    125. end

    126. if nargin < 3

    127. l=1;

    128. end

    129. if nargin < 4

    130. aff = 0;

    131. end

    132. if min(size(x)) == 1
    133. if size(x,2) == 1
    134. x = x';
    135. if nargin < 2
    136. t = 1:size(x,2);
    137. end
    138. end
    139. Nmodes = 1;
    140. else
    141. Nmodes = size(x,1);
    142. end

    143. lt=length(t);

    144. tt=t((l+1):(lt-l));

    145. for i=1:Nmodes

    146. an(i,:)=hilbert(x(i,:)')';
    147. f(i,:)=instfreq(an(i,:)',tt,l)';
    148. A=abs(an(:,l+1:end-l));

    149. if aff
    150. disprog(i,Nmodes,max(Nmodes,100))
    151. end

    152. end


    153. function disp_hhs(im,t,inf)

    154. % DISP_HHS(im,t,inf)
    155. % displays in a new figure the spectrum contained in matrix "im"
    156. % (amplitudes in log).
    157. %
    158. % inputs : - im : image matrix (e.g., output of "toimage")
    159. % - t (optional) : time instants (e.g., output of "toimage") 
    160. % - inf (optional) : -dynamic range in dB (wrt max)
    161. % default : inf = -20
    162. %
    163. % utilisation : disp_hhs(im) ; disp_hhs(im,t) ; disp_hhs(im,inf) 
    164. % disp_hhs(im,t,inf)

    165. figure
    166. colormap(bone)
    167. colormap(1-colormap);

    168. if nargin==1
    169. inf=-20;
    170. t = 1:size(im,2);

    171. end

    172. if nargin == 2
    173. if length(t) == 1
    174. inf = t;
    175. t = 1:size(im,2);
    176. else
    177. inf = -20;
    178. end
    179. end

    180. if inf >= 0
    181. error('inf doit etre < 0')
    182. end

    183. M=max(max(im));

    184. im = log10(im/M+1e-300);

    185. inf=inf/10;


    186. imagesc(t,fliplr((1:size(im,1))/(2*size(im,1))),im,[inf,0]);
    187. set(gca,'YDir','normal')
    188. xlabel(['time'])
    189. ylabel(['normalized frequency'])
    190. title('Hilbert-Huang spectrum')
    191. function [f,z]=fftfenxi(t,y)
    192. L=length(t);N=2^nextpow2(L);
    193. %fft默认计算的信号是从0开始的
    194. t=linspace(t(1),t(L),N);deta=t(2)-t(1);
    195. m=0:N-1;
    196. f=1./(N*deta)*m;
    197. %下面计算的Y就是x(t)的傅里叶变换数值
    198. %Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值
    199. Y=fft(y);
    200. z=sqrt(Y.*conj(Y));



    复制代码




    4.总结。

    (1)边际谱与傅里叶谱的比较:
              意义不同:边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布,而傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分。
             作用不同:边际谱可以处理非平稳信号,如果信号中存在某一频率的能量出现,就表示一定有该频率的振动波出现,也就是说,边际谱能比较准确地反映信号的实际频率成分。而傅里叶变换只能处理平稳信号。


    (2) HHT与Hilbert变换的比较:  
              Hilbert变换只是单纯地求信号的瞬时振幅,频率和相位,有可能出现没有意义的负频率;HHT变换先将信号进行EMD分解,得到的是各个不同尺度的分量,对每一个分量进行Hilbert变换后得到的是有实际意义的瞬时频率。
    展开全文
    heifan2014 2017-05-21 02:41:09
  • 3星
    458KB u011019894 2018-03-27 12:59:19
  • 1.45MB qq_42422055 2018-08-24 23:54:37
  • 712KB weixin_42242017 2018-05-19 10:52:04
  • 441KB weixin_38657353 2020-07-08 01:07:34
  • 浅谈傅里叶变换、小波变换、HHT变换一、傅里叶变换1.1傅里叶变换介绍二、小波变换2.1小波变换正反变换公式2.2小波变换适应场景及其优缺点2.3小波变换的应用三、HHT变换3.1HHT产生的背景3.1 HHT变换介绍3.2 HHT对信号...

    一、傅里叶变换

    1.1傅里叶变换介绍

    \quad 我们生活中常见的信息的描述基本上都是在时域空间内进行描述的,如下图1所示;但如果当我们碰到一些杂乱无章的信号需要处理时,如图二所示,我们就很难在时域空间内分析出任何有用信息。于是伟大的傅里叶提出了傅里叶变换理论,将时域空间内的信息可以转换到频域空间,并且将两个空间通过一套完整的转换公式联系起来。于是我们可以对图二的时域信号进行傅里叶变换,我们则会得到像图三(此处图三并不代表图二的频域显示图,我只是为了讲解时域到频域这一变换,还望理解)所示的信号在频域空间的分布图。
    在这里插入图片描述在这里插入图片描述在这里插入图片描述
    \quad 1822年,法国工程师傅里叶指出:一个“任意”的周期函数 x ( t ) x(t) x(t)都可以分解为无穷个不同频率正弦信号的和,即傅里叶级数。其中求解傅里叶系数的过程就是傅里叶变换。如下所示,第一个公式我们称之为傅里叶变换,将时域信号 f ( t ) f(t) f(t)在整个区间 R R R内进行积分,转换为频域信号 F ( w ) F(w) F(w)。第二个式子为傅里叶反变换,实现将频域信号 F ( w ) F(w) F(w)转换为时域信号 f ( t ) f(t) f(t).
    F ( w ) = ∫ − ∞ ∞ f ( t ) e − j w t d t F(w)=\int_{-\infty}^{\infty}f(t)e^{-jwt}dt F(w)=f(t)ejwtdt
    f ( t ) = 1 2 π ∫ − ∞ ∞ F ( w ) e j w t d w f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(w)e^{jwt}dw f(t)=2π1F(w)ejwtdw
    当然傅里叶变换的使用也是有一定限制的,比如:
    \qquad 1. 无穷区间的正弦波基函数所定义的一种整体变换
    \qquad 2. 仅适用于对信号进行全局分析,不适用于对信号进行局部分析。
    \qquad 3. 仅适用于分析频率不随时间变化的平稳信号,但不适用于分析频率随时间变化的非平稳信号。

    二、小波变换

    2.1小波变换正反变换公式

    \quad 以下我将会列出小波基变换的公式,其中第一个公式属于限制条件,第二个和第三个公式就像傅里叶正反变换的关系,可以相互转化。
    W f ( a , b ) < = f W_f(a,b)<=f Wf(a,b)<=f
    ψ a b > = 1 ∣ a ∣ ∫ R f ( t ) ψ ( t − b a ) ‾ d t \psi_{ab}>=\frac{1}{\sqrt{|a|}}\int_Rf(t)\overline{\psi({\frac{t-b}{a}})}dt ψab>=a 1Rf(t)ψ(atb)dt
    f ( t ) = 1 C ψ ∫ R + ∫ R 1 a 2 W f ( a , b ) ψ ( t − b a ) d a d b f(t)=\frac{1}{C_{\psi}}\int_{R^+}\int_R\frac{1}{a^2}W_f(a,b)\psi(\frac{t-b}{a})dadb f(t)=Cψ1R+Ra21Wf(a,b)ψ(atb)dadb
    其中式中的变量代表意思如下:
    \qquad\qquad a , b a,b a,b: \qquad 伸缩和平移因子,
    \qquad\qquad t − b t-b tb : :\quad :时间上或者空间上的平移,
    \qquad\qquad t a \frac{t}{a} at: \qquad 尺度上或频率上的伸缩。

    还需要注意的是 ψ ( t − b a ) {\psi({\frac{t-b}{a}})} ψ(atb)被称为小波基,它是可变的,但是它对信号没有自适应性 ψ ( t − b a ) ‾ \overline{\psi({\frac{t-b}{a}})} ψ(atb)是小波基的复数求共轭,它和小波基一样,都可以变化,但是对信号也没有自适应性。

    2.2小波变换适应场景及其优缺点

    1. 适用于非平稳信号的分析
    2. 小波基函数的选取十分重要,基函数的选取的不同可能会造成分析结果的不一致,分析结果的准确性取决于选取合适的小波基函数。
    3. 最优小波基的选取方法研究。现在国内外已经有一些最优基选取方法,但是缺乏系统规范的最佳小波基的选取方法,即针对不同的问题能最优的选择不同的小波基,以实现好的应用效果。
    4. 不存在一种小波基能使适用所有的情况。因此小波基的最优化选择始终是小波理论研究的重要内容。

    2.3小波变换的应用

    小波变换的应用广泛,目前主要应用领域有以下所列:

    1. 地震信号分析
    2. 连续小波用于漩涡研究
    3. 二进制小波用于图像的边缘检测、图像压缩和重构。
    4. 小波包用于图像压缩
    5. 噪声的未知瞬态信号。
    6. 语音信号处理
    7. 时频分析
    8. 正交小波用于算子和微分算子的简化。

    三、HHT变换

    3.1HHT产生的背景

    \quad 1. 频率的表示方法是建立在傅里叶变换的基础上的,由于傅里叶变换是一种全局的变换,要么完全在时域,要么完全在频域,因此无法表述信号的时频局部性质,而时频局部性质恰好是非平稳信号最基本和最关键的性质
    \quad 2. 虽然小波变换却具有自动改变窗口长短的功能,可以很好的把信号在时间上和频域上局部化,从而呈现了信号的局部奇异性。从分辨率看,小波变换很好的解决了时间和频率的分辨的矛盾,它在低频段用高的频率分辨率和低的时间分辨率;而在高频段使用低的频率分辨率和高的时间分辨率。(事实上,它也只能这样。)小波变换的窗宽变换是自适应的。(注意此处是窗宽,也就是矩形窗的宽度是可以自适应的,与之前说的小波基不能对信号的自适应是不矛盾的。)
    \qquad 3. 就这样HHT(Hilbert-Huang Transfrom),也称作希尔伯特-黄,在结合了傅里叶变换和小波变换的基础上诞生了。它是由美籍华裔科学家黄锷,创立的一种研究方法。它既吸取了小波变换的多分辨率的优势,又克服了小波变换中需要选择小波基的困难。

    3.2 HHT变换介绍

    \quad 1. HHT是用于过程采样、可描述和仿真非平稳过程的一种非线性分析新方法。它通过EMD经验模式分解,将信号分解成有限个数的 I M F IMF IMF信号,并对每个 I M F IMF IMF信号进行Hilert变换,就可以获得有意义的瞬时功率,从而给出频率变化的精确表达。
    HHT自适应的利用了信号的局部信息,从而获得信号某一时刻的顺势状态。
    \quad 2. 对于HHT方法,它由 EMD 分解和 Hilbert 变换两部分构成,核心 是 EMD分解。EMD方法可以提取信号的瞬时频率和瞬时能量参数,这是实现信号瞬时分析的有效方法。
    \quad 3. 接下来我将分别介绍HHT对于实现对信号瞬时分析的分析框图和EMD经验模式分解的基本原理。

    3.3 HHT对信号分析的框图

    在这里插入图片描述

    3.4 EMD经验模式分解的基本原理

    \qquad 为了得到有意义的瞬时频率,黄锷等人提出了在物理上得到一个有意义的瞬时频率的必要条件是:
    \qquad 函数对称于局部零均值且有相同的极值和过零点。并且根据这些条件进一步定义了满足以下两个条件的函数称为固有模态函数( I M F IMF IMF),且这类函数的任一点都存在一个有意义的瞬时频率:
    1. 信号上任意一点,由局部极大值点确定的包络线与由局部极小值点确定的包络线的平均值都为0,也就是说信号关于时间轴局部对称。
    2. 在整个离散信号序列中,极值点的个数与过零点的个数相等或最多相差1。

    在这里插入图片描述
    \qquad 以上图来说,数字从1-5中,有三个过零点,分别为1,2,5;有两个极值点,分别为极大值点3和极小值点4。按照上边的约定我们发现零点个数减去极值点个数为3-2=1,正好符合上边的极值点的个数与过零点的个数相等或最多相差1。并且条件1也满足。所以我们可以说,数字1-5所包含的信号则可以被认为成一个周期。 I M F IMF IMF存 在有意义的瞬时频率,可通过 对该分离出来的周期信号进行Hilbert 变换求得IMF。
    \qquad 由于 EMD 分解的分解基来自原信号本身,因此它具有自适应性,而不象傅里叶变换那样把信号分解为固定的正弦或余弦函数之和的形式。由此可见,EMD 分解法是一种基于信号本身的时间尺度特征的时域处理方法,它 把复杂的信号分解为不同尺度特征的本征模态函数( I M F IMF IMF)之和的形式,并且每个模函数分量上任一点都存在有意义的瞬时频率(通过 Hilbert 变换可求出)。
    EMD过程的具体算法如下:
    \qquad 对一原始信号 X ( t ) X(t) X(t),首先要找出 X ( t ) X(t) X(t)上所有的极值点,然后用三次样条函数曲线对所有的极大值点进行插值,从而拟合出原始信号 X ( t ) X(t) X(t)上的包络线 X m a x ( t ) X_{max}(t) Xmax(t),同理得到下包络线 X m i n ( t ) X_{min}(t) Xmin(t)。上、下两条包络线包含了所有的信号数据。按顺序连接上、下两条包络线的均值即可得到一条均值线 m 1 ( t ) m_1(t) m1(t):
    m 1 ( t ) = X m a x ( t ) + X m i n ( t ) 2 m_1(t)=\frac{X_{max}(t)+X_{min}(t)}{2} m1(t)=2Xmax(t)+Xmin(t)
    再用 X ( t ) X(t) X(t)减去 m 1 ( t ) m_1(t) m1(t)得到 h 1 ( t ) h_1(t) h1(t)
    h 1 ( t ) = X ( t ) − m 1 ( t ) h_1(t)=X(t)-m_1(t) h1(t)=X(t)m1(t)
    \qquad 对于不同的信号, h 1 ( t ) h_1(t) h1(t)可能是一个 I M F IMF IMF分量,也可能不是,一般来说它并不满足 I M F IMF IMF所需的条件,此时将 h 1 ( t ) h_1(t) h1(t)当作原信号,重复上述步骤,则有:
    h 1 k ( t ) = h 1 ( k − 1 ) ( t ) − m 1 ( k − 1 ) ( t ) h_{1k}(t)=h_{1{(k-1)}}(t)-m_{1(k-1)}(t) h1k(t)=h1(k1)(t)m1(k1)(t)
    \qquad 对于 h 1 k h_{1k} h1k是不是一个 I M F IMF IMF分量,我们必须要有一个筛选过程终止的原则,它可以利用两个连续的处理结果之间的标准差SD作为判断依据:
    S D = ∑ t = 0 T ∣ ∣ h 1 ( k − 1 ) ( t ) − h 1 k ( t ) ∣ 2 h 1 ( k − 1 ) 2 ( t ) ∣ SD=\sum_{t=0}^T\mid{\frac{{|h_{1(k-1)}(t)-h_{1k}(t)|}^2}{h_{1(k-1)}^2(t)}}\mid SD=t=0Th1(k1)2(t)h1(k1)(t)h1k(t)2
    \qquad 决定筛选过程是否终止,SD取值一定要谨慎,既要避免过于严格的准则,以导致 I M F IMF IMF分量变成纯粹的频率调制信号,造成幅值恒定;又要避免过于宽松的准则,从而产生与要求的IMF分量相差太远的分量。实际过程中可以通过对信号反复用筛选过程而取不同的SD值来最终确定,经验表明,SD值取0.2-0.3时为宜。这样既可以保证 I M F IMF IMF的线性和稳定性,又可以IMF具有相应的物理意义。
    h 1 k ( t ) h_{1k}(t) h1k(t)满足SD的值要求时,则称 h 1 k ( t ) h_{1k}(t) h1k(t)为第一阶IMF,记为 c 1 ( t ) c_1(t) c1(t):
    c 1 ( t ) = h 1 k ( t ) c_1(t)=h_{1k}(t) c1(t)=h1k(t)
    从原信号 X ( t ) X(t) X(t)中减去 c 1 ( t ) c_1(t) c1(t)得到剩余信号,即残差 r 1 ( t ) r_1(t) r1(t)
    r 1 ( t ) = X ( t ) − c 1 ( t ) r_1(t)=X(t)-c_1(t) r1(t)=X(t)c1(t)
    \qquad 然后将 r 1 ( t ) r_1(t) r1(t)看作一组新的“原信号”,重复上述的模态分解过程,这样经过多次运算即可得到全部的残差 r i ( t ) r_i(t) ri(t):
    r i ( t ) = r i − 1 ( t ) − c i ( t ) i = 2 , 3 , 4... n r_i(t)=r_{i-1}(t)-c_i(t)\qquad\qquad{i=2,3,4...n} ri(t)=ri1(t)ci(t)i=2,3,4...n
    r i ( t ) ( i = 1 , 2 , 3 , . . . n ) r_i(t)\quad({i=1,2,3,...n)} ri(t)(i=1,2,3,...n)满足条件:
    \qquad \qquad 1. c n ( t ) c_n(t) cn(t) r n ( t ) r_n(t) rn(t)小于预定的误差;
    \qquad \qquad 2.残差 r n ( t ) r_n(t) rn(t)成为一单调函数,即不能再从中提取出IMF分量时

    满足以上两条件之一的,就可以终止模态分解过程。
    \qquad 应该注意的是,该条件的选取也应该适中,若条件太严格,则得到的最后几个 I M F IMF IMF分量没有太大意义,并且还消耗时间;若条件太松,则会丢失有用信号分量。具体终止条件的选取可以通过对信号的反复分解并依据对于原始信号的知识来最终确定。因为我们只是对原始信号 X ( t ) X(t) X(t)进行拆解,而从未抛弃原始信号的数据。因此我们可以通过累加的形式再次复原出我们的原始信号 X ( t ) X(t) X(t),它可由n阶 I M F IMF IMF分量以及残差构成:
    X ( t ) = ( ∑ i = 1 n c i ( t ) ) + r n ( t ) X(t)=(\sum_{i=1}^nc_i(t))+r_n(t) X(t)=(i=1nci(t))+rn(t)
    \qquad 最后得到的频率成分从高到低的一系列本征模态分量( I M F IMF IMF)以及最后的残留分量。由于每个IMF分量代表一组时间特征尺度的数据序列,并且以不同的分辨率显示信号特征,所以整个分解过程实际上是将原始信号分解为各个不同时间特征尺度波动的叠加。

    致谢

    此次博客内容的书写,主要是在听取了昨日的昆明理工大学信自学院的刘增力教授的关于傅里叶变换、小波变换以及HHT变换讲座的内容,在课后我就将听到的以及理解的东西发到了博客上,做一下讲座报告的记录吧。如果博客中我理解的有什么偏差甚至是错误的地方,还请各位大佬指正,在这里小kingback先行谢谢大家了。

    展开全文
    weixin_38468077 2019-08-29 12:18:42
  • 压缩包 : f914a6a90d345a26f732d9223e682699.rar 列表复件 HHT变换的三种方法 Matlab/G Rilling/document.doc复件 HHT变换的三种方法 Matlab/G Rilling/emd.ppt复件 HHT变换的三种方法 Matlab/G Rilling/EMD程序使用...

    压缩包 : f914a6a90d345a26f732d9223e682699.rar 列表

    复件 HHT变换的三种方法 Matlab/G Rilling/document.doc

    复件 HHT变换的三种方法 Matlab/G Rilling/emd.ppt

    复件 HHT变换的三种方法 Matlab/G Rilling/EMD程序使用.doc

    复件 HHT变换的三种方法 Matlab/G Rilling/pack_emd.tar.gz

    复件 HHT变换的三种方法 Matlab/G Rilling/pack_emd.zip

    复件 HHT变换的三种方法 Matlab/G Rilling/tftb-0.2.tar.gz

    复件 HHT变换的三种方法 Matlab/G Rilling/~$cument.doc

    复件 HHT变换的三种方法 Matlab/G Rilling/~WRL2980.tmp

    复件 HHT变换的三种方法 Matlab/G Rilling/~WRL3223.tmp

    复件 HHT变换的三种方法 Matlab/G Rilling/~WRL3695.tmp

    复件 HHT变换的三种方法 Matlab/G Rilling/归一化频率.mht

    复件 HHT变换的三种方法 Matlab/G Rilling/最原始的EMD程序.mht

    复件 HHT变换的三种方法 Matlab/G Rilling/正弦插值法解决emd的端点.mht

    复件 HHT变换的三种方法 Matlab/G Rilling/镜像延拓.mht

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/MATLAB Central - File detail - Hilbert-Huang Transform.mht

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/emd.m

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/findpeaks.m

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/HHT.pdf

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/Hum.wav

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/plot_hht.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/dist_value.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/eemd.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help/gsta.dat

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help/help.doc

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/extrema.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/HHT MATLAB PROGRAM.mht

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/ifndq.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/readme.mht

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/significance.m

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/state of art.mht

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/tutorial.mht

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/中央大學數據分析中心.mht

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help

    复件 HHT变换的三种方法 Matlab/G Rilling

    复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center

    复件 HHT变换的三种方法 Matlab/National Taiwan Central University

    复件 HHT变换的三种方法 Matlab

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/emd.ppt

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/pack_emd.tar.gz

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/pack_emd.zip

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/tftb-0.2.tar.gz

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/~$cument.doc

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/~$MD程序使用.doc

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/~WRL2980.tmp

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/~WRL3223.tmp

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/~WRL3695.tmp

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/归一化频率.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/最原始的EMD程序.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/正弦插值法解决emd的端点.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling/镜像延拓.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/MATLAB Central - File detail - Hilbert-Huang Transform.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/emd.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/findpeaks.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/HHT.pdf

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/Hum.wav

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht/plot_hht.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/dist_value.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/eemd.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help/gsta.dat

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help/~$help.doc

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/extrema.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/HHT MATLAB PROGRAM.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/ifndq.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/readme.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/significance.m

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/state of art.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/tutorial.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/中央大學數據分析中心.mht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center/plot_hht

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University/example and help

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/G Rilling

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/Matlab File Exchange Center

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab/National Taiwan Central University

    HHT变换的三种方法+Matlab/复件 HHT变换的三种方法 Matlab

    HHT变换的三种方法+Matlab

    展开全文
    weixin_32751907 2021-04-21 02:41:58
  • qq_37385720 2020-10-13 11:00:13
  • 146B sinat_15156907 2014-05-20 19:44:10
  • 5KB weixin_42681774 2021-10-02 16:10:57
  • 4星
    53KB jhkg0000 2012-06-04 15:16:34
  • 1.64MB u014243064 2014-03-21 16:19:42
  • 11.54MB csdn14512 2014-07-26 16:07:20
  • weixin_42356073 2021-04-21 22:07:08
  • 5.94MB feijiaogu7393 2021-05-15 11:27:24
  • 475KB u013883025 2021-08-15 11:15:41
  • 97KB weixin_39840387 2019-08-13 06:36:01
  • 746KB weixin_38529293 2020-05-30 23:20:57
  • weixin_35560840 2021-04-20 06:30:02
  • 4星
    995KB lightsalt 2012-02-06 14:39:13
  • 4星
    6MB wyf198403 2012-02-16 20:43:47
  • 362KB jiebing2020 2021-09-27 21:47:13

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 623
精华内容 249
关键字:

hht变换

友情链接: Clustering-master.rar