【日光倾城整理】小波变换中的能量泄漏问题的研究

日光傾城

【小波分析】 db10的 6层分解 重构的小波能量谱分布

一、案例分析

采用 db10的 6层小波分析，经过分解、重构观察小波能量谱分布情概况。

分析的信号频率为16.7HZ。采样频率1000HZ，采样点N=1024。

由于小波分解遵循Nyquist标准，因此最大可测频率为采样频率的一半，即500HZ
对信号进行六层分解,由Mallat多分辨率算法， 频带划分：


d1:500--250
d2:250--125
d3:125--62.5
d4:62.5--31.25
d5:31.25--15.625

d6:15.625-7.8125

a6:7.8125--DC

出现问题：1）16.7HZ信号应该落在d5中，但是，能量直方图显示的a6重构的成分更多，即DC--7.8125之间。
                  2）原信号改成50HZ和200HZ 靠近频段的中心时，就不会发生这种现象了。


问题原因：小波变换中的能量泄漏问题


解决办法：改变采样频率（或者分解层数N），使感兴趣频率靠近频带中心位置。



未解决问题：由于谐波的存在，16.7HZ的信号会被分解1/2，1,2倍频。基频16.7的幅值比实际幅值有所衰减。

这样很难判断风机转速1000r/min，即频率16.7HZ，

当出现y=a*sin(2*pi*(1/2)*16.7)+b*sin(2*pi*16.7)+c*sin(2*pi*2*16.7)+...+f*sin(2*pi*5*16.7)的叠加情况时，实际的幅值会与小波分析得到的重构时间序列的幅值有所差别。期待大家能解决这个问题。

同样的问题也会出现FFT变换中。

例如y=sin(x),用离散点来描述时，1）离散的点越多，波形趋势越清晰。即采样点越多，频谱泄露越少。

                              2）采样频率，必须大于待测信号的2倍（最低标准）。

t=(0,2pi)时间段内，采样周期T=1，采样频率F=1，可以采集6点。采样频率F=2，可以采集12点。

因此，每个周期内采集20个点计算，采样频率选择20*（1/2pi）=3.2HZ.过大的采样长度会影响计算的速度。

日光傾城

仿真图：
图一 原始波形图：



图二小波分解系数：




图三 重构的新时间序列：


图四 能量直方图：

日光傾城

程序代码：

1. clear all;clc;
   
2. f1=16.7;%  频率1
   
3. f2=33.4;
   
4. f3=50.1;
   
5. Fs=1000;   %  采样频率
   
6. Ts=1/Fs;   %  采样间隔
   
7. N=1024;    %  采样点数
   
8. t=[0:1/Fs:(N-1)/Fs];   %采样时刻
   
9. lev=5;
   
10. x=sin(2*pi*f1*t);
    
11. figure(1);
    
12. plot(t,x);
    
13. %%小波分解
    
14. [C,L]=wavedec(x,lev,'db10');   %C由[cAj,cDj,cDj-1,...,cD1]     
    
15.   %提取多尺度小波变换的高频系数
    
16. figure(2);
    
17. lev_1=lev+1;
    
18. for i=1:lev
    
19.             cD=detcoef(C,L,i);
    
20.         subplot(lev_1,1,i)
    
21.         plot(cD);
    
22.             ylabel(['cD',num2str(i)]);
    
23.         title(['Detail cD',num2str(i)]);
    
24. end
    
25.   %提取多尺度小波变换的低频系数
    
26. figure(2);
    
27.    
    
28.        cA=appcoef(C,L,'db10',lev);
    
29.        subplot(lev_1,1,lev_1)
    
30.        plot( cA);
    
31.        ylabel(['cA',num2str(lev)]);
    
32.        title(['Approximation cA',num2str(lev)])
    
33. 
    
34. %%小波重构
    
35. figure(3);
    
36. for i=1:lev
    
37.             D=wrcoef('d',C,L,'db10',i);
    
38.         subplot(lev_1,1,i)
    
39.         plot(D);
    
40.             ylabel(['D',num2str(i)]);
    
41.         title(['Detail D',num2str(i)]);
    
42. end
    
43. 
    
44. figure(3);
    
45.    
    
46.        A=wrcoef('a',C,L,'db10',lev);
    
47.        subplot(lev_1,1,lev_1)
    
48.        plot( A);
    
49.        ylabel(['A',num2str(lev)]);
    
50.        title(['Approximation A',num2str(lev)]);
    
51.       
    
52. %%计算能量谱
    
53.    figure(4);
    
54. [Ea,Ed] = wenergy(C,L);
    
55. EE1=sum(Ed)+Ea;
    
56. for i=1:lev
    
57.       Ed1(lev_1-i)= Ed(i);
    
58. end
    
59. 
    
60. E1=[Ea Ed1]/EE1;
    
61. x2=1:lev_1;   %p=rand(6,1);
    
62. bar(x2,E1);
    
63. 
    
64. str = cell(1,lev_1); %建立单元数组1*10的空数组
    
65. for i = 0:lev
    
66.     str{i+1} = ['E',num2str(i)]; %num2str(i)数字转为字符串输出。
    
67. end
    
68. set(gca,'xticklabel',str)  

复制代码

日光傾城

【日光倾城整理】小波变换中的能量泄漏问题
离散小波分解通过Mallat算法来实现。Cn(j+1)=sum(H)C(j);Dn(j+1)=sum(G)C(j)
h(n),g(n)为正交镜像滤波器，由于h(n),g(n)始终不是理想的滤波器，在某一频率经过h(n),g(n)后，会像其他频率泄露。这就解释了信号中出现了原来并不含有的频率成份而使结果难以解释。
在频带中心位置，信号衰减的慢；靠近截止频率时，信号衰减很快。
选择小波分析时，要考虑待分析信号的频率分布，将感兴趣的信号尽量靠近各自的频带中心。
同时小波基的选择也会影响能量泄露，db4和db20比较，db20边缘泄露要小。具体具体可以通过改变小波基逐一比较。

只有真正的理解Mallat多分辨率的思想才能真正的理解小波变换。
只有理解了小波基函数的构成，才能选择适合的频带分析方式。

“通过MRA构造正交小波的本质就在于把整个L2空间正交分解，分解的依据就是MRA，然后再在每个子空间中构造正交基。有了MRA，构造正交基的过程就大大简化了，而且更重要的是保证有一个快速的算法来分解和重构函数。”
摘自百思论坛--- 谈谈多尺度分析的理论

改进 平方可积空间L^2(R)中的正交小波基，可以更好的解决实际问题。

日光傾城

忙了一个星期，终于解决一部分问题。下一阶段任务，小波滤波和奇异点检测。

woshiqiao

感觉楼主写的很有技术啊。。想问一下楼主小波分解时为何选用DB10？是试出来的吗？还是你只实验了DB小波。其他小波是怎么排除的？