设为首页收藏本站联系我们
切换到宽版

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛»论坛 › 工程技术 › 信号处理方法 › [分享]个人收集的一些关于小波分析的matlab程序
1234下一页
返回列表 发新帖
查看: 55049|回复: 51

[小波] [分享]个人收集的一些关于小波分析的matlab程序

 关闭 [复制链接]
simon21
发表于 2006-3-31 08:34 | 显示全部楼层 |阅读模式

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要 登录 才可以下载或查看，没有账号？我要加入

x
都是从网上收集来的，由于时间比较久，处处都忘记了，如果是谁的原创请和我联系，我在帖子中标出来的
内容比较多，将会逐步贴出来

提升法97经典程序 （二楼）
2代小波示意程序 （三楼）
二代小波漫谈 （四楼）
小波滤波器构造和消噪程序（2个） （五楼）
小波谱分析mallat算法经典程序 （六楼）
2维小波变换经典程序 （七楼）
基于LeventCodes平台的小波去噪程序包 （十一楼）
连续小波和离散小波分析的应用实例（十二楼）
小波插值与小波构造（3个程序）（十三楼）
采用多孔trous算法(undecimated wavelet transform)实现小波变换（十四楼）
Daubechies小波基的构造（十五楼）
消失矩作用的程序（二十三楼）
平移变换平移法(cycle_spinning)消除gibbs效应 (二十四楼)
使用小波包变换分析信号的MATLAB程序（五十四楼）
基于小波消噪的雷达回波检测，可以检测雷达回波的有无及其准确的位置（五十五楼）
二维小波变换（正和逆变换）（五十六楼）
第二代小波变换源码（五十七楼）
利用小波变换实现对电能质量检测的算法实现（五十八楼）
基于小波变换的图象去噪 Normalshrink算法（五十九楼）
基于小波变换模极大的多尺度图像边缘检测（六十楼）
利用小波变根据二进制数（水印）来改变图片，提取其中一个子带的直方图（六十一楼）
用小波函数构建神经网络的源程序（六十二楼）
利用小波和霍夫曼对单声道文件进行压缩编码 并解码输出（六十三楼）
用小波神经网络来对时间序列进行预测（六十四楼）
基于小波特征提取方法的图象匹配算法（六十五楼）

今天（2007年4月4日）先贴到这里

[ 本帖最后由 simon21 于 2007-4-4 07:39 编辑 ]

本帖被以下淘专辑推荐:

回复
分享到：

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:36 | 显示全部楼层

