用希尔伯特黄变换(HHT)求时频谱和边际谱

精品文档 1欢迎下载 用希尔伯特黄变换 HHT 求时频谱和边际谱 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 对求上述信号的幅频和相频曲线 function fftfenxi clear clc N 2048 fft 默认计算的信号是从 0 开始的 t linspace 1 2 N deta t 2 t 1 1 deta x 5 sin 2 pi 10 t 5 sin 2 pi 35 t N1 256 N2 512 w1 0 2 2 pi w2 0 3 2 pi w3 0 4 2 pi x t 200 y x m 0 N 1 f 1 N deta m 可以查看课本就是这样定义横坐标频率范围的 下面计算的 Y 就是 x t 的傅里叶变换数值 Y exp i 4 pi f fft y 将计算出来的频谱乘以 exp i 4 pi f 得到频移后 2 2 之间 的频谱值 Y fft y 精品文档 2欢迎下载 z sqrt Y conj Y plot f 1 100 z 1 100 title 幅频曲线 xiangwei angle Y figure 2 plot f xiangwei title 相频曲线 figure 3 plot t y r axis 2 2 0 1 2 title 原始信号 复制代码 2 用 Hilbert 变换直接求该信号的瞬时频率 clear clc clf 假设待分析的函数是 z t 3 N 2048 fft 默认计算的信号是从 0 开始的 t linspace 1 2 N deta t 2 t 1 fs 1 deta x 5 sin 2 pi 10 t 5 sin 2 pi 35 t z x hx hilbert z xr real hx xi imag hx 计算瞬时振幅 sz sqrt xr 2 xi 2 计算瞬时相位 精品文档 3欢迎下载 sx angle hx 计算瞬时频率 dt diff t dx diff sx sp dx dt plot t 1 N 1 sp title 瞬时频率 复制代码 小结 傅里叶变换不能得到瞬时频率 即不能得到某个时刻的频率值 Hilbert 变换是求 取瞬时频率的方法 但如果只用 Hilbert 变换求出来的瞬时频率也不准确 出现负频 实 际上负频没有意义 3 用 HHT 求取信号的时频谱与边际谱 function HHT clear clc clf N 2048 fft 默认计算的信号是从 0 开始的 t linspace 1 2 N deta t 2 t 1 fs 1 deta x 5 sin 2 pi 10 t 5 sin 2 pi 35 t z x c emd z 计算每个 IMF 分量及最后一个剩余分量 residual 与原始信号的相关性 m n size c for i 1 m a corrcoef c i z xg i a 1 2 end 精品文档 4欢迎下载 xg for i 1 m 1 计算各 IMF 的方差贡献率 定义 方差为平方的均值减去均值的平方 均值的平方 imfp2 mean c i 2 2 平方的均值 imf2p mean c i 2 2 各个 IMF 的方差 mse i mean c i 2 2 mean c i 2 2 end mmse sum mse for i 1 m 1 mse i mean c i 2 2 mean c i 2 2 方差百分比 也就是方差贡献率 mseb i mse i mmse 100 显示各个 IMF 的方差和贡献率 end 画出每个 IMF 分量及最后一个剩余分量 residual 的图形 figure 1 for i 1 m 1 disp imf int2str i disp mse i mseb i end subplot m 1 1 1 精品文档 5欢迎下载 plot t z set gca fontname times New Roman set gca fontsize 14 0 ylabel signal Amplitude for i 1 m 1 subplot m 1 1 i 1 set gcf color w plot t c i k set gca fontname times New Roman set gca fontsize 14 0 ylabel imf int2str i end subplot m 1 1 m 1 set gcf color w plot t c m k set gca fontname times New Roman set gca fontsize 14 0 ylabel r int2str m 1 画出每个 IMF 分量及剩余分量 residual 的幅频曲线 figure 2 subplot m 1 1 1 set gcf color w f z fftfenxi t z plot f z k set gca fontname times New Roman 精品文档 6欢迎下载 set gca fontsize 14 0 ylabel initial signal int2str m 1 Amplitude for i 1 m 1 subplot m 1 1 i 1 set gcf color w f z fftfenxi t c i plot f z k set gca fontname times New Roman set gca fontsize 14 0 ylabel imf int2str i Amplitude end subplot m 1 1 m 1 set gcf color w f z fftfenxi t c m plot f z k set gca fontname times New Roman set gca fontsize 14 0 ylabel r int2str m 1 Amplitude hx hilbert z xr real hx xi imag hx 计算瞬时振幅 sz sqrt xr 2 xi 2 计算瞬时相位 sx angle hx 计算瞬时频率 精品文档 7欢迎下载 dt diff t dx diff sx sp dx dt figure 6 plot t 1 N 1 sp title 瞬时频率 计算 HHT 时频谱和边际谱 A fa tt hhspectrum c E tt1 toimage A fa tt length tt figure 3 disp hhs E tt1 二维图显示 HHT 时频谱 E 是求得的 HHT 谱 pause figure 4 for i 1 size c 1 faa fa i FA TT1 meshgrid faa tt1 三维图显示 HHT 时频图 surf FA TT1 E title HHT 时频谱三维显示 hold on end hold off E flipud E for k 1 size E 1 bjp k sum E k 1 fs end 精品文档 8欢迎下载 f 1 N 2 N fs 2 figure 5 plot f bjp xlabel 频率 Hz ylabel 信号幅值 title 信号边际谱 要求边际谱必须先对信号进行 EMD 分解 function A f tt hhspectrum x t l aff error nargchk 1 4 nargin if nargin 2 t 1 size x 2 end if nargin 3 l 1 end if nargin 4 aff 0 end if min size x 1 if size x 2 1 x x if nargin 0 error inf doit etre 0 end M max max im im log10 im M 1e 300 inf inf 10 imagesc t fliplr 1 size im 1 2 size im 1 im inf 0 set gca YDir normal xlabel time ylabel normalized frequency 精品文档 11欢迎下载 title Hilbert Huang spectrum function f z fftfenxi t y L length t N 2 nextpow2 L fft 默认计算的信号是从 0 开始的 t linspace t 1 t L N deta t 2 t 1 m 0 N 1 f 1 N deta m 下面计算的 Y 就是 x t 的傅里叶变换数值 Y exp i 4 pi f fft y 将计算出来的频谱乘以 exp i 4 pi f 得到频移后 2 2 之间 的频谱值 Y fft y z sqrt Y conj Y 复制代码 4 总结 1 边际谱与傅里叶谱的比较 意义不同 边际谱从统计意义上表征了整组数据每个频率点的累积幅值分布 而 傅里叶频谱的某一点频率上的幅值表示在整个信号里有一个含有此频率的三角函数组分 作用不同 边际谱可以处理非平稳信号 如果信号中存在某一频率的能量出现 就表示一定有该频率的振动波出现 也就是说 边际谱能比较准确地反映信号的实际频率 成分 而傅里叶变换只能处理平稳信号 2 HHT 与 Hilbert 变换的比较 Hilbert 变换只是单纯地求