HHT边际谱的定义以及理解

绿荷雨丝

（1)我对HHT边际谱的定义理解为就是将时频分布对时间求积分，具体地说就是将频率相同的时刻所对应的幅值进行累加。这样理解正确吗？
（2) 如果我的理解是正确的，那么写程序时关键的一步就是判断在各种IMF下哪些时刻的频率是相同的。在实际计算中，计算机的精度导致了各种IMF下各时刻的频率值几乎不可能是完全一致的，那么我要用什么样的算法来做边际谱呢？哪位高手能提供一下具体的算法或者有现成的code可以下载？

[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]

AaronSpark

Marginal spectrum是对Hilbert spectrum即H(A,w,t)进行时间积分得到的。Hilbert spectrum是指H(A,w,t)时、频、幅三维图。

simon21

Fourier变换是将任何信号分解为正弦信号的加权和,而每一个正弦信号对应着一个固定的频率(Fourier频率)和固定的幅值,因此,用Fourier 变换分析频率不随时间变化的平稳信号是十分有效的。但对于频率随时间变化的非平稳信号,Fourier 变换就无能为力了。

HHT是历史上首次对Fourier变换的基本信号和频率定义作的创造性的改进。他们不再认为组成信号的基本信号是正弦信号,而是一种称为固有模态函数的信号,也就是满足以下两个条件的信号: (1) 整个信号中,零点数与极点数相等或至多相差1 ; (2) 信号上任意一点,由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零,即信号关于时间轴局部对称。

无论Hilbert谱中的频率还是边际谱中的频率(即瞬时频率) ,其意义都与Fourier分析中的频率(即Fourier 频率) 完全不同,但在Fourier分析中,某一频率处能量的存在,代表一个正弦或余弦波在整个时间轴上的存在,而边际谱h中某一频率处能量的存在仅代表在整个时间轴上可能有这样一个频率的振动波在局部出现过,h越大,代表该频率出现的可能性越大。

[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]

readj

目前最好的贴

追风鸟

我正急需emd分解的matlab的程序呢，meil32的回贴借了我的燃眉之急。

MVH

HHT问题很多，下面是别人提出来的
(1)因为分解的不稳定，分解所得的固有模态函数并非是系统的固有(一定)模态。
(2)且正因为分解的不稳定，分解结果一般也就很难用确切的物理过程来解释了。
(3)此外，再对一般情况下不稳定分解所得的本来就无明确物理意义的分量进行诸如波内调频(intra-wave frequency modulation)的细致分析，无疑亦是在做多此一举的无用功了。

[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]

liuyun

to：meil32，你的源程序里好像缺一个函数吧，extr（m），缺了这个不能正常运行啊，麻烦再传一下了，谢谢！

wangxiaoyu

对你的程序里却两个函数，是不能正常运行的，你要是有的话请传上来，或发到邮箱wugeren@126.com

[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]

chinawang

HHT著作已经出版，那位高人去搞个电子版的 ！

http://www.worldscibooks.com/mathematics/5862.html
http://www.amazon.com/gp/product ... s&v=glance&n=283155

[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]

wangxiaoyu

这两本书那位高手能介绍买的地方，谢谢了！！！

梦里花落

非常感谢你的程序

[ 本帖最后由 zhlong 于 2007-10-15 22:29 编辑 ]

电光石火

谢谢，我正在求这个程序，谢谢！

simon21

其中缺少的extr试试下面的代码

1. %extr
   
2. function [indmin,indmax,indzer]=extr(x)
   
3. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
   
4. watchon;
   
5. 
   
6. FILENAME=[DATApath,DATAfile];
   
7. x=load(FILENAME);
   
8. 
   
9. numzer=1;
   
10. nummin=1;
    
11. nummax=1;
    
12. for g=2:(length(x)-1)
    
13. 
    
14. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
    
15. indzer(numzer)=g-1;
    
16. numzer=numzer+1;
    
17. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
    
18. indzer(numzer)=g;
    
19. numzer=numzer+1;
    
20. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
    
21. indzer(numzer)=(2*g-1)/2;
    