提升法97经典程序


  2. %% 本程序实现任意偶数大小图像第二代双正交97提升小波变换 %% 注1： 采用标准正交方法，对行列采用不同矩阵（和matlab里不同）
  3. %% 注2： 为了保证正交，所有边界处理，全部采用循环处理
  4. %% 注3： 正交性验证，将单位阵带入函数，输出仍是单位阵（matlab不具有此性质）
  5. %% 注4： 此程序是矩阵实现，所以图像水平分量和垂直分量估计被交换位置
  6. %% 注5： 此程序实现的是类小波（wavelet-like）变换，是介于小波包变换与小波变换之间的变换
  7. %% 注6： 此程序每层变换相对原图像矩阵，产生的矩阵都是正交阵，这和小波包一致
  8. %% 注7： 但小波变换每层产生的矩阵，是相对每个待分解子块的正交矩阵，而不是原图像的正交矩阵
  9. %% 注8： 且小波变换产生的正交矩阵维数，随分解层数2分减少
  10. %% 注9: 提升系数可以在MATLAB7.0以上版本，用liftwave('9.7')获取，这里直接给出，考虑兼容性
  11. %% 注10：由于MATLAB数组下标从1开始，所以注意奇偶序列的变化
  12. %% 注11：d为对偶上升，即预测；p为原上升，即更新 %% 编程人 沙威 安徽大学
  13. %% 编程时间 2004年12月18日 %% x输入图像,y输出图像
  14. %% flag_trans为正变换或反变换标志，0执行正变换，1执行反变换
  15. %% flag_max,是否最大层数变换标志,0执行用户设定层数,1执行最大层数变换
  16. %% layer,用户层数设置（小于最大层） function y=db97(x,flag_trans,flag_max,layer); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  17. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  19. % 1.输入参数检查 % 矩阵维数判断
  20. [sa,sb]=size(x); if (sa~=sb) % 防止非图像数据
  21. errordlg('非图像数据！');
  22. error('非图像数据！');
  23. end; % 变换标志判断
  24. [sa,sb]=size(flag_trans);
  25. if ((sa~=1) | (sb~=1)) % 变换标志错误
  26. errordlg('变换标志错误！');
  27. error('变换标志错误！');
  28. end; if ((flag_trans~=1) & (flag_trans~=0)) % 变换标志错误
  29. errordlg('变换标志错误！');
  30. error('变换标志错误！');
  31. end; % 最大层数标志判断
  32. [sa,sb]=size(flag_max);
  33. if ((sa~=1) | (sb~=1)) % 最大层数标志错误
  34. errordlg('最大层数标志错误！');
  35. error('最大层数标志错误！');
  36. end; if ((flag_max~=1) & (flag_max~=0)) % 最大层数标志错误
  37. errordlg('最大层数标志错误！');
  38. error('最大层数标志错误！');
  39. end; % 用户设置层数判断
  40. if (flag_max~=1) [sa,sb]=size(layer);
  41. if ((sa~=1) | (sb~=1)) % 层数设置错误
  42. errordlg('层数设置错误！');
  43. error('层数设置错误！');
  44. end; if (flag_max<0) % 层数设置错误
  45. errordlg('层数设置错误！');
  46. error('层数设置错误！');
  47. end;
  48. end;
  49. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  50. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  52. % 2.提升系数确定
  53. % t1=liftwave('9.7'); % 获取提升系数(MATLAB7.0以后) d1=[-1.586100000000000e+000,-1.586134342069360e+000];
  54. p1=[1.079600000000000e+000,-5.298011857188560e-002];
  55. d2=[-8.829110755411875e-001,-8.829110755411875e-001];
  56. p2=[4.435068520511142e-001,1.576123746148364e+000];
  57. d3=-8.698644516247808e-001;
  58. p3=-1.149604398860242e+000;

  60. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  61. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  63. % 3.分解层数确定
  64. % 采用用户输入和自动给出最大层数两种方法 N=length(x); % 矩阵大小
  65. S=N; % 变量
  66. s=log2(N); % 最大循环次数
  67. n1=N/2; % 初始一半矩阵大小
  68. n2=N; % 初始矩阵大小
  69. u=0; % 初始值 % 对非2的整数幂大小图像确定最大分解层数
  70. for ss=1:s
  71. if (mod(S,2)==0)
  72. u=u+1;
  73. S=S/2;
  74. end;
  75. end;
  76. u=u-1; % 分解最大层数减1(后面的边界处理造成) % 最大层数确定
  77. if (flag_max==0) % 手动输入
  78. T=layer; % 用户输入值
  79. else % 自动确定最大层数
  80. T=u; % 分解最大层数
  81. end
  82. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  83. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  85. % 4.最大层数和图像大小检查 if (T>u) % 防止用户层数越界
  86. errordlg('已超过最大分解层数！或者非偶数大小图像！');
  87. error('已超过最大分解层数！或者非偶数大小图像！');
  88. end; if (mod(N,2)~=0) % 防止图像大小错误
  89. errordlg('非偶数大小图像！');
  90. error('非偶数大小图像！');
  91. end;
  92. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  93. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  95. % 5.提升法正变换 if (flag_trans==0)
  96. for time=1:T; % 行正变换

  98. % d;
  99. x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);
  100. x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:);

  102. % p;
  103. x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:);
  104. x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)*x1([2:n1],:);
  105. x([n1+1:n2],:)=x1([1:n1],:);

  107. % d;
  108. x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:);
  109. x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1)*x([2:n1],:);

  111. % p;
  112. x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:);
  113. x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:);
  114. x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1)*x([n1+3:n2],:);

  116. % 归一
  117. x([1:n1],:)=p3*x([1:n1],:);
  118. x([n1+1:n2],:)=d3*x([n1+1:n2],:); clear x1;

  120. % 列正变换

  122. % d;
  123. x1(:,[1:n1])=x(:,[2:2:n2]);

  125. % p;
  126. x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);
  127. x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1(:,[2:n1]);
  128. x(:,[n1+1:n2])=x1(:,[1:n1]);

  130. % d;
  131. x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);
  132. x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p1(2)*x(:,[2:n1]);

  134. % p;
  135. x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);
  136. x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2)*x(:,[n1+2:n2]);

  138. % d;
  139. x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1);
  140. x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);
  141. x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2)*x(:,[2:n1-1]);

  143. % 归一
  144. x(:,[1:n1])=d3*x(:,[1:n1]);
  145. x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]); clear x1;

  147. n2=n2/2; % 原大小
  148. n1=n2/2; % 一半大小
  149. end;
  150. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  151. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  153. % 6.提升法反变换 else
  154. n2=N/(2.^(T-1)); % 分解最小子块维数
  155. n1=n2/2;
  156. for time=1:T; % 行反变换

  158. % 去归一
  159. x([1:n1],:)=x([1:n1],:)/p3;
  160. x([n1+1:n2],:)=x([n1+1:n2],:)/d3; % 反p;
  161. x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:);
  162. x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);
  163. x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1)*x([n1+3:n2],:);

  165. % 反d;
  166. x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);
  167. x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1)*x([2:n1],:);

  169. % 反p;
  170. x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:);
  171. x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x([n1+2:n2],:);

  173. % 反d;
  174. x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);
  175. x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:);

  177. % 偶数
  178. x([1:2:n2-1],:)=x1([1:n1],:);

  180. clear x1;

  182. % 列反变换

  184. % 归一
  185. x(:,[1:n1])=x(:,[1:n1])/d3;
  186. x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3; % 反d;
  187. x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1);
  188. x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1);
  189. x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2)*x(:,[2:n1-1]);

  191. % 反p;
  192. x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1);
  193. x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2)*x(:,[n1+2:n2]);

  195. % 反d;
  196. x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);
  197. x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+p1(2)*x(:,[2:n1]);

  199. % 反p;
  200. x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);
  201. x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)*x(:,[n1+2:n2]); % 奇偶
  202. x(:,[2:2:n2])=x(:,[n1+1:n2]);
  203. x(:,[1:2:n2-1])=x1(:,[1:n1]); clear x1;

  205. n2=n2*2; % 原大小
  206. n1=n2/2; % 一半大小 end;
  207. end;

  209. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  210. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  212. % 7.结果输出 y=x;
  213. % 传输最后结果 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  214. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  216. % 8.内存清理 clear x;
  217. clear flag_max;
  218. clear layer;
  219. clear flag_trans;
  220. clear N;
  221. clear n1;
  222. clear n2;
  223. clear s;
  224. clear ss;
  225. clear u;
  226. clear d1;
  227. clear d2;
  228. clear d3;
  229. clear p1;
  230. clear p2;
  231. clear p3;
  232. clear sa;
  233. clear sb;

  235. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  236. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

