小波分解A'Trous算法Matlab代码示例

zsh1975

　　采用多孔trous算法(undecimated wavelet transform)实现小波变换

1. 　　clear;clc;
   
2. 
   
3. 　　%% 1.生成信号
   
4. 
   
5. 　　f=50; % 频率
   
6. 
   
7. 　　fs=800; % 采样率
   
8. 
   
9. 　　T=128; % 信号长度
   
10. 
    
11. 　　n=1:T;
    
12. 
    
13. 　　%y=yuanshishuju
    
14. 
    
15. 　　y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 信号
    
16. 
    
17. 　　% y=circshift(y.',3).';
    
18. 
    
19. 　　%% 2.正变换
    
20. 
    
21. 　　l1=wfilters('db4','l')*sqrt(2)/2; % 参考低通滤波器
    
22. 
    
23. 　　l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通滤波器1
    
24. 
    
25. 　　h1=wfilters('db4','h')*sqrt(2)/2; % 参考高通滤波器
    
26. 
    
27. 　　h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通滤波器1
    
28. 
    
29. 　　low1=ifft(fft(y).*fft(l1_zeros)); % 低频分量1
    
30. 
    
31. 　　high1=ifft(fft(y).*fft(h1_zeros)); % 高频分量1
    
32. 
    
33. 　　l2=dyadup(l1); % 原滤波器插值
    
34. 
    
35. 　　l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通滤波器2
    
36. 
    
37. 　　h2=dyadup(h1); % 原滤波器插值
    
38. 
    
39. 　　h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通滤波器2
    
40. 
    
41. 　　low2=ifft(fft(low1).*fft(l2_zeros)); % 低频分量2
    
42. 
    
43. 　　high2=ifft(fft(low1).*fft(h2_zeros)); % 高频分量2
    
44. 
    
45. 　　%% 3.反变换
    
46. 
    
47. 　　lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重构低通滤波器2
    
48. 
    
49. 　　hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重构高通滤波器2
    
50. 
    
51. 　　lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重构低通滤波器1
    
52. 
    
53. 　　hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重构高通滤波器1
    
54. 
    
55. 　　lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重构低频分量1(lowr=low1)
    
56. 
    
57. 　　r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重构源信号(r_s=y)
    
58. 
    
59. 　　%% 4.绘图
    
60. 
    
61. 　　figure(1);
    
62. 
    
63. 　　plot(y);
    
64. 
    
65. 　　title('源信号');
    
66. 
    
67. 　　figure(2);
    
68. 
    
69. 　　plot(low1,'r');
    
70. 
    
71. 　　hold on;
    
72. 
    
73. 　　plot(low2,'b');
    
74. 
    
75. 　　legend('第一层低频','第二层低频');
    
76. 
    
77. 　　figure(3);
    
78. 
    
79. 　　plot(high1,'r');
    
80. 
    
81. 　　hold on;
    
82. 
    
83. 　　plot(high2,'b');
    
84. 
    
85. 　　legend('第一层高频','第二层高频');
    
86. 
    
87. 　　figure(4);
    
88. 
    
89. 　　plot(low1,'r');
    
90. 
    
91. 　　hold on;
    
92. 
    
93. 　　plot(lowr,'b.');
    
94. 
    
95. 　　legend('第一层低频','重构第一层低频');
    
96. 
    
97. 　　figure(5);
    
98. 
    
99. 　　plot(y,'r');
    
100. 
     
101. 　　hold on;
     
102. 
     
103. 　　plot(r_s,'b.');
     
104. 
     
105. 　　legend('源信号','重构信号');
     
106. 
     
107. 　　disp(norm(low1-lowr))
     
108. 
     
109. 　　disp(norm(y-r_s))
     

复制代码

转自：http://blog.sina.com.cn/s/blog_4d7c97a00101a0dd.html