22. numzer=numzer+1;
    
23. end
    
24. 
    
25. if (x(g-1)>x(g))&(x(g)>x(g+1))
    
26. continue;
    
27. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
    
28. nummin=nummin+1;
    
29. elseif (x(g-1) continue;
    
30. elseif (x(g-1)x(g+1))
    
31. indmax(nummax)=g;
    
32. nummax=nummax+1;
    
33. end
    
34. end
    
35. %subplot(2,3,1)
    
36. plot(indzer,x(indzer),'r-+')
    
37. title('indzer')
    
38. pause
    
39. % subplot(2,3,2)
    
40. plot(indmin,x(indmin),'b-o')
    
41. title('indmin')
    
42. pause
    
43. %subplot(2,3,3)
    
44. plot(indmax,x(indmax),'g-*')
    
45. title('indmax')
    
46. pause
    
47. %subplot(2,3,4)
    
48. plot(1:length(x),x,'k-s')
    
49. title('x')
    
50. pause
    
51. 
    
52. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
    
53. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
    
54. %subplot(2,3,5)
    
55. plot(1:length(x),envmin,'y-d')
    
56. title('envmin')
    
57. pause
    
58. %subplot(2,3,6)
    
59. plot(1:length(x),envmax,'c-h')
    
60. title('envmax')
    
61. pause
    
62. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
    
63. title('mean((envmax+envmin)/2)')%extr
    
64. function [indmin,indmax,indzer]=extr(x)
    
65. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
    
66. watchon;
    
67. 
    
68. FILENAME=[DATApath,DATAfile];
    
69. x=load(FILENAME);
    
70. 
    
71. numzer=1;
    
72. nummin=1;
    
73. nummax=1;
    
74. for g=2:(length(x)-1)
    
75. 
    
76. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
    
77. indzer(numzer)=g-1;
    
78. numzer=numzer+1;
    
79. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
    
80. indzer(numzer)=g;
    
81. numzer=numzer+1;
    
82. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
    
83. indzer(numzer)=(2*g-1)/2;
    
84. numzer=numzer+1;
    
85. end
    
86. 
    
87. if (x(g-1)>x(g))&(x(g)>x(g+1))
    
88. continue;
    
89. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
    
90. nummin=nummin+1;
    
91. elseif (x(g-1) continue;
    
92. elseif (x(g-1)x(g+1))
    
93. indmax(nummax)=g;
    
94. nummax=nummax+1;
    
95. end
    
96. end
    
97. %subplot(2,3,1)
    
98. plot(indzer,x(indzer),'r-+')
    
99. title('indzer')
    
100. pause
     
101. % subplot(2,3,2)
     
102. plot(indmin,x(indmin),'b-o')
     
103. title('indmin')
     
104. pause
     
105. %subplot(2,3,3)
     
106. plot(indmax,x(indmax),'g-*')
     
107. title('indmax')
     
108. pause
     
109. %subplot(2,3,4)
     
110. plot(1:length(x),x,'k-s')
     
111. title('x')
     
112. pause
     
113. 
     
114. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
     
115. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
     
116. %subplot(2,3,5)
     
117. plot(1:length(x),envmin,'y-d')
     
118. title('envmin')
     
119. pause
     
120. %subplot(2,3,6)
     
121. plot(1:length(x),envmax,'c-h')
     
122. title('envmax')
     
123. pause
     
124. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
     
125. title('mean((envmax+envmin)/2)')%extr
     
126. function [indmin,indmax,indzer]=extr(x)
     
127. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
     
128. watchon;
     
129. 
     
130. FILENAME=[DATApath,DATAfile];
     
131. x=load(FILENAME);
     
132. 
     
133. numzer=1;
     
134. nummin=1;
     
135. nummax=1;
     
136. for g=2:(length(x)-1)
     
137. 
     
138. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
     
139. indzer(numzer)=g-1;
     
140. numzer=numzer+1;
     
141. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
     
142. indzer(numzer)=g;
     
143. numzer=numzer+1;
     
144. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
     
145. indzer(numzer)=(2*g-1)/2;
     
146. numzer=numzer+1;
     
147. end
     
148. 
     
149. if (x(g-1)>x(g))&(x(g)>x(g+1))
     
150. continue;
     
151. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
     
152. nummin=nummin+1;
     
153. elseif (x(g-1) continue;
     
154. elseif (x(g-1)x(g+1))
     
155. indmax(nummax)=g;
     
156. nummax=nummax+1;
     
157. end
     
158. end
     
159. %subplot(2,3,1)
     
160. plot(indzer,x(indzer),'r-+')
     
161. title('indzer')
     
162. pause
     
163. % subplot(2,3,2)
     
164. plot(indmin,x(indmin),'b-o')
     
165. title('indmin')
     
166. pause
     
167. %subplot(2,3,3)
     
168. plot(indmax,x(indmax),'g-*')
     
169. title('indmax')
     
170. pause
     
171. %subplot(2,3,4)
     
172. plot(1:length(x),x,'k-s')
     
173. title('x')
     
174. pause
     
175. 
     
176. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
     
177. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
     
178. %subplot(2,3,5)
     
179. plot(1:length(x),envmin,'y-d')
     
180. title('envmin')
     
181. pause
     
182. %subplot(2,3,6)
     
183. plot(1:length(x),envmax,'c-h')
     
184. title('envmax')
     
185. pause
     
186. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
     
187. title('mean((envmax+envmin)/2)')%extr
     
188. function [indmin,indmax,indzer]=extr(x)
     
189. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
     
190. watchon;
     
191. 
     
192. FILENAME=[DATApath,DATAfile];
     
193. x=load(FILENAME);
     
194. 
     
195. numzer=1;
     
196. nummin=1;
     
197. nummax=1;
     
198. for g=2:(length(x)-1)
     
199. 
     
200. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
     
201. indzer(numzer)=g-1;
     
202. numzer=numzer+1;
     
203. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
     
204. indzer(numzer)=g;
     
205. numzer=numzer+1;
     
206. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
     
207. indzer(numzer)=(2*g-1)/2;
     
208. numzer=numzer+1;
     
209. end
     
210. 
     
211. if (x(g-1)>x(g))&(x(g)>x(g+1))
     
212. continue;
     
213. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
     
214. nummin=nummin+1;
     
215. elseif (x(g-1) continue;
     
216. elseif (x(g-1)x(g+1))
     
217. indmax(nummax)=g;
     
218. nummax=nummax+1;
     
219. end
     
220. end
     
221. %subplot(2,3,1)
     
222. plot(indzer,x(indzer),'r-+')
     
223. title('indzer')
     
224. pause
     
225. % subplot(2,3,2)
     
226. plot(indmin,x(indmin),'b-o')
     
227. title('indmin')
     
228. pause
     
229. %subplot(2,3,3)
     
230. plot(indmax,x(indmax),'g-*')
     
231. title('indmax')
     
232. pause
     
233. %subplot(2,3,4)
     
234. plot(1:length(x),x,'k-s')
     
235. title('x')
     
236. pause
     
237. 
     
238. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
     
239. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
     
240. %subplot(2,3,5)
     
241. plot(1:length(x),envmin,'y-d')
     
242. title('envmin')
     
243. pause
     
244. %subplot(2,3,6)
     
245. plot(1:length(x),envmax,'c-h')
     
246. title('envmax')
     
247. pause
     
248. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
     
249. title('mean((envmax+envmin)/2)')%extr
     
250. function [indmin,indmax,indzer]=extr(x)
     
251. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
     
252. watchon;
     
253. 
     
254. FILENAME=[DATApath,DATAfile];
     
255. x=load(FILENAME);
     
256. 
     
257. numzer=1;
     
258. nummin=1;
     
259. nummax=1;
     
260. for g=2:(length(x)-1)
     
261. 
     
262. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
     
263. indzer(numzer)=g-1;
     
264. numzer=numzer+1;
     
265. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
     
266. indzer(numzer)=g;
     
267. numzer=numzer+1;
     
268. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
     
269. indzer(numzer)=(2*g-1)/2;
     
270. numzer=numzer+1;
     
271. end
     
272. 
     
273. if (x(g-1)>x(g))&(x(g)>x(g+1))
     
274. continue;
     
275. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
     
276. nummin=nummin+1;
     
277. elseif (x(g-1) continue;
     
278. elseif (x(g-1)x(g+1))
     
279. indmax(nummax)=g;
     
280. nummax=nummax+1;
     
281. end
     
282. end
     
283. %subplot(2,3,1)
     
284. plot(indzer,x(indzer),'r-+')
     
285. title('indzer')
     
286. pause
     
287. % subplot(2,3,2)
     
288. plot(indmin,x(indmin),'b-o')
     
289. title('indmin')
     
290. pause
     
291. %subplot(2,3,3)
     
292. plot(indmax,x(indmax),'g-*')
     
293. title('indmax')
     
294. pause
     
295. %subplot(2,3,4)
     
296. plot(1:length(x),x,'k-s')
     
297. title('x')
     
298. pause
     
299. 
     
300. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
     
301. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
     
302. %subplot(2,3,5)
     
303. plot(1:length(x),envmin,'y-d')
     
304. title('envmin')
     
305. pause
     
306. %subplot(2,3,6)
     
307. plot(1:length(x),envmax,'c-h')
     
308. title('envmax')
     
309. pause
     
310. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
     
311. title('mean((envmax+envmin)/2)')%extr
     
312. function [indmin,indmax,indzer]=extr(x)
     
313. [DATAfile DATApath]=uigetfile('*.txt','输入信号');
     
314. watchon;
     
315. 
     
316. FILENAME=[DATApath,DATAfile];
     
317. x=load(FILENAME);
     
318. 
     
319. numzer=1;
     
320. nummin=1;
     
321. nummax=1;
     
322. for g=2:(length(x)-1)
     
323. 
     
324. if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
     
325. indzer(numzer)=g-1;
     
326. numzer=numzer+1;
     
327. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
     
328. indzer(numzer)=g;
     
329. numzer=numzer+1;
     
330. elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
     
331. indzer(numzer)=(2*g-1)/2;
     
332. numzer=numzer+1;
     
333. end
     
334. 
     
335. if (x(g-1)>x(g))&(x(g)>x(g+1))
     
336. continue;
     
337. elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
     
338. nummin=nummin+1;
     
339. elseif (x(g-1) continue;
     
340. elseif (x(g-1)x(g+1))
     
341. indmax(nummax)=g;
     
342. nummax=nummax+1;
     
343. end
     
344. end
     
345. %subplot(2,3,1)
     
346. plot(indzer,x(indzer),'r-+')
     
347. title('indzer')
     
348. pause
     
349. % subplot(2,3,2)
     
350. plot(indmin,x(indmin),'b-o')
     
351. title('indmin')
     
352. pause
     
353. %subplot(2,3,3)
     
354. plot(indmax,x(indmax),'g-*')
     
355. title('indmax')
     
356. pause
     
357. %subplot(2,3,4)
     
358. plot(1:length(x),x,'k-s')
     
359. title('x')
     
360. pause
     
361. 
     
362. envmax = interp1(indmax,x(indmax),1:length(x),'spline');
     
363. envmin = interp1(indmin,x(indmin),1:length(x),'spline');
     
364. %subplot(2,3,5)
     
365. plot(1:length(x),envmin,'y-d')
     
366. title('envmin')
     
367. pause
     
368. %subplot(2,3,6)
     
369. plot(1:length(x),envmax,'c-h')
     
370. title('envmax')
     
371. pause
     
372. plot(1:length(x),mean((envmax+envmin)/2),'m-p')
     
373. title('mean((envmax+envmin)/2)')

复制代码

[ 本帖最后由 yejet 于 2006-9-17 12:58 编辑 ]

zhishan

本人东南大学无线电系大三学生,正准备做个SRTP,各位大侠多多指教.MEIL32的程序我收下了,缺的两个函数有时间发到syyu514@126.com.谢谢各位。

[ 本帖最后由 yejet 于 2006-9-17 12:58 编辑 ]

sealily

我刚刚下载
多谢了！