[ 本帖最后由 yejet 于 2006-8-31 20:32 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:38 | 显示全部楼层

2代小波示意程序

  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %% 此程序用提升法实现第二代小波变换
  3. %% 我用的是非整数阶小波变换
  4. %% 采用时域实现,步骤先列后行
  5. %% 正变换：分裂，预测，更新；
  6. %% 反变换：更新，预测，合并
  7. %% 只做一层（可以多层，而且每层的预测和更新方程不同） clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  8. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  9. %%%%

  11. % 1.调原始图像矩阵 load wbarb; % 下载图像
  12. f=X; % 原始图像
  13. % f=[0 0 0 0 0 0 0 0 ;...
  14. % 0 0 0 1 1 0 0 0 ;...
  15. % 0 0 2 4 4 2 0 0 ;...
  16. % 0 1 4 8 8 4 1 0 ;...
  17. % 0 1 4 8 8 4 1 0 ;...
  18. % 0 0 2 4 4 2 0 0 ;...
  19. % 0 0 0 1 1 0 0 0 ;...
  20. % 0 0 0 0 0 0 0 0 ;]; % 原始图像矩阵 N=length(f); % 图像维数
  21. T=N/2;

  23. % 子图像维数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  24. %%%%%%%%%%%%%%%%%正变换%%%%%%%%%%%%%%%%%%%%%%%%%
  25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  26. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  27. %%%% 1.列变换

  29. % A.分裂（奇偶分开） f1=f([1:2:N-1],:); % 奇数
  30. f2=f([2:2:N],:); % 偶数 % f1(:,T+1)=f1(:,1); % 补列
  31. % f2(T+1,:)=f2(1,:); % 补行 % B.预测 for i_hc=1:T;
  32. high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
  33. end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 补行 % C.更新 for i_lc=1:T;
  34. low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
  35. end; % D.合并
  36. f_column([1:1:T],:)=low_frequency_column([1:T],:);
  37. f_column([T+1:1:N],:)=high_frequency_column([1:T],:);


  40. figure(1)
  41. colormap(map);
  42. image(f); figure(2)
  43. colormap(map);
  44. image(f_column);

  46. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  47. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  48. %%%% 2.行变换

  50. % A.分裂（奇偶分开） f1=f_column(:,[1:2:N-1]); % 奇数
  51. f2=f_column(:,[2:2:N]); % 偶数
  52. % f2(:,T+1)=f2(:,1); % 补行 % B.预测 for i_hr=1:T;
  53. high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
  54. end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行 % C.更新 for i_lr=1:T;
  55. low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
  56. end; % D.合并
  57. f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);
  58. f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);
  59. figure(3)
  60. colormap(map);
  61. image(f_row);
  62. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63. %%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%%
  64. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  65. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  66. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  68. %%%% 1.行变换
  69. % A.提取（低频高频分开） f1=f_row(:,[T+1:1:N]); % 奇数
  70. f2=f_row(:,[1:1:T]); % 偶数
  71. % f2(:,T+1)=f2(:,1); % 补行 % B.更新 for i_lr=1:T;
  72. low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
  73. end; % C.预测 for i_hr=1:T;
  74. high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
  75. end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行
  76. % D.合并(奇偶分开合并)
  77. f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);
  78. f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);
  79. figure(4)
  80. colormap(map);
  81. image(f_row);

  83. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  84. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  85. %%%% 2.列变换

  87. % A.提取（低频高频分开） f1=f_row([T+1:1:N],:); % 奇数
  88. f2=f_row([1:1:T],:); % 偶数 % f1(:,T+1)=f1(:,1); % 补列
  89. % f2(T+1,:)=f2(1,:); % 补行 % B.更新 for i_lc=1:T;
  90. low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
  91. end; % C.预测 for i_hc=1:T;
  92. high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
  93. end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 补行 % D.合并(奇偶分开合并）
  94. f_column([2:2:N],:)=low_frequency_column([1:T],:);
  95. f_column([1:2:N-1],:)=high_frequency_column([1:T],:);


  98. figure(5)
  99. colormap(map);
  100. image(f_column);
复制代码

[ 本帖最后由 yejet 于 2006-8-31 20:32 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:40 | 显示全部楼层

二代小波漫谈

现在我就举例，对一个8点序列，怎样实现第二代小波变换。

1. 奇偶分开。
非常简单，就是[2,4,6,8]组成一列向量，[1,3,5,7]组成一列向量。

2. 预测。
用[2,4,6,8]来预测[1,3,5,7]。比如说1，3估计2； 3，5估计4； 5，7估计6； 7，1估计8。（边缘处理，我采用循环方法）。估计公式可以用别人的，也可以自己做。举一个线性的例子：2=1*a+3*b,4=3*a+ 5*b,...，其他的都一样。这样我们就可找到最优的a,b，使得(2-(1*a+3*b)).^2+(4-(3*a+5*b)).^2+...最小化。就是最小均方准则。若正好为零，说明偶可以完全预测奇，也就是我们只要存储偶数列向量，和a,b就可以了，压缩也就是实现了。对于信号很长序列，就等于压缩了一半。当然，我们可以采用更复杂的立方差值预测，多项式预测，或其它的准则，来使其最小，这样我们的压缩也就得到了最优。

3. 提升。
我们总希望，均方为零，但可望不可及。于是，提升就需要了。我们经过预测后，要存储的是偶数序列[2,4,6,8]，新的奇数序列[n1,n3,n5, n7]=[2-(1*a+3*b),4-(3*a+5*b),...]和线性变换系数(a,b)。这里新的奇数序列就是高频分量。但偶数序列是不能完全代表信号的性质的，有所差距。所以我们要对偶数序列进行修正。即所谓的提升。我们这次用个简单的提升吧。[n2,n4,n6,n8]=[2,4,6,8]+ k*[n1,n3,n5,n7]。[n2,n4,n6,n8],就是要分解的低频分量。那k怎么求呢？因为要保持n2,n4,n6,n8和原始信号 [1,2,3,4,5,6,7,8]一样的性质。一般就是均值和高阶矩。这里就一个未知数k,所以用均值相等，就行了。1/8*(1+2+3+..8)= 1/4*(n2+n4+n6+n8)。k很容易就求出来了。我们最终存储的就是[n1,n3,n5,n7]和[n2,n4,n6,n8]以及a,b,k。

现在，所谓的第二代就完了。再说几句。
1.反变换，就是3->2->1。

2.二维。先行提升，再列提升。（我置顶的贴子里有harr二维提升的源代码）。

3.整数阶。就是加一个取整。

4.多层或小波包提升，就是在对序列[n1,n3,n5,n7]或[n2,n4,n6,n8]，再做1->2->3。

5.灵活。不一定是a,b,也可能就一个a,或a,b,c;不一定是一个k,也可能是k1,k2。但越多计算量太大。最好是用大师们做好的CDF,5/3,7/9等。

6.最重要的，任何一代小波，总可以通过一次或多次提升实现。它和一代小波没有本质区别。

7.优势。文献都有，我随便谈谈。时域实现，最优压缩，无边缘效应，灵活多变，无损压缩，编程方便，速度快。

文章写完了，希望对大家有帮助。最主要的，动手编，不要依赖MATLABM，这样才有所体会。希望和大家多交流。

给 simon21 加一点人气

[ 本帖最后由 simon21 于 2007-4-4 07:05 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:45 | 显示全部楼层

小波滤波器构造和消噪程序（2个）

1．重构

  1. % mallet_wavelet.m

  3. % 此函数用于研究Mallet算法及滤波器设计

  5. % 此函数仅用于消噪

  7. a=pi/8; %角度赋初值

  9. b=pi/8;

  11. %低通重构FIR滤波器h0(n)冲激响应赋值

  13. h0=cos(a)*cos(b);

  15. h1=sin(a)*cos(b);

  17. h2=-sin(a)*sin(b);

  19. h3=cos(a)*sin(b);

  21. low_construct=[h0,h1,h2,h3];

  23. L_fre=4; %滤波器长度

  25. low_decompose=low_construct(end:-1:1); %确定h0(-n)，低通分解滤波器

  27. for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器

  29. if(mod(i_high,2)==0);

  31. coefficient=-1;

  33. else

  35. coefficient=1;

  37. end

  39. high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;

  41. end

  43. high_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n)

  45. L_signal=100; %信号长度

  47. n=1:L_signal; %信号赋值

  49. f=10;

  51. t=0.001;

  53. y=10*cos(2*pi*50*n*t).*exp(-20*n*t);

  55. figure(1);

  57. plot(y);

  59. title('原信号');

  61. check1=sum(high_decompose); %h0(n)性质校验

  63. check2=sum(low_decompose);

  65. check3=norm(high_decompose);

  67. check4=norm(low_decompose);

  69. l_fre=conv(y,low_decompose); %卷积

  71. l_fre_down=dyaddown(l_fre); %抽取,得低频细节

  73. h_fre=conv(y,high_decompose);

  75. h_fre_down=dyaddown(h_fre); %信号高频细节

  77. figure(2);

  79. subplot(2,1,1)

  81. plot(l_fre_down);

  83. title('小波分解的低频系数');

  85. subplot(2,1,2);

  87. plot(h_fre_down);

  89. title('小波分解的高频系数');

  91. l_fre_pull=dyadup(l_fre_down); %0差值

  93. h_fre_pull=dyadup(h_fre_down);

  95. l_fre_denoise=conv(low_construct,l_fre_pull);

  97. h_fre_denoise=conv(high_construct,h_fre_pull);

  99. l_fre_keep=wkeep(l_fre_denoise,L_signal); %取结果的中心部分,消除卷积影响

  101. h_fre_keep=wkeep(h_fre_denoise,L_signal);

  103. sig_denoise=l_fre_keep+h_fre_keep; %信号重构

  105. compare=sig_denoise-y; %与原信号比较

  107. figure(3);

  109. subplot(3,1,1)

  111. plot(y);

  113. ylabel('y'); %原信号

  115. subplot(3,1,2);

  117. plot(sig_denoise);

  119. ylabel('sig\_denoise'); %重构信号

  121. subplot(3,1,3);

  123. plot(compare);

  125. ylabel('compare'); %原信号与消噪后信号的比较
复制代码



2．消噪

% mallet_wavelet.m

% 此函数用于研究Mallet算法及滤波器设计

% 此函数用于消噪处理

%角度赋值

%此处赋值使滤波器系数恰为db9

%分解的高频系数采用db9较好,即它的消失矩较大

%分解的有用信号小波高频系数基本趋于零

%对于噪声信号高频分解系数很大，便于阈值消噪处理

[l,h]=wfilters('db10','d');

low_construct=l;

L_fre=20; %滤波器长度

low_decompose=low_construct(end:-1:1); %确定h0(-n)，低通分解滤波器

for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器

if(mod(i_high,2)==0);

coefficient=-1;

else

coefficient=1;

end

high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;

end

high_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n)

L_signal=100; %信号长度

n=1:L_signal; %原始信号赋值

f=10;

t=0.001;

y=10*cos(2*pi*50*n*t).*exp(-30*n*t);

zero1=zeros(1,60); %信号加噪声信号产生

zero2=zeros(1,30);

noise=[zero1,3*(randn(1,10)-0.5),zero2];

y_noise=y+noise;

figure(1);

subplot(2,1,1);

plot(y);

title('原信号');

subplot(2,1,2);

plot(y_noise);

title('受噪声污染的信号');

check1=sum(high_decompose); %h0(n),性质校验

check2=sum(low_decompose);

check3=norm(high_decompose);

check4=norm(low_decompose);

l_fre=conv(y_noise,low_decompose); %卷积

l_fre_down=dyaddown(l_fre); %抽取,得低频细节

h_fre=conv(y_noise,high_decompose);

h_fre_down=dyaddown(h_fre); %信号高频细节

figure(2);

subplot(2,1,1)

plot(l_fre_down);

title('小波分解的低频系数');

subplot(2,1,2);

plot(h_fre_down);

title('小波分解的高频系数');

% 消噪处理

for i_decrease=31:44;

if abs(h_fre_down(1,i_decrease))>=0.000001

h_fre_down(1,i_decrease)=(10^-7);

end

end

l_fre_pull=dyadup(l_fre_down); %0差值

h_fre_pull=dyadup(h_fre_down);

l_fre_denoise=conv(low_construct,l_fre_pull);

h_fre_denoise=conv(high_construct,h_fre_pull);

l_fre_keep=wkeep(l_fre_denoise,L_signal); %取结果的中心部分,消除卷积影响

h_fre_keep=wkeep(h_fre_denoise,L_signal);

sig_denoise=l_fre_keep+h_fre_keep; %消噪后信号重构

%平滑处理

for j=1:2

for i=60:70;

sig_denoise(i)=sig_denoise(i-2)+sig_denoise(i+2)/2;

end;

end;

compare=sig_denoise-y; %与原信号比较

figure(3);

subplot(3,1,1)

plot(y);

ylabel('y'); %原信号

subplot(3,1,2);

plot(sig_denoise);

ylabel('sig\_denoise'); %消噪后信号

subplot(3,1,3);

plot(compare);

ylabel('compare'); %原信号与消噪后信号的比较

[ 本帖最后由 simon21 于 2007-4-4 07:04 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:46 | 显示全部楼层

小波谱分析mallat算法经典程序

  1. clc;clear;
  2. %% 1.正弦波定义
  3. f1=50; % 频率1
  4. f2=100; % 频率2
  5. fs=2*(f1+f2); % 采样频率
  6. Ts=1/fs; % 采样间隔
  7. N=120; % 采样点数
  8. n=1:N;
  9. y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合
  10. figure(1)
  11. plot(y);
  12. title('两个正弦信号')
  13. figure(2)
  14. stem(abs(fft(y)));
  15. title('两信号频谱')
  16. %% 2.小波滤波器谱分析
  17. h=wfilters('db30','l'); % 低通
  18. g=wfilters('db30','h'); % 高通
  19. h=[h,zeros(1,N-length(h))]; % 补零（圆周卷积，且增大分辨率变于观察）
  20. g=[g,zeros(1,N-length(g))]; % 补零（圆周卷积，且增大分辨率变于观察）
  21. figure(3);
  22. stem(abs(fft(h)));
  23. title('低通滤波器图')
  24. figure(4);
  25. stem(abs(fft(g)));
  26. title('高通滤波器图')
  27. %% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)
  28. sig1=ifft(fft(y).*fft(h)); % 低通(低频分量)
  29. sig2=ifft(fft(y).*fft(g)); % 高通(高频分量)
  30. figure(5); % 信号图
  31. subplot(2,1,1)
  32. plot(real(sig1));
  33. title('分解信号1')
  34. subplot(2,1,2)
  35. plot(real(sig2));
  36. title('分解信号2')
  37. figure(6); % 频谱图
  38. subplot(2,1,1)
  39. stem(abs(fft(sig1)));
  40. title('分解信号1频谱')
  41. subplot(2,1,2)
  42. stem(abs(fft(sig2)));
  43. title('分解信号2频谱')
  44. %% 4.MALLET重构算法
  45. sig1=dyaddown(sig1); % 2抽取
  46. sig2=dyaddown(sig2); % 2抽取
  47. sig1=dyadup(sig1); % 2插值
  48. sig2=dyadup(sig2); % 2插值
  49. sig1=sig1(1,[1:N]); % 去掉最后一个零
  50. sig2=sig2(1,[1:N]); % 去掉最后一个零
  51. hr=h(end:-1:1); % 重构低通
  52. gr=g(end:-1:1); % 重构高通
  53. hr=circshift(hr',1)'; % 位置调整圆周右移一位
  54. gr=circshift(gr',1)'; % 位置调整圆周右移一位
  55. sig1=ifft(fft(hr).*fft(sig1)); % 低频
  56. sig2=ifft(fft(gr).*fft(sig2)); % 高频
  57. sig=sig1+sig2; % 源信号
  58. %% 5.比较
  59. figure(7);
  60. subplot(2,1,1)
  61. plot(real(sig1));
  62. title('重构低频信号');
  63. subplot(2,1,2)
  64. plot(real(sig2));
  65. title('重构高频信号');
  66. figure(8);
  67. subplot(2,1,1)
  68. stem(abs(fft(sig1)));
  69. title('重构低频信号频谱');
  70. subplot(2,1,2)
  71. stem(abs(fft(sig2)));
  72. title('重构高频信号频谱');
  73. figure(9)
  74. plot(real(sig),'r','linewidth',2);
  75. hold on;
  76. plot(y);
  77. legend('重构信号','原始信号')
  78. title('重构信号与原始信号比较')
复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:04 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-3-31 08:46 | 显示全部楼层

2维小波变换经典程序


  2. % FWT_DB.M;
  3. % 此示意程序用DWT实现二维小波变换
  4. % 编程时间2004-4-10，编程人沙威
  5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6. %%%%%%%%%%%%%%%%%%%%%%%%%%%
  7. clear;clc;
  8. T=256; % 图像维数
  9. SUB_T=T/2; % 子图维数
  10. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  11. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  12. % 1.调原始图像矩阵
  13. load wbarb; % 下载图像
  14. f=X; % 原始图像
  15. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  16. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  17. % 2.进行二维小波分解
  18. l=wfilters('db10','l'); % db10（消失矩为10)低通分解滤波器冲击响应（长度为20）
  19. L=T-length(l);
  20. l_zeros=[l,zeros(1,L)]; % 矩阵行数与输入图像一致，为2的整数幂
  21. h=wfilters('db10','h'); % db10（消失矩为10)高通分解滤波器冲击响应（长度为20）
  22. h_zeros=[h,zeros(1,L)]; % 矩阵行数与输入图像一致，为2的整数幂
  23. for i=1:T; % 列变换
  24. row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积<->FFT
  25. row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积<->FFT
  26. end;
  27. for j=1:T; % 行变换
  28. line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ); % 圆周卷积<->FFT
  29. line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ); % 圆周卷积<->FFT
  30. end;
  31. decompose_pic=line; % 分解矩阵
  32. % 图像分为四块
  33. lt_pic=decompose_pic(1:SUB_T,1:SUB_T); % 在矩阵左上方为低频分量--fi(x)*fi(y)
  34. rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T); % 矩阵右上为--fi(x)*psi(y)
  35. lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T); % 矩阵左下为--psi(x)*fi(y)
  36. rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T); % 右下方为高频分量--psi(x)*psi(y)

  38. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  39. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  40. % 3.分解结果显示
  41. figure(1);
  42. colormap(map);
  43. subplot(2,1,1);
  44. image(f); % 原始图像
  45. title('original pic');
  46. subplot(2,1,2);
  47. image(abs(decompose_pic)); % 分解后图像
  48. title('decomposed pic');
  49. figure(2);
  50. colormap(map);
  51. subplot(2,2,1);
  52. image(abs(lt_pic)); % 左上方为低频分量--fi(x)*fi(y)
  53. title('\Phi(x)*\Phi(y)');
  54. subplot(2,2,2);
  55. image(abs(rt_pic)); % 矩阵右上为--fi(x)*psi(y)
  56. title('\Phi(x)*\Psi(y)');
  57. subplot(2,2,3);
  58. image(abs(lb_pic)); % 矩阵左下为--psi(x)*fi(y)
  59. title('\Psi(x)*\Phi(y)');
  60. subplot(2,2,4);
  61. image(abs(rb_pic)); % 右下方为高频分量--psi(x)*psi(y)
  62. title('\Psi(x)*\Psi(y)');


  65. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  66. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  67. % 5.重构源图像及结果显示
  68. % construct_pic=decompose_matrix'*decompose_pic*decompose_matrix;
  69. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  70. l_re=l_zeros(end:-1:1); % 重构低通滤波
  71. l_r=circshift(l_re',1)'; % 位置调整
  72. h_re=h_zeros(end:-1:1); % 重构高通滤波
  73. h_r=circshift(h_re',1)'; % 位置调整

  75. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  76. top_pic=[lt_pic,rt_pic]; % 图像上半部分
  77. t=0;
  78. for i=1:T; % 行插值低频

  80. if (mod(i,2)==0)
  81. topll(i,:)=top_pic(t,:); % 偶数行保持
  82. else
  83. t=t+1;
  84. topll(i,:)=zeros(1,T); % 奇数行为零
  85. end
  86. end;
  87. for i=1:T; % 列变换
  88. topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圆周卷积<->FFT
  89. end;

  91. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  92. bottom_pic=[lb_pic,rb_pic]; % 图像下半部分
  93. t=0;
  94. for i=1:T; % 行插值高频
  95. if (mod(i,2)==0)
  96. bottomlh(i,:)=bottom_pic(t,:); % 偶数行保持
  97. else
  98. bottomlh(i,:)=zeros(1,T); % 奇数行为零
  99. t=t+1;
  100. end
  101. end;
  102. for i=1:T; % 列变换
  103. bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圆周卷积<->FFT
  104. end;

  106. construct1=bottomch_re+topcl_re; % 列变换重构完毕

  108. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  109. left_pic=construct1(:,1:SUB_T); % 图像左半部分
  110. t=0;
  111. for i=1:T; % 列插值低频

  113. if (mod(i,2)==0)
  114. leftll(:,i)=left_pic(:,t); % 偶数列保持
  115. else
  116. t=t+1;
  117. leftll(:,i)=zeros(T,1); % 奇数列为零
  118. end
  119. end;
  120. for i=1:T; % 行变换
  121. leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圆周卷积<->FFT
  122. end;


  125. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  126. right_pic=construct1(:,SUB_T+1:T); % 图像右半部分

  128. t=0;
  129. for i=1:T; % 列插值高频
  130. if (mod(i,2)==0)
  131. rightlh(:,i)=right_pic(:,t); % 偶数列保持
  132. else
  133. rightlh(:,i)=zeros(T,1); % 奇数列为零
  134. t=t+1;
  135. end
  136. end;
  137. for i=1:T; % 行变换
  138. rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圆周卷积<->FFT
  139. end;

  141. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  142. construct_pic=rightch_re+leftcl_re; % 重建全部图像

  144. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  145. % 结果显示
  146. figure(3);
  147. colormap(map);
  148. subplot(2,1,1);
  149. image(f); % 源图像显示
  150. title('original pic');

  152. subplot(2,1,2);
  153. image(abs(construct_pic)); % 重构源图像显示
  154. title('reconstructed pic');
  155. error=abs(construct_pic-f); % 重构图形与原始图像误值
  156. figure(4);
  157. mesh(error); % 误差三维图像
  158. title('absolute error display');
复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:05 编辑 ]
回复 支持 反对

使用道具 举报

eros
发表于 2006-4-2 09:40 | 显示全部楼层
请问能不能指教一下二次样条小波呢？matlab里面如何实现呢？
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-4-6 14:26 | 显示全部楼层

回复：（simon21）[分享]个人收集的一些关于小波分析的...

基于LeventCodes平台的小波去噪程序

包括以下方法：

BivaShrink方法、模型1、模型2、模型3（TrivaShrink方法）、BayesShrink方法、 LAWMLShrink方法的DWT实现和DT_CWT实现。
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-4-6 14:30 | 显示全部楼层

回复：（simon21）[分享]个人收集的一些关于小波分析的...

连续小波和离散小波分析的应用实例
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-4-6 14:31 | 显示全部楼层

小波插值与小波构造（3个程序）

  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. % 小波构造 function casade
  3. clear;clc;
  4. t=3;
  5. phi=[0,1,0]; h=wfilters('db7','r');
  6. h=h*sqrt(2); h_e=h(1,[2:2:14]);
  7. h_o=h(1,[1:2:13]); for m=1:15;
  8. stem(phi);
  9. drawnow;
  10. pause(1);
  11. ee=conv(h_e,phi);
  12. oo=conv(h_o,phi);
  13. phi(1,[2:2:2*length(ee)])=ee;
  14. phi(1,[1:2:2*length(oo)-1])=oo;


  17. end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  18. % cubic_average（立方b样条）
  19. % 均值插值 % 初始化
  20. s=[0 0 1 0 0] % 正弦波
  21. % f=50;
  22. % ts=1/200;
  23. % n=0:16;
  24. % s=sin(2*pi*f*n*ts); % 系数
  25. se=[1/8,6/8,1/8];
  26. so=[4/8,4/8] % 循环
  27. for p=1:10;
  28. t=length(s)-1;
  29. o(1:t)=s(1:t)*so(1)+s(2:t+1)*so(2);
  30. e(1)=s(t+1)*se(1)+s(1)*se(2)+s(2)*se(3);
  31. e(2:t)=s(1:t-1)*se(1)+s(2:t)*se(2)+s(3:t+1)*se(3);
  32. e(t+1)=s(t)*se(1)+s(t+1)*se(2)+s(1)*se(3);
  33. s([1:2:2*t+1])=e([1:t+1]);
  34. s([2:2:2*t])=o([1:t]);
  35. plot(s);
  36. drawnow;
  37. end; % 抽取
  38. t=length(s); % 总长度
  39. p=128; % 需要点数 % 间隔
  40. d=(t-1)/p; % 最终尺度函数
  41. r=s(2:d:t-1); % 画图
  42. figure(2);
  43. plot(r);
  44. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  45. % cubic_subdivision（立方插值）
  46. % 细分插值 % %% 初始化（尺度函数）
  47. % s=[0,0,1,0,0]; % 正弦函数
  48. n=1:20;
  49. f=50;
  50. ts=1/200;
  51. s=sin(2*pi*f*n*ts); % 指数函数
  52. % n=0:16;
  53. % s=exp(n); % % 系数
  54. a=[-1/16,9/16,9/16,-1/16]; % 循环
  55. for p=1:4;
  56. t=length(s)-1;
  57. o(1)=s(4)*a(1)+s(1)*a(2)+s(2)*a(3)+s(3)*a(4);
  58. o(2:t-1)=s(1:t-2)*a(1)+s(2:t-1)*a(2)+s(3:t)*a(3)+s(4:t+1)*a(4);
  59. o(t)=s(t-2)*a(4)+s(t+1)*a(3)+s(t)*a(2)+s(t-1)*a(1);
  60. s([1:2:2*t+1])=s([1:t+1]);
  61. s([2:2:2*t])=o([1:t]);
  62. plot(s);
  63. drawnow;
  64. end; % % 抽取
  65. % t=length(s); % 总长度
  66. % p=128; % 需要点数
  67. %
  68. % % 间隔
  69. % d=(t-1)/p;
  70. %
  71. % % 最终尺度函数
  72. % r=s(2:d:t-1);
  73. %
  74. % % 画图
  75. % figure(2);
  76. % plot(r);
复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:06 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-4-6 14:33 | 显示全部楼层

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

  1. clear;clc;
  2. % 1.生成信号
  3. f=50; % 频率
  4. fs=800; % 采样率
  5. T=128; % 信号长度
  6. n=1:T;
  7. y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 信号
  8. % y=circshift(y.',3).'; %%
  9. 2.正变换 l1=wfilters('db4','l')*sqrt(2)/2;
  10. % 参考低通滤波器
  11. l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通滤波器1
  12. h1=wfilters('db4','h')*sqrt(2)/2; % 参考高通滤波器
  13. h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通滤波器1 low1=ifft(fft(y).*fft(l1_zeros)); % 低频分量1
  14. high1=ifft(fft(y).*fft(h1_zeros)); % 高频分量1 l2=dyadup(l1); % 原滤波器插值
  15. l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通滤波器2
  16. h2=dyadup(h1); % 原滤波器插值
  17. h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通滤波器2 low2=ifft(fft(low1).*fft(l2_zeros)); % 低频分量2
  18. high2=ifft(fft(low1).*fft(h2_zeros)); % 高频分量2
  19. %% 3.反变换 lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重构低通滤波器2
  20. hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重构高通滤波器2
  21. lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重构低通滤波器1
  22. hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重构高通滤波器1 lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重构低频分量1(lowr=low1)
  23. r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重构源信号(r_s=y)
  24. %% 4.绘图 figure(1);
  25. plot(y);
  26. title('源信号'); figure(2);
  27. plot(low1,'r');
  28. hold on;
  29. plot(low2,'b');
  30. legend('第一层低频','第二层低频'); figure(3);
  31. plot(high1,'r');
  32. hold on;
  33. plot(high2,'b');
  34. legend('第一层高频','第二层高频'); figure(4);
  35. plot(low1,'r');
  36. hold on;
  37. plot(lowr,'b.');
  38. legend('第一层低频','重构第一层低频'); figure(5);
  39. plot(y,'r');
  40. hold on;
  41. plot(r_s,'b.');
  42. legend('源信号','重构信号'); disp(norm(low1-lowr))
  43. disp(norm(y-r_s))
复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:07 编辑 ]
回复 支持 反对

使用道具 举报

simon21
 楼主| 发表于 2006-4-6 14:35 | 显示全部楼层

此程序实现构造小波基

  1. % 此程序实现构造小波基
  2. % periodic_wavelet.m function ss=periodic_wavelet; clear;clc; % global MOMENT; % 消失矩阶数
  3. % global LEFT_SCALET; % 尺度函数左支撑区间
  4. % global RIGHT_SCALET; % 尺度函数右支撑区间
  5. % global LEFT_BASIS; % 小波基函数左支撑区间
  6. % global RIGHT_BASIS; % 小波基函数右支撑区间
  7. % global MIN_STEP; % 最小离散步长
  8. % global LEVEL; % 计算需要的层数(离散精度）
  9. % global MAX_LEVEL; % 周期小波最大计算层数
  10. [s2,h]=scale_integer;
  11. [test,h]=scalet_stretch(s2,h);
  12. wave_base=wavelet(test,h);
  13. ss=periodic_waveletbasis(wave_base);
  14. function [s2,h]=scale_integer; % 本函数实现求解小波尺度函数离散整数点的值
  15. % sacle_integer.m MOMENT=10; % 消失矩阶数
  16. LEFT_SCALET=0; % 尺度函数左支撑区间
  17. RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
  18. LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
  19. RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
  20. MIN_STEP=1/512; % 最小离散步长
  21. LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
  22. MAX_LEVEL=8; % 周期小波最大计算层数
  23. h=wfilters('db10','r'); % 滤波器系数 h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1; for i=LEFT_SCALET+1:RIGHT_SCALET-1
  24. for j=LEFT_SCALET+1:RIGHT_SCALET-1
  25. k=2*i-j+1;
  26. if (k>=1&k<=RIGHT_SCALET+1)
  27. a(i,j)=h(k); % 矩阵系数矩阵
  28. else
  29. a(i,j)=0;
  30. end
  31. end
  32. end [s,w]=eig(a); % 求特征向量，解的基
  33. s1=s(:,1);
  34. s2=[0;s1/sum(s1);0]; % 根据条件SUM(FI(T))=1,求解; % 本函数实现尺度函数经伸缩后的离散值
  35. % scalet_stretch.m function [s2,h]=scalet_stretch(s2,h); MOMENT=10; % 消失矩阶数
  36. LEFT_SCALET=0; % 尺度函数左支撑区间
  37. RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
  38. LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
  39. RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
  40. MIN_STEP=1/512; % 最小离散步长
  41. LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
  42. MAX_LEVEL=8; % 周期小波最大计算层数
  43. for j=1:LEVEL % 需要计算到尺度函数的层数
  44. t=0;
  45. for i=1:2:2*length(s2)-3 % 需要计算的离散点取值（0，1，2，3 -> 1/2, 3/2, 5/2)
  46. t=t+1;
  47. fi(t)=0;
  48. for n=LEFT_SCALET:RIGHT_SCALET; % 低通滤波器冲击响应紧支撑判断
  49. if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) % 小波尺度函数紧支撑判断
  50. fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1); % 反复应用双尺度方程求解
  51. end
  52. end
  53. end
  54. clear s
  55. n1=length(s2);
  56. n2=length(fi);
  57. for i=1:length(s2)+length(fi) % 变换后的矩阵长度
  58. if (mod(i,2)==1)
  59. s(i)=s2((i+1)/2); % 矩阵奇数下标为小波上一层(0,1,2,3)离散值
  60. else
  61. s(i)=fi(i/2); % 矩阵偶数下标为小波下一层(1/2,3/2,5/2)(经过伸缩变换后)的离散值
  62. end
  63. end
  64. s2=s;
  65. end
  66. % 采用双尺度方程求解小波基函数 PSI(T)
  67. % wavelet.m function wave_base=wavelet(test,h); MOMENT=10; % 消失矩阶数
  68. LEFT_SCALET=0; % 尺度函数左支撑区间
  69. RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
  70. LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
  71. RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
  72. MIN_STEP=1/512; % 最小离散步长
  73. LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
  74. MAX_LEVEL=8; % 周期小波最大计算层数
  75. i=0;
  76. for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS; % 小波基支撑长度
  77. s=0;
  78. for n=1-RIGHT_SCALET:1-LEFT_SCALET % g(n)取值范围
  79. if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET) % 尺度函数判断
  80. s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1); % 计算任意精度的小波基函数值
  81. end
  82. end
  83. i=i+1;
  84. wave_base(i)=s;
  85. end
复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:07 编辑 ]
回复 支持 反对

使用道具 举报

infant
发表于 2006-4-26 16:18 | 显示全部楼层
本帖最后由 wdhd 于 2016-3-17 14:21 编辑

　　请问楼主有对称小波的程序吗?谢谢了

　　[em04]
回复 支持 反对

使用道具 举报

ruola
发表于 2006-4-30 10:50 | 显示全部楼层
请问有没有小波系数模极大值的程序啊！
回复 支持 反对

使用道具 举报

下一页 »
1234下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-28 19:12 , Processed in 0.149622 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表