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

simon21

二维小波变换（正和逆变换），用二维卷积实现的，其中有个子函数用于产生小波系数的（小波分解和小波重构的系数都有）！你自己也可以往里面加系数！

1. function [varargout]=wavefilter(wname,type)
   
2. error(nargchk(1,2,nargin));
   
3. if(nargin==1 & nargout ~=4) | (nargin==2 & nargout~=2)
   
4.     error('Invalid number of output arguments.');
   
5. end
   
6. if nargin==1 & ~ischar(wname)
   
7.     error('WNAME must be a string.');
   
8. end
   
9. if nargin==2 & ~ischar(type)
   
10.     error('TYPE must be a string.');
    
11. end
    
12. 
    
13. switch lower(wname)
    
14.     case {'haar','db1'}
    
15.         ld=[1 1]/sqrt(2); hd=[-1 1]/sqrt(2);
    
16.         lr=ld;            hr=-hd;
    
17.     case 'db4'
    
18.         ld=[-1.059740178499728e-002 3.288301166698295e-002 3.084138183598697e-002 -1.870348117188811e-001 ...
    
19.             -2.798376941698385e-002 6.308807679295904e-001 7.148465705525415e-001 2.303778133088552e-001];
    
20.         t=[0:7];
    
21.         hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
    
22.         lr=ld; lr(end:-1:1)=ld;
    
23.         hr=cos(pi*t).*ld;
    
24.         
    
25.     case 'sym4'
    
26.         ld=[-7.576571478927333e-002 -2.963552764599851e-002 4.976186676320155e-001 8.037387518059161e-001 ...
    
27.             2.978577956052774e-001  -9.921954357684722e-002 -1.260396726203783e-002 3.222310060404270e-002];
    
28.         t=[0:7];
    
29.         hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
    
30.         lr=ld; lr(end:-1:1)=ld;
    
31.         hr=cos(pi*t).*ld;
    
32.     case 'bior6.8'
    
33.         ld=[0 1.908831736481291e-003 -1.9142861290088767e-003 -1.699063986760234e-002 1.1934565279772926e-002 ...
    
34.             4.973290349094079e-002 -7.726317316720414e-002 -9.405920349573646e-002 4.207962846098268e-001 ...
    
35.             8.259229974584023e-001 4.207962846098268e-001 -9.405920349573646e-002 -7.726317316720414e-002 ...
    
36.             4.973290349094079e-002 1.193456527972926e-002 -1.699063986760234e-002 -1.914286129088767e-003 ...
    
37.             1.908831736481291e-003];
    
38.         hd=[0 0 0 1.442628250562444e-002 -1.446750489679015e-002 -7.872200106262882e-002 4.036797903033992e-002 ...
    
39.             4.178491091502746e-001 -7.589077294536542e-001 4.178491091502746e-001 4.036797903033992e-002 ...
    
40.             -7.872200106262882e-002 -1.446750489679015e-002 1.442628250562444e-002 0 0 0 0];
    
41.         t=[0:17];
    
42.         lr=cos(pi*(t+1)).*hd;
    
43.         hr=cos(pi*t).*ld;
    
44.         
    
45.     case 'jpeg9.7'
    
46.         ld=[0 0.02674875741080976 -0.01686411844287495 -0.07822326652898785 0.2668641184428723 ...
    
47.             0.6029490182363579 0.2668641184428723 -0.07822326652898785 -0.01686411844287495  ...
    
48.             0.02674875741080976];
    
49.         hd=[0 -0.09127176311424948 0.05754352622849957 0.5912717631142470 -1.115087052456994 ...
    
50.             0.5912717631142470 0.05754352622849957 -0.09127176311424948 0 0];
    
51.         t=[0:9];
    
52.         lr=cos(pi*(t+1)).*hd;
    
53.         hr=cos(pi*t).*ld;
    
54.     otherwise
    
55.         error('Unrecongizable wavelet name (WNAME).');
    
56. end
    
57. if(nargin==1)
    
58.     varargout(1:4)={ld,hd,lr,hr};
    
59. else
    
60.     switch lower(type(1))
    
61.         case 'd'
    
62.             varargout={ld,hd};
    
63.         case 'r'
    
64.             varargout={lr,hr};
    
65.         otherwise
    
66.             error('Unrecongizable filter TYPE.');
    
67.     end
    
68. end

复制代码

1. function [c,s]=wavefast(x,n,varargin)
   
2. error(nargchk(3,4,nargin));
   
3. if nargin==3
   
4.     if ischar(varargin{1})
   
5.         [lp,hp]=wavefilter(varargin{1},'d');
   
6.     else
   
7.         error('Missing wavelet name.');
   
8.     end
   
9. else
   
10.     lp=varargin{1}; hp=varargin{2};
    
11. end
    
12. fl=length(lp);sx=size(x);
    
13. if (ndims(x)~=2) | (min(sx)<2) | ~isreal(x) | ~isnumeric(x)
    
14.     error('X must be a real ,numeric matric.');
    
15. end
    
16. if(ndims(lp)~=2) | ~ isreal(lp) | ~isnumeric(lp) | (ndims(hp) ~=2) | ~ isreal(hp) | ~isnumeric(hp) ...
    
17.         | (fl~=length(hp)) | rem(fl,2)~=0
    
18.     error(['LP and HP must be ever and equal length real numeric filter vector.']);
    
19. end
    
20. if ~isreal(n) | ~isnumeric(n) |(n<1) |(n>log2(max(sx)))
    
21.     error(['N must be a real scalar between 1 and log2(max(size(X))).']);
    
22. end
    
23. c=[];s=sx;app=double(x);
    
24. for i=1:n
    
25.     [app,keep]=symextend(app,fl);
    
26.     rows=symconv(app,hp,'row',fl,keep);
    
27.     coefs=symconv(rows,hp,'col',fl,keep);
    
28.     c=[coefs(:)' c]; s=[size(coefs);s];
    
29.     coefs=symconv(rows,lp,'col',fl,keep);
    
30.     c=[coefs(:)' c];
    
31.     rows=symconv(app,lp,'row',fl,keep);
    
32.     coefs=symconv(rows,hp,'col',fl,keep);
    
33.     c=[coefs(:)' c];
    
34.     app=symconv(rows,lp,'col',fl,keep);
    
35. end
    
36. c=[app(:)' c]; s=[size(app);s];
    
37. function [y, keep]=symextend(x,fl)
    
38. keep=floor((fl+size(x)-1)/2);
    
39. y=padarray(x,[(fl-1) (fl-1)],'symmetric','both');
    
40. function y=symconv(x,h,type,fl,keep)
    
41. if strcmp(type, 'row')
    
42.     y=conv2(x,h);
    
43.     y=y(:,1:2:end);
    
44.     y=y(:,fl/2+1:fl/2+keep(2));
    
45. else
    
46.     y=conv2(x,h');
    
47.     y=y(1:2:end,:);
    
48.     y=y(fl/2+1:fl/2+keep(2),:);
    
49. end

复制代码

1. function [varargout]=waveback(c,s,varargin)
   
2. error(nargchk(3,5,nargin));
   
3. error(nargchk(1,2,nargout));
   
4. if(ndims(c) ~=2) | (size(c,1) ~=1)
   
5.     error('C must be a row vector .');
   
6. end
   
7. if (ndims(s) ~=2) | ~isreal(s) | ~isnumeric(s) | (size(s,2) ~=2)
   
8.     error('S must be a real,numeric two-column array.');
   
9. end
   
10. elements=prod(s,2);
    
11. if (length(c) <elements(end)) | ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))
    
12.     error(['[C S] must be a standard wavelet decomposition structure.']);
    
13. end
    
14. nmax=size(s,1)-2;
    
15. wname=varargin{1};filterchk=0;nchk=0;
    
16. switch nargin
    
17.     case 3
    
18.         if ischar(wname)
    
19.             [lp,hp]=wavefilter(wname,'r');
    
20.             n=nmax;
    
21.         else
    
22.                 error('Undefined filter.');
    
23.         end
    
24.         if nargout~=1
    
25.             error('Wrong number of output arguments.');
    
26.         end
    
27.     case 4
    
28.         if ischar(wname)
    
29.             [lp,hp]=wavefilter(wname,'r');
    
30.             n=varargin{2};nchk=1;
    
31.         else
    
32.             lp=varargin{1};
    
33.             hp=varargin{2};
    
34.             filterchk=1;n=nmax;
    
35.             if nargout ~=1
    
36.                 error('Wrong number of output arguments.');
    
37.             end
    
38.         end
    
39.     case 5
    
40.         lp=varargin{1};hp=varargin{2};filterchk=1;
    
41.         n=varargin{3};nchk=1;
    
42.     otherwise
    
43.         error('Improper number of input arguments.');
    
44. end
    
45. fl=length(lp);
    
46. if filterchk
    
47.     if (ndims(lp) ~=2) | ~isreal(lp) | ~isnueric(lp) |(ndims(hp) ~=2) | ~isreal(hp) ...
    
48.             | ~isnumeric(hp) |(fl ~=length(hp)) | rem(fl,2) ~=0
    
49.         error(['LP and HP must be even and equal length real,numeric filter vectors.']);
    
50.     end
    
51. end
    
52. if nchk & (~isnumeric(n) | ~isreal(n))
    
53.     error('N must be a real numeric.');
    
54. end
    
55. if(n~=nmax) & (nargout ~=2)
    
56.     error('Not enough output arguments.');
    
57. end
    
58. nc=c;ns=s;nnmax=nmax;
    
59. for i=1:n
    
60.     a=symconvup(wavecopy('a',nc,ns),lp,lp,fl,ns(3,:))+ ...
    
61.         symconvup(wavecopy('h',nc,ns,nnmax),hp,lp,fl,ns(3,:))+ ...
    
62.         symconvup(wavecopy('v',nc,ns,nnmax),lp,hp,fl,ns(3,:))+ ...
    
63.         symconvup(wavecopy('d',nc,ns,nnmax),hp,hp,fl,ns(3,:));
    
64.     nc=nc(4*prod(ns(1,:))+1:end);
    
65.     nc=[a(:)' nc];
    
66.     ns=ns(3:end,:);
    
67.     ns=[ns(1,:);ns];
    
68.     nnmax=size(ns,1)-2;
    
69. end
    
70. if nargout ==1
    
71.     a=nc; nc=repmat(0,ns(1,:)); nc(:)=a;
    
72. end
    
73. varargout{1}=nc;
    
74. if nargout==2
    
75.     varargout{2}=ns;
    
76. end
    
77. 
    
78. function z=symconvup(x,f1,f2,fln,keep)
    
79. y=zeros([2 1].*size(x));
    
80. y(1:2:end,:)=x;
    
81. y=conv2(y,f1');
    
82. z=zeros([1 2].*size(y));
    
83. z(:,1:2:end)=y;
    
84. z=conv2(z,f2);
    
85. z=z(fln-1:fln+keep(1)-2,fln-1:fln+keep(2)-2);

复制代码

simon21

1. %%第二代小波变换源码
   
2. %%  此程序用提升法实现第二代小波变换
   
3. %%  我用的是非整数阶小波变换
   
4. %%  采用时域实现,步骤先列后行
   
5. %%  正变换：分裂，预测，更新；
   
6. %%  反变换：更新，预测，合并
   
7. %%  只做一层（可以多层，而且每层的预测和更新方程不同）
   
8. 
   
9. clear;clc;
   
10. 
    
11. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
13. %%%%%  1.调原始图像矩阵
    
14. 
    
15. load wbarb;  %  下载图像
    
16. f=X;         %  原始图像
    
17. % f=[0 0 0 0 0 0 0 0 ;...
    
18. %    0 0 0 1 1 0 0 0 ;...
    
19. %    0 0 2 4 4 2 0 0 ;...
    
20. %    0 1 4 8 8 4 1 0 ;...
    
21. %    0 1 4 8 8 4 1 0 ;...
    
22. %    0 0 2 4 4 2 0 0 ;...
    
23. %    0 0 0 1 1 0 0 0 ;...
    
24. %    0 0 0 0 0 0 0 0 ;];  %  原始图像矩阵
    
25. 
    
26. N=length(f);         %  图像维数
    
27. T=N/2;               %  子图像维数
    
28. 
    
29. 
    
30. 
    
31. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
32. %%%%%%%%%%%%%%%%%正变换%%%%%%%%%%%%%%%%%%%%%%%%%
    
33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
34. 
    
35. 
    
36. 
    
37. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
38. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
39. %%%%   1.列变换
    
40. 
    
41. %  A.分裂（奇偶分开）
    
42. 
    
43. f1=f([1:2:N-1],:);  %  奇数
    
44. f2=f([2:2:N],:);    %  偶数
    
45. 
    
46. % f1(:,T+1)=f1(:,1);  %  补列
    
47. % f2(T+1,:)=f2(1,:);  %  补行
    
48. 
    
49. %  B.预测
    
50. 
    
51. for i_hc=1:T;
    
52.     high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
    
53. end;
    
54. 
    
55. % high_frequency_column(T+1,:)=high_frequency_column(1,:);  %  补行
    
56. 
    
57. %  C.更新
    
58. 
    
59. for i_lc=1:T;
    
60.     low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
    
61. end;
    
62. 
    
63. %  D.合并
    
64. f_column([1:1:T],:)=low_frequency_column([1:T],:);
    
65. f_column([T+1:1:N],:)=high_frequency_column([1:T],:);
    
66.    
    
67.    
    
68. figure(1)
    
69. colormap(map);
    
70. image(f);
    
71. 
    
72. figure(2)
    
73. colormap(map);
    
74. image(f_column);
    
75. 
    
76. 
    
77. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
78. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
79. %%%%   2.行变换
    
80. 
    
81. %  A.分裂（奇偶分开）
    
82. 
    
83. f1=f_column(:,[1:2:N-1]);  %  奇数
    
84. f2=f_column(:,[2:2:N]);    %  偶数
    
85. 
    
86. 
    
87. % f2(:,T+1)=f2(:,1);    %  补行
    
88. 
    
89. %  B.预测
    
90. 
    
91. for i_hr=1:T;
    
92.     high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
    
93. end;
    
94. 
    
95. % high_frequency_row(:,T+1)=high_frequency_row(:,1);  %  补行
    
96. 
    
97. %  C.更新
    
98. 
    
99. for i_lr=1:T;
    
100.     low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
     
101. end;
     
102. 
     
103. %  D.合并
     
104. f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);
     
105. f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);
     
106.    
     
107. 
     
108. figure(3)
     
109. colormap(map);
     
110. image(f_row);
     
111. 
     
112. 
     
113. 
     
114. 
     
115. 
     
116. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
117. %%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%%
     
118. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
119. 
     
120. 
     
121. 
     
122. 
     
123. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
124. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
125. %%%%   1.行变换
     
126. 
     
127. %   A.提取（低频高频分开）
     
128. 
     
129. f1=f_row(:,[T+1:1:N]);  %  奇数
     
130. f2=f_row(:,[1:1:T]);    %  偶数
     
131. 
     
132. 
     
133. % f2(:,T+1)=f2(:,1);    %  补行
     
134. 
     
135. %  B.更新
     
136. 
     
137. for i_lr=1:T;
     
138.     low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
     
139. end;
     
140. 
     
141. %  C.预测
     
142. 
     
143. for i_hr=1:T;
     
144.     high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
     
145. end;
     
146. 
     
147. % high_frequency_row(:,T+1)=high_frequency_row(:,1);  %  补行
     
148. 
     
149. 
     
150. %  D.合并(奇偶分开合并)
     
151. f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);
     
152. f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);
     
153.    
     
154. 
     
155. figure(4)
     
156. colormap(map);
     
157. image(f_row);
     
158. 
     
159. 
     
160. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
161. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
162. %%%%   2.列变换
     
163. 
     
164. %  A.提取（低频高频分开）
     
165. 
     
166. f1=f_row([T+1:1:N],:);  %  奇数
     
167. f2=f_row([1:1:T],:);    %  偶数
     
168. 
     
169. % f1(:,T+1)=f1(:,1);  %  补列
     
170. % f2(T+1,:)=f2(1,:);  %  补行
     
171. 
     
172. %  B.更新
     
173. 
     
174. for i_lc=1:T;
     
175.     low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
     
176. end;
     
177. 
     
178. %  C.预测
     
179. 
     
180. for i_hc=1:T;
     
181.     high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
     
182. end;
     
183. 
     
184. % high_frequency_column(T+1,:)=high_frequency_column(1,:);  %  补行
     
185. 
     
186. 
     
187. 
     
188. %  D.合并(奇偶分开合并）
     
189. f_column([2:2:N],:)=low_frequency_column([1:T],:);
     
190. f_column([1:2:N-1],:)=high_frequency_column([1:T],:);
     
191.    
     
192.    
     
193. figure(5)
     
194. colormap(map);
     
195. image(f_column);

复制代码

simon21

1. N=10000;
   
2. s=zeros(1,N);
   
3. for n=1:N        
   
4.     if n<0.4*N||n>0.8*N      
   
5.         s(n)=31.1*sin(2*pi*50/10000*n);   
   
6.     else        
   
7.         s(n)=22.5*sin(2*pi*50/10000*n);     
   
8.     end
   
9. end
   
10. l=length(s);
    
11. [c,l]=wavedec(s,6,'db5'); %用db5小波分解信号到第六层
    
12. subplot(8,1,1);
    
13. plot(s);
    
14. title('用db5小波分解六层：s=a6+d6+d5+d4+d3+d2+d1');
    
15. Ylabel('s');%对分解结构【c，l】中第六层低频部分进行重构
    
16. a6=wrcoef('a',c,l,'db5',6);
    
17. subplot(8,1,2);
    
18. plot(a6);
    
19. Ylabel('a6');%对分解结构【c，l】中各层高频部分进行重构
    
20. for i=1:6   
    
21.     decmp=wrcoef('d',c,l,'db5',7-i);   
    
22.     subplot(8,1,i+2);   
    
23.     plot(decmp);   
    
24.     Ylabel(['d',num2str(7-i)]);
    
25. end
    
26. %-----------------------------------------------------------
    
27. rec=zeros(1,300);
    
28. rect=zeros(1,300);
    
29. ke=1;
    
30. u=0;
    
31. d1=wrcoef('d',c,l,'db5',1);
    
32. figure(2);
    
33. plot(d1);
    
34. si=0;
    
35. N1=0;
    
36. N0=0;
    
37. sce=0;
    
38. for n=20:N-30
    
39.   rect(ke)=s(n);
    
40.   ke=ke+1;
    
41.   if(ke>=301)
    
42.         if(si==2)
    
43.             rec=rect;
    
44.             u=2;
    
45.         end;
    
46.         si=0;
    
47.         ke=1;
    
48.     end;
    
49. if(d1(n)>0.01)                       % the condition of abnormal append.
    
50.         N1=n;
    
51.         if(N0==0)
    
52.             N0=n;
    
53.             si=si+1;
    
54.         end;
    
55.         if(N1>N0+30)
    
56.             Nlen=N1-N0;
    
57.             Tab=Nlen/10000;
    
58.         end;
    
59.     end;
    
60. 
    
61. if(si==1)
    
62.     for k=N0:N0+99                       %testing of 1/4 period signals to
    
63.        sce=sce+s(k)*s(k)/10000;
    
64.     end;
    
65.     re=sqrt(sce*200)                   %re indicate the pike value of .
    
66.     sce=0;
    
67.     si=si+1;
    
68.   end;
    
69. end;
    
70.   Nlen
    
71.   N0
    
72.   n=1:300;
    
73.   figure(3)
    
74.   plot(n,rec);

复制代码

simon21

1. function [T_img,Sub_T]=threshold_2_N(img,levels)
   
2. 
   
3. % reference :image denoising using wavelet thresholding
   
4. 
   
5. [xx,yy]=size(img);
   
6. HH=img((xx/2+1):xx,(yy/2+1):yy);
   
7. 
   
8. delt_2=(std(HH(:)))^2;%(median(abs(HH(:)))/0.6745)^2;%
   
9. 
   
10. T_img=img;
    
11. 
    
12. for i=1:levels
    
13.       
    
14.     temp_x=xx/2^i;
    
15.     temp_y=yy/2^i;
    
16. %     belt=1.0*(log(temp_x/(2*levels)))^0.5;
    
17.      belt=1.0*(log(temp_x/(2*levels)))^0.5;  %2.5 0.8
    
18.    %HL
    
19.     HL=img(1:temp_x,(temp_y+1):2*temp_y);
    
20.    
    
21.     delt_y=std(HL(:));
    
22.     T_1=belt*delt_2/delt_y;
    
23.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
24.     T_HL=sign(HL).*max(0,abs(HL)-T_1);
    
25. 
    
26.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
27.     T_img(1:temp_x,(temp_y+1):2*temp_y)=T_HL;
    
28.    
    
29.     Sub_T(3*(i-1)+1)=T_1;
    
30.       
    
31.    %LH
    
32.     LH=img((temp_x+1):2*temp_x,1:temp_y);
    
33.    
    
34.     delt_y=std(LH(:));
    
35.     T_2=belt*delt_2/delt_y;
    
36.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
37.     T_LH=sign(LH).*max(0,abs(LH)-T_2);
    
38.    
    
39.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
40.     T_img((temp_x+1):2*temp_x,1:temp_y)=T_LH;
    
41.    
    
42.     Sub_T(3*(i-1)+2)=T_2;
    
43.    
    
44.    %HH
    
45.     HH=img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y);
    
46.    
    
47.     delt_y=std(HH(:));
    
48.     T_3=belt*delt_2/delt_y;
    
49.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
50.     T_HH=sign(HH).*max(0,abs(HH)-T_3);
    
51.    
    
52.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
53.     T_img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y)=T_HH;
    
54.    
    
55.     Sub_T(3*(i-1)+3)=T_3;
    
56. end

复制代码

simon21

1. clear all;
   
2. load wbarb;
   
3. I = ind2gray(X,map);imshow(I);
   
4. I1 = imadjust(I,stretchlim(I),[0,1]);figure;imshow(I1);
   
5. [N,M] = size(I);
   
6. 
   
7. h = [0.125,0.375,0.375,0.125];
   
8. g = [0.5,-0.5];
   
9. delta = [1,0,0];
   
10. 
    
11. J = 3;
    
12. 
    
13. a(1:N,1:M,1,1:J+1) = 0;
    
14. dx(1:N,1:M,1,1:J+1) = 0;
    
15. dy(1:N,1:M,1,1:J+1) = 0;
    
16. d(1:N,1:M,1,1:J+1) = 0;
    
17. 
    
18. a(:,:,1,1) = conv2(h,h,I,'same');
    
19. dx(:,:,1,1) = conv2(delta,g,I,'same');
    
20. dy(:,:,1,1) = conv2(g,delta,I,'same');
    
21. 
    
22. x = dx(:,:,1,1);
    
23. y = dy(:,:,1,1);
    
24. d(:,:,1,1) = sqrt(x.^2+y.^2);
    
25. I1 = imadjust(d(:,:,1,1),stretchlim(d(:,:,1,1)),[0 1]);figure;imshow(I1);
    
26. 
    
27. lh = length(h);
    
28. lg = length(g);
    
29. 
    
30. for j = 1:J+1
    
31.   lhj = 2^j*(lh-1)+1;
    
32.   lgj = 2^j*(lg-1)+1;
    
33.   hj(1:lhj)=0;
    
34.   gj(1:lgj)=0;
    
35.   for n = 1:lh
    
36.     hj(2^j*(n-1)+1)=h(n);
    
37.   end
    
38. 
    
39.   for n = 1:lg
    
40.     gj(2^j*(n-1)+1)=g(n);
    
41.   end
    
42.   
    
43.   a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
    
44.   dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same');
    
45.   dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same');
    
46. 
    
47.   x = dx(:,:,1,j+1);
    
48.   y = dy(:,:,1,j+1);
    
49.   dj(:,:,1,j+1) = sqrt(x.^2+y.^2);
    
50. 
    
51.   I1 = imadjust(dj(:,:,1,j+1),stretchlim(dj(:,:,1,j+1)),[0 1]);figure;imshow(I1);
    
52. end  

复制代码

simon21

1. close all
   
2. clear all
   
3. clc;
   
4. 
   
5. Original_image = imread('lena512_marked.bmp');
   
6. Original_image = double(Original_image);
   
7. [ll1,h1,v1]=ADWT2(Original_image,1);
   
8. HL1 = v1{1};
   
9. HL1_1D = sort(HL1(:)).';
   
10. bin_size = 2;
    
11. range_value = 81;
    
12. step =[-range_value:bin_size:range_value];
    
13. Extract_histogram1 = histc(HL1_1D,step);  % 提取的 histogram 步长为1
    
14. % Extract_histogram1 = histc(HL1_1D,[floor(min(HL1_1D)):2:(floor(max(HL1_1D))+1)]);
    
15. [m1,n1] = size(Extract_histogram1);
    
16. relation1 = zeros(1,n1);
    
17. for i = 2:3:(n1/2)
    
18.     relation1(1,i) = 2*Extract_histogram1(1,i)/(Extract_histogram1(1,i+1)+Extract_histogram1(1,i-1));
    
19. end
    
20. for i = ((n1/2)+4):3:(n1-2)
    
21.     relation1(1,i) = 2*Extract_histogram1(1,i)/(Extract_histogram1(1,i+1)+Extract_histogram1(1,i-1));
    
22. end
    
23. 
    
24. %%% Extract watermark W
    
25. W_number1 = 24;         % 水印的个数 W_number
    
26. H1 = zeros(1,3*W_number1);
    
27. H1(1:(3*W_number1/2)) = Extract_histogram1(1:(3*W_number1/2));  % H的第一个直方图的起点是-range_value
    
28. H1((3*W_number1/2+1):3*W_number1) = Extract_histogram1((n1/2+3):(n1/2+2+3*W_number1/2));% H的第（2*W_number+1）个直方图的起点是-range_value+bin_size*n/2
    
29. a = zeros(1,W_number1);
    
30. b = zeros(1,W_number1);
    
31. c = zeros(1,W_number1);
    
32. H_matrix1 = reshape(H1,3,W_number1);
    
33. a = H_matrix1(1,:);
    
34. b = H_matrix1(2,:);
    
35. c = H_matrix1(3,:);    % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
    
36. extract_W1=zeros(1,W_number1);
    
37. for i=1:W_number1
    
38.     if (2*b(1,i) >= a(1,i)+c(1,i))
    
39.         extract_W1(1,i)=1;
    
40.     else if (2*b(1,i) < (a(1,i)+c(1,i)))
    
41.         extract_W1(1,i)=-1;
    
42.         end
    
43.     end
    
44. end
    
45. extract_watermark1=extract_W1
    
46. 
    
47. %%% 缩小
    
48. Original_image = imread('lena512_marked.bmp');
    
49. % scaled_image = imresize(Original_image,[450,450]);
    
50. scaled_image = imresize(Original_image,.9);
    
51. scaled_image = double(scaled_image);
    
52. [ll2,h2,v2]=ADWT2(scaled_image,1);
    
53. HL2 = v2{1};
    
54. HL2_1D = sort(HL2(:)).';    % 将HL子带变为一维向量
    
55. Extract_histogram2 = histc(HL2_1D,step);    % 提取的 histogram 步长为1
    
56. [m2,n2] = size(Extract_histogram2);
    
57. % hist(sort(HL(:)),n);
    
58. relation2 = zeros(1,n1);
    
59. for i = 2:3:(n1/2)
    
60.     relation2(1,i) = 2*Extract_histogram2(1,i)/(Extract_histogram2(1,i+1)+Extract_histogram2(1,i-1));
    
61. end
    
62. for i = ((n1/2)+4):3:(n1-2)
    
63.     relation2(1,i) = 2*Extract_histogram2(1,i)/(Extract_histogram2(1,i+1)+Extract_histogram2(1,i-1));
    
64. end
    
65. 
    
66. %%% Extract watermark W
    
67. W_number2 = 24;         % 水印的个数 W_number
    
68. H2 = zeros(1,3*W_number2);
    
69. H2(1:(3*W_number2/2)) = Extract_histogram2(1:(3*W_number2/2));  % H的第一个直方图的起点是-range_value
    
70. H2((3*W_number2/2+1):3*W_number2) = Extract_histogram2((n2/2+3):(n2/2+2+3*W_number2/2));% H的第（2*W_number+1）个直方图的起点是-range_value+bin_size*n/2
    
71. a = zeros(1,W_number2);
    
72. b = zeros(1,W_number2);
    
73. c = zeros(1,W_number2);
    
74. H_matrix2 = reshape(H2,3,W_number2);
    
75. a = H_matrix2(1,:);
    
76. b = H_matrix2(2,:);
    
77. c = H_matrix2(3,:);    % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
    
78. extract_W2=zeros(1,W_number2);
    
79. for i=1:W_number2
    
80.     if (2*b(1,i) >= a(1,i)+c(1,i))
    
81.         extract_W2(1,i)=1;
    
82.     else if (2*b(1,i) < (a(1,i)+c(1,i)))
    
83.         extract_W2(1,i)=-1;
    
84.         end
    
85.     end
    
86. end
    
87. extract_watermark2=extract_W2
    
88. 
    
89. 
    
90. %%% 放大
    
91. Original_image = imread('lena512_marked.bmp');
    
92. scaled_image = imresize(Original_image,1.3);
    
93. scaled_image = double(scaled_image);
    
94. [ll3,h3,v3]=ADWT2(scaled_image,1);
    
95. HL3 = v3{1};
    
96. HL3_1D = sort(HL3(:)).';    % 将HL子带变为一维向量
    
97. Extract_histogram3 = histc(HL3_1D,step);    % 提取的 histogram 步长为1
    
98. [m3,n3] = size(Extract_histogram3);
    
99. % hist(sort(HL(:)),n);
    
100. relation3 = zeros(1,n1);
     
101. for i = 2:3:(n1/2)
     
102.     relation3(1,i) = 2*Extract_histogram3(1,i)/(Extract_histogram3(1,i+1)+Extract_histogram3(1,i-1));
     
103. end
     
104. for i = ((n1/2)+4):3:(n1-2)
     
105.     relation3(1,i) = 2*Extract_histogram3(1,i)/(Extract_histogram3(1,i+1)+Extract_histogram3(1,i-1));
     
106. end
     
107. 
     
108. %%% Extract watermark W
     
109. W_number3 = 24;         % 水印的个数 W_number
     
110. H3 = zeros(1,3*W_number3);
     
111. H3(1:(3*W_number3/2)) = Extract_histogram3(1:(3*W_number3/2));  % H的第一个直方图的起点是-range_value
     
112. H3((3*W_number3/2+1):3*W_number3) = Extract_histogram3((n3/2+3):(n3/2+2+3*W_number3/2));% H的第（2*W_number+1）个直方图的起点是-range_value+bin_size*n/2
     
113. a = zeros(1,W_number3);
     
114. b = zeros(1,W_number3);
     
115. c = zeros(1,W_number3);
     
116. H_matrix3 = reshape(H3,3,W_number3);
     
117. a = H_matrix3(1,:);
     
118. b = H_matrix3(2,:);
     
119. c = H_matrix3(3,:);    % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
     
120. extract_W3=zeros(1,W_number3);
     
121. for i=1:W_number3
     
122.     if (2*b(1,i) >= a(1,i)+c(1,i))
     
123.         extract_W3(1,i)=1;
     
124.     else if (2*b(1,i) < (a(1,i)+c(1,i)))
     
125.         extract_W3(1,i)=-1;
     
126.         end
     
127.     end
     
128. end
     
129. extract_watermark3=extract_W3
     
130. 
     
131. %%% 旋转
     
132. Original_image = imread('lena512_marked.bmp');
     
133. Rescale_image = imrotate(Original_image,10);
     
134. Rescale_image = double(Rescale_image);
     
135. [ll4,h4,v4]=ADWT2(Rescale_image,1);
     
136. HL4 = v4{1};
     
137. HL4_1D = sort(HL4(:)).';    % 将HL子带变为一维向量
     
138. Extract_histogram4 = histc(HL4_1D,step);    % 提取的 histogram 步长为1
     
139. [m4,n4] = size(Extract_histogram4);
     
140. % hist(sort(HL(:)),n);
     
141. relation4 = zeros(1,n1);
     
142. for i = 2:3:(n1/2)
     
143.     relation4(1,i) = 2*Extract_histogram4(1,i)/(Extract_histogram4(1,i+1)+Extract_histogram4(1,i-1));
     
144. end
     
145. for i = ((n1/2)+4):3:(n1-2)
     
146.     relation4(1,i) = 2*Extract_histogram4(1,i)/(Extract_histogram4(1,i+1)+Extract_histogram4(1,i-1));
     
147. end
     
148. 
     
149. %%% Extract watermark W
     
150. W_number4 = 24;         % 水印的个数 W_number4
     
151. H4 = zeros(1,3*W_number4);
     
152. H4(1:(3*W_number4/2)) = Extract_histogram4(1:(3*W_number4/2));  % H的第一个直方图的起点是-range_value
     
153. H4((3*W_number4/2+1):3*W_number4) = Extract_histogram4((n4/2+3):(n4/2+2+3*W_number4/2));% H的第（2*W_number+1）个直方图的起点是-range_value+bin_size*n/2
     
154. a = zeros(1,W_number4);
     
155. b = zeros(1,W_number4);
     
156. c = zeros(1,W_number4);
     
157. H_matrix4 = reshape(H4,3,W_number4);
     
158. a = H_matrix4(1,:);
     
159. b = H_matrix4(2,:);
     
160. c = H_matrix4(3,:);    % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
     
161. extract_W4=zeros(1,W_number4);
     
162. for i=1:W_number4
     
163.     if (2*b(1,i) >= a(1,i)+c(1,i))
     
164.         extract_W4(1,i)=1;
     
165.     else if (2*b(1,i) < (a(1,i)+c(1,i)))
     
166.         extract_W4(1,i)=-1;
     
167.         end
     
168.     end
     
169. end
     
170. extract_watermark4=extract_W4
     
171. 
     
172. X1=2:1:(n1-2);
     
173. Y1=relation1(1,X1);
     
174. X2=2:1:(n1-2);
     
175. Y2=relation2(1,X2);
     
176. X3=2:1:(n1-2);
     
177. Y3=relation3(1,X3);
     
178. X4=2:1:(n1-2);
     
179. Y4=relation4(1,X4);
     
180. figure;plot(X1,Y1,'r*',X2,Y2,'b+',X3,Y3,'kd',X4,Y4,'gd');
     
181. title('The estimation of bin relation');
     

复制代码

simon21

该程序是用小波函数构建神经网络的源程序。用以分析心电信号、脑电信号等等。

1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
2. clear all
   
3. load data
   
4. M1=20;
   
5. epo=15;
   
6. A=4;
   
7. B=18;
   
8. B2=B/2+1
   
9. N=500;
   
10. M=(A+1)*(B+1);
    
11. for  a0=1:A+1;
    
12. for  b0=1:B+1;
    
13. i=(B+1)*(a0-1)+b0;
    
14. b_init(i)=((b0-B2)/10)/(2^(-A)); a_init(i)=1/(2^(-A));
    
15. c_init(i)=(20-A)/2;
    
16. end
    
17. end
    
18. w0=ones(1,M);
    
19.           for i=1:N            
    
20.               for j=1:M           
    
21.                   t=x(i);
    
22.                   t= a_init(j)*t-b_init(j);
    
23.                  %P0(i,j)= (cos(1.75*t)*exp(-t*t/2))/2^c_init(j);
    
24.                  P0(i,j)= ((1-t*t)*exp(-t*t/2))/2^c_init(j);
    
25.              end
    
26.           end
    
27. %calculation of output of network
    
28.            for i=1:N
    
29.                  u=0;
    
30.                  for j=1:M
    
31.                      u=u+w0(j)*P0(i,j);%w0?aè¨?μ
    
32.                  end
    
33.                  y0(i)= u;% y(p)= u=??W(j)*phi(p,j)= ??W(j)* |μj(t)
    
34.            end
    
35. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
36. for k=1:M
    
37.     W(:,k)=P0(:,k);
    
38. end
    
39. for k=2:M
    
40.     u=0;
    
41.    
    
42.      for      i=1:k-1   
    
43.               aik(i)=(P0(:,k)'*W(:,i))/(W(:,i)'*W(:,i));
    
44.               u=u+aik(i) *W(:,i);
    
45.               
    
46.      end
    
47.      W(:,k)=P0(:,k)-u;
    
48. end
    
49. 
    
50. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
51. 
    
52. for   i=1:M                        
    
53.      g(i)= (W(:,i)'*d')/( W(:,i)'* W(:,i));
    
54. end
    
55. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
56. u=0;
    
57. for    i=1:M
    
58.        u=u+g(i)*(W(:,i)'*W(:,i));
    
59. end
    
60. 
    
61. DD=u;
    
62. 
    
63. for    i=1:M
    
64.        Erro(i)=(g(i)^2)*(W(:,i)'*W(:,i))/DD;
    
65. end
    
66. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
67. k=1;
    
68. 
    
69. while(k<=M1)
    
70.     u=1;
    
71.     for  i=2:M
    
72.          if   abs(Erro(u))<abs(Erro(i));
    
73.               u=i
    
74.          else u=u
    
75.          end
    
76.     end
    
77.     I(k)=u;
    
78.     Erro(u)=0
    
79. k=k+1;  
    
80. end
    
81. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
82. for k=1:M1
    
83.     u=I(k);
    
84.     a(k)=a_init(u);
    
85.     b(k)=b_init(u);
    
86.     c(k)=c_init(u);
    
87.     w1(k)=w0(u);
    
88. end
    
89. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
90. epoch=1;
    
91. 
    
92. error=0.1;
    
93. err=0.0001;
    
94. lin=0.5;
    
95. 
    
96. while (error>=err & epoch<=epo)
    
97.           for i=1:N            
    
98.               for j=1:M1        
    
99.                   t=x(i);
    
100.                   t= a(j)*t-b(j);
     
101.                  %P1(i,j)= (cos(1.75*t)*exp(-t*t/2))/2^c(j);
     
102.                  P1(i,j)= ((1-t*t)*exp(-t*t/2))/2^c(j);
     
103.              end
     
104.           end
     
105. %calculation of output of network
     
106.            for i=1:N
     
107.                  u=0;
     
108.                  for j=1:M1
     
109.                      u=u+w1(j)*P1(i,j);
     
110.                  end
     
111.                  y1(i)= u;
     
112.            end
     
113. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
114.    u=0;
     
115.      for i=1:N
     
116.             u=u+(d(i)-y1(i))^2;
     
117.      end
     
118.     u=u/2;%u=1/2??(d-p)^2
     
119. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
120. if   u>error
     
121. lin=lin*0.8;
     
122. else
     
123. lin=lin*1.2;
     
124. end
     
125. error=u; %error=u=1/2??(d-p)^2
     
126. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
127.          for j=1:M1   
     
128.              u=0;
     
129.              for i=1:N
     
130.                  u=u+(d(i)-y1(i))*P1(i,j);
     
131.              end
     
132.              EW(j)=-u;
     
133.           end
     
134.    if    epoch==1
     
135.        SW=-EW;
     
136.        w1_=w1;
     
137.       
     
138.    else
     
139.    SW=-EW+((EW*EW')*SW_)/(EW_*EW_');
     
140.    end
     
141.    EW_=EW;
     
142.    SW_=SW;
     
143.     w1=w1_+SW*lin;
     
144.     w1_=w1;
     
145.   
     
146.     %number of epoch increase by 1
     
147. epoch=epoch+1;
     
148. end
     
149. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
150. subplot(2,1,1);plot(x,d);
     
151. subplot(2,1,2);plot(x,y1);
     

复制代码

simon21

1. x=wavread('vqdeout.wav');
   
2. subplot(3,2,1);
   
3. plot(x);
   
4. [c,l]=wavedec(x,3,'db1');
   
5. cA3=appcoef(c,l,'db1',3);
   
6. [cD1,cD2,cD3] = detcoef(c,l,[1 2 3]);
   
7. subplot(3,2,2);
   
8. plot(cA3);
   
9. subplot(3,2,3);
   
10. plot(cD1);
    
11. subplot(3,2,4);
    
12. plot(cD2);
    
13. subplot(3,2,5);
    
14. plot(cD3);
    
15. 
    
16. %quantize
    
17. [partition,codebook]=py_quantize(cA3);                  
    
18. [index,quants]=quantiz(cA3,partition,codebook);                 
    
19. [partition1,codebook1]=py_quantize(cD1);
    
20. [index1,quants1]=quantiz(cD1,partition1,codebook1);               
    
21. [partition2,codebook2]=py_quantize(cD2);
    
22. [index2,quants2]=quantiz(cD2,partition2,codebook2);               
    
23. [partition3,codebook3]=py_quantize(cD3);
    
24. [index3,quants3]=quantiz(cD3,partition3,codebook3);               
    
25. 
    
26. %probability
    
27. g=py_probability(quants,l(1,1),codebook);
    
28. g1=py_probability(quants1,l(4,1),codebook1);
    
29. g2=py_probability(quants2,l(3,1),codebook2);
    
30. g3=py_probability(quants3,l(2,1),codebook3);
    
31. 
    
32. 
    
33. 
    
34. %huffman编解码
    
35. [dict,avglen] = huffmandict(codebook,g);
    
36. comp = huffmanenco(quants,dict); % Encode the data.
    
37. dsig = huffmandeco(comp,dict);
    
38. isequal(quants,dsig)
    
39. 
    
40. %huffman编解码
    
41. [dict1,avglen1] = huffmandict(codebook1,g1);
    
42. comp1 = huffmanenco(quants1,dict1); % Encode the data.
    
43. dsig1 = huffmandeco(comp1,dict1);
    
44. isequal(quants1,dsig1)
    
45. 
    
46. %huffman编解码
    
47. [dict2,avglen2] = huffmandict(codebook2,g2);
    
48. comp2 = huffmanenco(quants2,dict2); % Encode the data.
    
49. dsig2 = huffmandeco(comp2,dict2);
    
50. isequal(quants2,dsig2)
    
51. 
    
52. %huffman编解码
    
53. [dict3,avglen3] = huffmandict(codebook3,g3);
    
54. comp3 = huffmanenco(quants3,dict3); % Encode the data.
    
55. dsig3 = huffmandeco(comp3,dict3);
    
56. isequal(quants3,dsig3)
    
57. 
    
58. %dsig3=dsig3';
    
59. %dsig2=dsig2';
    
60. %dsig1=dsig1';
    
61. %dsig=dsig';
    
62. %重构
    
63. cA2=idwt(dsig,dsig3,'db1');
    
64. icA2(1:13950)=cA2(1:13950);
    
65. cA1=idwt(icA2,dsig2,'db1');
    
66. icA1(1:27899)=cA1(1:27899);
    
67. XX=idwt(icA1,dsig1,'db1');
    
68. XX=XX';
    
69. subplot(3,2,6);
    
70. plot(XX);
    
71. 
    
72. MSE=sum((x-XX).^2)/length(XX);
    
73. PSNR=10*log10((max(max(XX)))^2/MSE);
    

复制代码

simon21

1. 
   
2. %File name :       nprogram.m
   
3. %Description   :   This   file   reads   the  data   from   its   source   into   their   respective  matrices   prior   to
   
4.   %                  performing wavelet decomposition.
   
5. 
   
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
7. %Clear command screen and variables
   
8. clc;
   
9. clear;
   
10. 
    
11. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. % user desired resolution level (Tested: resolution = 2 is best)
    
13. level = menu('Enter desired resolution level: ', '1',...
    
14.                   '2 (Select this for testing)', '3', '4');
    
15. switch level
    
16.       case 1, resolution = 1;
    
17.       case 2, resolution = 2;
    
18.       case 3, resolution = 3;
    
19.       case 4, resolution = 4;
    
20. end
    
21. 
    
22. msg = ['Resolution level to be used is ', num2str(resolution)];
    
23. disp(msg);
    
24. 
    
25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
26. % user desired amount of data to use
    
27. data = menu('Choose amount of data to use: ', '1 day', '2 days', '3 days', '4 days',...
    
28.                  '5 days', '6 days', '1 week (Select this for testing)');
    
29. switch data
    
30.       case 1, dataPoints = 48;        %1 day = 48 points
    
31.       case 2, dataPoints = 96;        %2 days = 96 points
    
32.       case 3, dataPoints = 144;       %3 days = 144 points
    
33.       case 4, dataPoints = 192;       %4 days = 192 points
    
34.       case 5, dataPoints = 240;       %5 days = 240 points
    
35.       case 6, dataPoints = 288;       %6 days = 288 points
    
36.       case 7, dataPoints = 336;       %1 weeks = 336 points
    
37. 
    
38. end
    
39. 
    
40. msg = ['No. of data points to be used is ', num2str(dataPoints)];
    
41. disp(msg);
    
42. 
    
43. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
44. %Menu for data set selection
    
45. select = menu('Use QLD data of: ', 'Jan02',...
    
46.                     'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02');
    
47. switch select
    
48.       case 1, demandFile = 'DATA200601_QLD1';
    
49. 
    
50.       case 2, demandFile = 'DATA200602_QLD1';
    
51. 
    
52.       case 3, demandFile = 'DATA200603_QLD1';
    
53. 
    
54.       case 4, demandFile = 'DATA200604_QLD1';
    
55. 
    
56.       case 5, demandFile = 'DATA200605_QLD1';
    
57. end
    
58. 
    
59. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
60. %Reading the historical DEMAND data into tDemandArray
    
61. selectedDemandFile=[demandFile,'.csv'];
    
62. [regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...
    
63. = textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
    
64. 
    
65. %Display no. of points in the selected time series demand data
    
66. [demandDataPoints, y] = size(tDemandArray);
    
67. msg = ['The no. of points in the selected Demand data is ', num2str(demandDataPoints)];
    
68. disp(msg);
    
69. 
    
70. %Decompose historical demand data signal
    
71. [dD, l] = swtmat(tDemandArray, resolution, 'db2');
    
72. approx = dD(resolution, :);
    
73. 
    
74. %Plot the original demand data signal
    
75. figure (1);
    
76. subplot(resolution + 2, 1, 1); plot(tDemandArray(1: dataPoints))
    
77. legend('Demand Original');
    
78. title('QLD Demand Data Signal');
    
79. 
    
80. %Plot the approximation demand data signal
    
81. for i = 1 : resolution
    
82.       subplot(resolution + 2, 1, i     + 1); plot(approx(1: dataPoints))
    
83.       legend('Demand Approximation');
    
84. end
    
85. 
    
86. %After displaying approximation signal, display detail x
    
87. for i = 1: resolution
    
88.     if( i > 1 )
    
89.             detail(i, :) = dD(i-1, :)- dD(i, :);
    
90.       else
    
91.             detail(i, :) = tDemandArray' - dD(1, :);
    
92.       end
    
93. 
    
94.            if i == 1
    
95.                  subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
    
96.                  legendName = ['Demand Detail ', num2str(i)];
    
97.                   legend(legendName);
    
98. 
    
99.             else
    
100.                  subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
     
101.                  legendName = ['Demand Detail ', num2str(i)];
     
102.                  legend(legendName);
     
103. 
     
104.             end
     
105. 
     
106.     i = i + 1;
     
107. end
     
108. 
     
109. %Normalising approximation demand data
     
110. maxDemand = max(approx'); %Find largest component
     
111. normDemand = approx ./ maxDemand; %Right divison
     
112. 
     
113. maxDemandDetail = [ ];
     
114. normDemandDetail = [, ];
     
115. 
     
116. detail = detail + 4000;
     
117. 
     
118. for i = 1: resolution
     
119.       maxDemandDetail(i) = max(detail(i, :));
     
120.       normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i);
     
121. end
     
122. 
     
123. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
124. %Reading the historical historical PRICE data into rrpArray
     
125. selectedPriceFile = [demandFile, '.csv'];
     
126. [regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...
     
127. = textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
     
128. 
     
129. %Display no. of points in Price data
     
130. [noOfDataPoints, y] = size(rrpArray);
     
131. msg = ['The no. of points in Price data data is ', num2str(noOfDataPoints)];
     
132. disp(msg);
     
133. 
     
134. %Decompose historical Price data signal
     
135. [dP, l] = swtmat(rrpArray, resolution, 'db2');
     
136. approxP = dP(resolution, :);
     
137. 
     
138. %Plot the original Price data signal
     
139. figure (2);
     
140. subplot(resolution + 3, 1, 1); plot(rrpArray(2: dataPoints))
     
141. legend('Price Original');
     
142. title('Price Data Signal');
     
143. 
     
144. %Plot the approximation Price data signal
     
145. for i = 1 : resolution
     
146.       subplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints))
     
147.       legend('Price Approximation');
     
148. end
     
149. 
     
150. %After displaying approximation signal, display detail x
     
151. for i = 1: resolution
     
152.      if( i > 1 )
     
153.             detailP(i, :) = dP(i-1, :)- dP(i, :);
     
154.       else
     
155.             detailP(i, :) = rrpArray' - dP(1, :);
     
156.       end
     
157. 
     
158. if i == 1
     
159.       [B,A]=butter(1,0.65,'low');
     
160.       result =filter(B,A, detailP(i, 1: dataPoints));
     
161. 
     
162.       subplot(resolution + 3, 1, resolution - i + 4);plot(result(i, 2: dataPoints))
     
163.       legendName = ['low pass filter', num2str(i)];
     
164.       legend(legendName);
     
165. 
     
166.       subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
     
167.       legendName = ['Price Detail ', num2str(i)];
     
168.       legend(legendName);
     
169. 
     
170. else
     
171.       subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
     
172.       legendName = ['Price Detail ', num2str(i)];
     
173.       legend(legendName);
     
174. 
     
175. end
     
176.      i = i + 1;
     
177. end
     
178. 
     
179. %Normalising approximation Price data
     
180. maxPrice = max(approxP'); %Find largest component
     
181. normPrice = approxP ./ maxPrice; %Right divison
     
182. 
     
183. maxPriceDetail = [ ];
     
184. normPriceDetail = [, ];
     
185. 
     
186. detailP = detailP + 40;
     
187. 
     
188. for i = 1: resolution
     
189.       maxPriceDetail(i) = max(detailP(i, :));
     
190.       normPriceDetail(i, :) = detailP(i, :) ./maxPriceDetail(i);
     
191. end
     
192. 
     
193. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
194. %Function here allows repetitive options to,
     
195. %      1) Create new NNs, 2) Retrain the existing NNs,
     
196. %      3) Perform load demand forecasting and 4) Quit session
     
197. 
     
198. while (1)
     
199. 
     
200.       choice = menu('Please select one of the following options: ',...
     
201.                           'CREATE new Neural Networks',...
     
202.                           'Perform FORECASTING of load demand', 'QUIT session...');
     
203.       switch choice
     
204.            case 1, scheme = '1';
     
205.            case 2, scheme = '2';
     
206.            case 3, scheme = '3';
     
207.            case 4, scheme = '4';
     
208. end
     
209. 
     
210.       %If scheme is 'c', call <nCreate> to create new NNs, train them then perform forecast
     
211.       if(scheme == '1')
     
212.            nCreate;
     
213.       end
     
214. 
     
215.       %If scheme is 'r', call <nRetrain> to retrain the existing NNs
     
216.       if(scheme == '2')
     
217.            nRetrain;
     
218.       end
     
219. 
     
220.       %If scheme is 'f', call <nForecast> to load the existing NN model
     
221.       if(scheme == '3')
     
222.            nForecast;
     
223.       end
     
224. 
     
225.       %If scheme is 'e', verifies and quit session if 'yes' is selected else continue
     
226.       if(scheme == '4')
     
227.            button = questdlg('Quit session?', 'Exit Dialog','Yes','No','No');
     
228.            switch button
     
229.                  case 'Yes', disp(' ');
     
230.                                    disp('Session has ended!!');
     
231.                                    disp(' ');
     
232.                                    break;
     
233.                  case 'No', quit cancel;
     
234.            end
     
235.       end
     
236. end
     
237. 
     
238. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
239. 
     
240. %File name :      ncreate.m
     
241. %Description : This file prepares the input &  output data for the NNs. It creates then trains the
     
242. %                  NNs to mature them.
     
243. 
     
244. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
245. %Clear command screen and set target demand ouput to start at point 2
     
246. clc;
     
247. targetStartAt = 2;
     
248. 
     
249. disp('Program will now CREATE a Neural Network for training and forecasting...');
     
250. disp(' ');
     
251. disp('To capture the pattern of the signal, the model is ')
     
252. disp('set to accept dataPoints x 2 sets of training examples; ');
     
253. disp('1 set of demand + 1 sets of price. ');
     
254. disp(' ');
     
255. disp('The normalised demand data <point 2>, is to be taken as the ')
     
256. disp('output value for the first iteration of training examples. ');
     
257. disp(' ');
     
258. disp('Press ENTER key to continue...');
     
259. pause;
     
260. 
     
261. %Clear command screen then prompt for no. of training cycles
     
262. %For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
     
263. clc;
     
264. cycle = input('Input the number of training cycles: ');
     
265. 
     
266. numOfTimes = resolution + 1;
     
267. %Creating and training the NNs for the respective
     
268. %demand and price coefficient signals
     
269. for x = 1: numOfTimes
     
270. 
     
271.       %Clearing variables
     
272.       clear targetDemand;
     
273.       clear inputs;
     
274.       clear output;
     
275.       clc;
     
276. 
     
277.       if(x == 1)
     
278.             neuralNetwork = ['Neural network settings for approximation level ',
     
279.      num2str(resolution)];
     
280.       else
     
281.             neuralNetwork = ['Neural network settings for detail level ', num2str(x - 1)];
     
282.       end
     
283.       disp(neuralNetwork);
     
284.       disp(' ');
     
285. 
     
286.       %Set no. of input nodes and hidden neurons for the
     
287.       %respective demand and price coefficient signal
     
288.       numOfInputs = 2;
     
289.       inputValue = ['Number of neural network INPUT units is set at ', num2str(numOfInputs)];
     
290.       disp(inputValue);
     
291.       disp(' ');
     
292.       numOfOutput = 1;
     
293.       outValue = ['Output is set to ', num2str(numOfOutput)];
     
294.       disp(outValue);
     
295.       disp(' ');
     
296.       numOfHiddens = input('Enter the no. of HIDDEN units for the NN hidden : ');
     
297.       hiddenValue = ['Number of neural network HIDDEN units is set at ',
     
298.      num2str(numOfHiddens)];
     
299.       disp(hiddenValue);
     
300.       disp(' ');
     
301. 
     
302.      %Setting no. of training examples
     
303.      trainingLength = dataPoints;
     
304. 
     
305.       %Set target outputs of the training examples
     
306.      if(x == 1)
     
307.            targetDemand = normDemand(targetStartAt: 1 + trainingLength);
     
308.      else
     
309.            targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + trainingLength);
     
310.      end
     
311. 
     
312.      %Preparing training examples
     
313.       %Setting training i/ps to be 2 with user defined no. of iterations (dataPoints)
     
314.      y = 0;
     
315.      while y < trainingLength
     
316.            if(x == 1)
     
317.                 inputs(1, y + 1) = normDemand(y + 1);
     
318.                 inputs(2, y + 1) = normPrice(y + 1);
     
319. 
     
320.            else
     
321.                 inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
     
322.                 inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
     
323. 
     
324.            end
     
325. 
     
326.            output(y + 1, :) = targetDemand(y + 1);
     
327. 
     
328.            y = y + 1;
     
329.      end
     
330. 
     
331.      inputs = (inputs');
     
332. 
     
333.      %Setting no. of training cycles
     
334.       [ni, np] = size(targetDemand); % <== [ni, np] tells the NN how long is 1 cycle;
     
335.      size(targetDemand)
     
336. 
     
337.      clear nn;
     
338. 
     
339.       %Create new neural network for respective coefficient signal
     
340.      %NET = MLP(NIN, NHIDDEN, NOUT, FUNC)
     
341.      nn = mlp(numOfInputs, numOfHiddens, numOfOutput, 'linear');
     
342. 
     
343.      %NN options
     
344.       options = zeros(1, 18);
     
345.       options(1) = 1; %Provides display of error values
     
346.       options(14) = cycle * ni * np;
     
347. 
     
348.      %Training the neural network
     
349. 
     
350.       %netopt(net, options, x, t, alg);
     
351.       nn = netopt(nn, options, inputs, output, 'scg');
     
352. 
     
353.       %Save the neural network
     
354.       filename = ['nn', num2str(x)];
     
355.       save(filename, 'nn');
     
356. 
     
357.       disp(' ');
     
358.       msg = ['Neural network successfully CREATED and saved as => ', filename];
     
359.       disp(msg);
     
360. 
     
361.       if(x < 3)
     
362.             disp(' ');
     
363.             disp('Press ENTER key to continue training the next NN...');
     
364.       else
     
365.             disp(' ');
     
366.             disp('Model is now ready to forecast!!');
     
367.             disp(' ');
     
368.             disp('Press ENTER key to continue...');
     
369.       end
     
370. 
     
371.       pause;
     
372. end
     
373. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
374. 
     
375. %File name :       nforecast.m
     
376. %Description : This file loads the trained NNs for load forcasting and %recombines the predicted
     
377. %                   outputs from the NNs to form the final predicted load series.
     
378. 
     
379. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
380. %Clear command screen and variables
     
381. clc;
     
382. clear forecastResult;
     
383. clear actualDemand;
     
384. clear predicted;
     
385. clear actualWithPredicted;
     
386. 
     
387. disp('Neural networks loaded, performing electricity demand forecast...');
     
388. disp(' ');
     
389. pointsAhead = input('Enter no. of points to predict (should be < 336): ');
     
390. 
     
391. numOfTimes = resolution + 1;
     
392. %Predict coefficient signals using respective NNs
     
393. for x = 1 : numOfTimes
     
394.       %Loading NN
     
395.       filename = ['nn', num2str(x)];
     
396.       clear nn
     
397.       load(filename);
     
398.       clear in;
     
399.       numOfInputs = nn.nin;
     
400.       y = 0;
     
401.       %Preparing details to forecast
     
402.       while y < pointsAhead
     
403.            if(x == 1)
     
404.                  propData(1, y + 1) = normDemand(y + 1);
     
405.                  propData(2, y + 1) = normPrice(y + 1);
     
406. 
     
407.            else
     
408.                  propData(1, y + 1) = normDemandDetail(x - 1, y + 1);
     
409.                  propData(2, y + 1) = normPriceDetail(x - 1, y + 1);
     
410. 
     
411.            end
     
412. 
     
413.            y = y + 1;
     
414.       end
     
415. 
     
416.       if(x == 1)
     
417.            forecast = mlpfwd(nn, propData') * maxDemand;
     
418.       else
     
419.             forecast = mlpfwd(nn, propData') * maxDemandDetail(x - 1) - 4000;
     
420.       end
     
421.       forecastResult(x, :) = forecast';
     
422. end
     
423. 
     
424. %For debugging
     
425. % forecastResult
     
426. 
     
427. actualDemand = tDemandArray(2: 1 + pointsAhead);
     
428. predicted = sum(forecastResult, 1)';
     
429. 
     
430. %Calculating Mean Absolute Error
     
431. AbsError = abs(actualDemand - predicted(1: pointsAhead)) ./ actualDemand;
     
432. msg = ['Mean Absolute Error = ', num2str(mean(AbsError(1: pointsAhead))), ' !!'];
     
433. disp(' ');
     
434. disp(msg);
     
435. 
     
436. %Plot actual time series against predicted result
     
437. figure(3)
     
438. actualWithPredicted(:, 1) = actualDemand;
     
439. actualWithPredicted(:, 2) = predicted(1: pointsAhead);
     
440. plot(actualWithPredicted);
     
441. graph = ['Mean Absolute Error = ', num2str(mean(AbsError))];
     
442. title(graph);
     
443. legend('Actual', 'Forecasted');
     
444. 
     
445. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
446. 
     
447. %File name :      nretrain.m
     
448. %Description : This file loads the existing NNs and trains them again.
     
449. 
     
450. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
451. 
     
452. clc;
     
453. %Prompt for the starting point for training
     
454. disp('Program will now RETRAIN the Neural Networks ')
     
455. disp('with the SAME intial data series again...');
     
456. disp(' ');
     
457. disp('To capture the pattern of the signal, the model is ')
     
458. disp('set to accept dataPoints x 2 sets of training examples; ');
     
459. disp('1 set of demand + 1 sets of price. ');
     
460. disp(' ');
     
461. disp('The normalised demand data <point 2>, is to be taken as the ')
     
462. disp('output value for the first iteration of training examples. ');
     
463. disp(' ');
     
464. msg = ['Data points to be used for reTraining the NNs is from 1 to ',...
     
465.             num2str(dataPoints)];
     
466. disp(msg);
     
467. disp(' ');
     
468. disp('Press ENTER key to continue...');
     
469. pause;
     
470. 
     
471. %Clear command screen
     
472. clc;
     
473. 
     
474. %Prompt for no. of training cycles
     
475. %For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
     
476. cycle = input('Input number of cycles to retrain the NNs: ');
     
477. 
     
478. numOfTimes = resolution + 1;
     
479. %Loading existing NNs for training
     
480. for x = 1: numOfTimes
     
481. 
     
482.       %Re-initialising variables
     
483.       clear targetDemand;
     
484.       clear inputs;
     
485.       clear output;
     
486.       clc;
     
487. 
     
488.       %Loading NN for the respective demand and temperature coefficient signals
     
489.       filename = ['nn', num2str(x)];
     
490.       clear nn
     
491.       load(filename);
     
492. 
     
493.       %Getting the size of NN
     
494.       numOfInputs = nn.nin;
     
495.       numOfHiddens = nn.nhidden;
     
496.       numOfOutput = 1;
     
497. 
     
498.       %Setting length of reTraining examples and target outputs
     
499.       reTrainLength = dataPoints;
     
500.       targetLength = reTrainLength;
     
501. 
     
502.       targetStartAt = 2;
     
503. 
     
504.       %Set target outputs of the training examples
     
505.       if(x == 1)
     
506.             targetDemand = normDemand(targetStartAt: 1 + targetLength);
     
507.       else
     
508.             targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + targetLength);
     
509.       end
     
510. 
     
511.       %Preparing training examples
     
512.       %Setting training i/ps to be 2 with user set no. of iterations (dataPoints)
     
513.       y = 0;
     
514.       while y < reTrainLength
     
515.             if(x == 1)
     
516.                   inputs(1, y + 1) = normDemand(y + 1);
     
517.                   inputs(2, y + 1) = normPrice(y + 1);
     
518. 
     
519.             else
     
520.                   inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
     
521.                   inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
     
522. 
     
523.             end
     
524. 
     
525.             output(y + 1, :) = targetDemand(y + 1);
     
526. 
     
527.             y = y + 1;
     
528.       end
     
529. 
     
530.       inputs = (inputs');
     
531. 
     
532.       %Setting no. of training cycles
     
533.       [ni, np] = size(targetDemand); % <== [ni, np] tells the NN how long is 1 cycle;
     
534.       size(targetDemand)                                %With reference to line 106
     
535. 
     
536.       %NN options
     
537.       options = zeros(1, 18);
     
538.       options(1) = 1; %Provides display of error values
     
539.       options(14) = cycle * ni * np;
     
540. 
     
541.       %Training the neural network
     
542.       %netopt(net, options, x, t, alg);
     
543.       nn = netopt(nn, options, inputs, output, 'scg');
     
544. 
     
545.       %Save the neural network
     
546.       filename = ['nn', num2str(x)];
     
547.       save(filename, 'nn');
     
548. 
     
549.       disp(' ');
     
550.       msg = ['Neural network => ', filename, ' <= successfully RETRAINED and saved!! '];
     
551. 
     
552.       disp(msg);
     
553. 
     
554.       if(x < 3)
     
555.             disp(' ');
     
556.             disp('Press ENTER key to continue training the next NN...');
     
557.       else
     
558.             disp(' ');
     
559.             disp('Model is now ready to forecast again!!');
     
560.             disp(' ');
     
561.             disp('Press ENTER key to continue...');
     
562.       end
     
563. 
     
564.       pause;
     
565. end

复制代码

simon21

1. function feature_wl=getf_wavelet(Input)
   
2.    
   
3.     A=double(Input);
   
4.     Xrgb=0.2990*A(:,:,1)+0.5870*A(:,:,2)+0.1140*A(:,:,3);
   
5.     NbColor=255;
   
6. %  WCODEMAT Extended pseudocolor matrix scaling. Y = WCODEMAT(X,NBCODES,OPT,ABSOL) returns a coded version of input
   
7. % matrix X if ABSOL=0, or ABS(X) if ABSOL is  nonzero, using the first NBCODES integers.
   
8. %  Coding can be done row-wise (OPT='row' or 'r'), columnwise (OPT='col' or 'c'), or globally (OPT='mat' or 'm').
   
9. %  Coding uses a regular grid between the minimum and the maximum values of each row (column or matrix,respectively).
   
10. %  Y = WCODEMAT(X,NBCODES) is equivalent to  Y = WCODEMAT(X,NBCODES,'mat',1).
    
11.     X1=wcodemat(Xrgb,NbColor);%伪彩矩阵压缩
    
12. 
    
13. % Function 'wavedec2' Multilevel 2-D wavelet decomposition.[C,S] = WAVEDEC2(X,N,'wname') returns the wavelet
    
14. % decomposition of the matrix X at level N, using the wavelet named in string 'wname' (see WFILTERS).
    
15. % Outputs are the decomposition vector C and the corresponding bookkeeping matrix S.
    
16. % N must be a strictly positive integer (see WMAXLEV).
    
17.    [c2,l2]=wavedec2(X1,4,'sym4');%二维小波分解
    
18. 
    
19. % Function 'detcoef2' Extract 2-D detail coefficients.
    
20. % D = detcoef2(O,C,S,N) extracts from the wavelet decomposition structure [C,S], the horizontal, vertical
    
21. % or diagonal detail coefficients for O = 'h' (or 'v' or 'd',respectively), at level N. N must
    
22. % be an integer such that 1 <= N <= size(S,1)-2.
    
23.    ch11=detcoef2('h',c2,l2,1);%提取小波分解细节系数
    
24.    ch12=detcoef2('h',c2,l2,2);
    
25.    ch13=detcoef2('h',c2,l2,3);
    
26.    ch14=detcoef2('h',c2,l2,4);
    
27.    cv11=detcoef2('v',c2,l2,1);
    
28.    cv12=detcoef2('v',c2,l2,2);
    
29.    cv13=detcoef2('v',c2,l2,3);
    
30.    cv14=detcoef2('v',c2,l2,4);
    
31.    cd11=detcoef2('d',c2,l2,1);
    
32.    cd12=detcoef2('d',c2,l2,2);
    
33.    cd13=detcoef2('d',c2,l2,3);
    
34.    cd14=detcoef2('d',c2,l2,4);
    
35. 
    
36. % Function 'appcoef2' Extract 2-D approximation coefficients.APPCOEF2 computes the approximation coefficients of a
    
37. % two-dimensional signal.A = APPCOEF2(C,S,'wname',N) computes the approximation coefficients at level N using the wavelet decomposition
    
38. % structure [C,S] (see WAVEDEC2). 'wname' is a string containing the wavelet name.Level N must be an integer such that 0 <= N <= size(S,1)-2.
    
39.    ca14=appcoef2(c2,l2,'sym4',4);%提取小波分解概貌系数
    
40. 
    
41. %==========compute mean & var(求分解后的每个图象的均值和方差)==================
    
42. %image x1's u & std
    
43.    Csize1=size(ch11);
    
44.    s1=0;
    
45.    for i=1:Csize1(1)
    
46.        for j=1:Csize1(2)
    
47.            s1=s1+abs(ch11(i,j));%abs求绝对值
    
48.        end
    
49.    end
    
50.    u1=(1/(Csize1(1)*Csize1(2)))*s1;
    
51. 
    
52.    st1=0;
    
53.    for i=1:Csize1(1)
    
54.        for j=1:Csize1(2)
    
55.            st1=st1+(abs(ch11(i,j))-u1)*(abs(ch11(i,j))-u1);
    
56.        end
    
57.    end
    
58.    z1=sqrt((1/(Csize1(1)*Csize1(2)-1))*st1);
    
59. %======================================================
    
60.    Csize2=size(ch12);
    
61.    s2=0;
    
62.    for i=1:Csize2(1)
    
63.        for j=1:Csize2(2)
    
64.            s2=s2+abs(ch12(i,j));
    
65.        end
    
66.    end
    
67.    u2=(1/(Csize2(1)*Csize2(2)))*s2;
    
68. 
    
69.    st2=0;
    
70.    for i=1:Csize2(1)
    
71.        for j=1:Csize2(2)
    
72.            st2=st2+(abs(ch12(i,j))-u2)*(abs(ch12(i,j))-u2);
    
73.        end
    
74.    end
    
75.    z2=sqrt((1/(Csize2(1)*Csize2(2)-1))*st2);
    
76. %===========================================================
    
77.    Csize3=size(ch13);
    
78.    s3=0;
    
79.    for i=1:Csize3(1)
    
80.        for j=1:Csize3(2)
    
81.            s3=s3+abs(ch13(i,j));
    
82.        end
    
83.    end
    
84.    u3=(1/(Csize3(1)*Csize3(2)))*s3;
    
85. 
    
86.    st3=0;
    
87.    for i=1:Csize3(1)
    
88.        for j=1:Csize3(2)
    
89.            st3=st3+(abs(ch13(i,j))-u3)*(abs(ch13(i,j))-u3);
    
90.        end
    
91.    end
    
92.    z3=sqrt((1/(Csize3(1)*Csize3(2)-1))*st3);
    
93. %================================================================
    
94.    Csize3=size(ch14);
    
95.    s3=0;
    
96.    for i=1:Csize3(1)
    
97.        for j=1:Csize3(2)
    
98.            s3=s3+abs(ch14(i,j));
    
99.        end
    
100.    end
     
101.    u4=(1/(Csize3(1)*Csize3(2)))*s3;
     
102. 
     
103.    st3=0;
     
104.    for i=1:Csize3(1)
     
105.        for j=1:Csize3(2)
     
106.            st3=st3+(abs(ch14(i,j))-u4)*(abs(ch14(i,j))-u4);
     
107.        end
     
108.    end
     
109.    z4=sqrt((1/(Csize3(1)*Csize3(2)-1))*st3);
     
110. %=================================================================
     
111.    Csize4=size(cv11);
     
112.    s4=0;
     
113.    for i=1:Csize4(1)
     
114.        for j=1:Csize4(2)
     
115.           s4=s4+abs(cv11(i,j));
     
116.        end
     
117.    end
     
118.    u5=(1/(Csize4(1)*Csize4(2)))*s4;
     
119. 
     
120.    st4=0;
     
121.    for i=1:Csize4(1)
     
122.        for j=1:Csize4(2)
     
123.            st4=st4+(abs(cv11(i,j))-u5)*(abs(cv11(i,j))-u5);
     
124.        end
     
125.    end
     
126.    z5=sqrt((1/(Csize4(1)*Csize4(2)-1))*st4);
     
127. %===============================================================
     
128.    Csize5=size(cv12);
     
129.    s5=0;
     
130.    for i=1:Csize5(1)
     
131.        for j=1:Csize5(2)
     
132.            s5=s5+abs(cv12(i,j));
     
133.        end
     
134.    end
     
135.    u6=(1/(Csize5(1)*Csize5(2)))*s5;
     
136. 
     
137.    st5=0;
     
138.    for i=1:Csize5(1)
     
139.        for j=1:Csize5(2)
     
140.            st5=st5+(abs(cv12(i,j))-u6)*(abs(cv12(i,j))-u6);
     
141.        end
     
142.    end
     
143.    z6=sqrt((1/(Csize5(1)*Csize5(2)-1))*st5);
     
144. %================================================================
     
145.    Csize6=size(cv13);
     
146.    s6=0;
     
147.    for i=1:Csize6(1)
     
148.        for j=1:Csize6(2)
     
149.            s6=s6+abs(cv13(i,j));
     
150.        end
     
151.    end
     
152.    u7=(1/(Csize6(1)*Csize6(2)))*s6;
     
153. 
     
154.    st6=0;
     
155.    for i=1:Csize6(1)
     
156.        for j=1:Csize6(2)
     
157.           st6=st6+(abs(cv13(i,j))-u7)*(abs(cv13(i,j))-u7);
     
158.        end
     
159.    end
     
160.    z7=sqrt((1/(Csize6(1)*Csize6(2)-1))*st6);
     
161. %==================================================================
     
162.    Csize6=size(cv14);
     
163.    s6=0;
     
164.    for i=1:Csize6(1)
     
165.        for j=1:Csize6(2)
     
166.            s6=s6+abs(cv14(i,j));
     
167.        end
     
168.    end
     
169.    u8=(1/(Csize6(1)*Csize6(2)))*s6;
     
170. 
     
171.    st6=0;
     
172.    for i=1:Csize6(1)
     
173.        for j=1:Csize6(2)
     
174.            st6=st6+(abs(cv14(i,j))-u8)*(abs(cv14(i,j))-u8);
     
175.        end
     
176.    end
     
177.    z8=sqrt((1/(Csize6(1)*Csize6(2)-1))*st6);
     
178. %==============================================================
     
179.    Csize7=size(cd11);
     
180.    s7=0;
     
181.    for i=1:Csize7(1)
     
182.        for j=1:Csize7(2)
     
183.           s7=s7+abs(cd11(i,j));
     
184.        end
     
185.    end
     
186.    u9=(1/(Csize7(1)*Csize7(2)))*s7;
     
187. 
     
188.    st7=0;
     
189.    for i=1:Csize7(1)
     
190.        for j=1:Csize7(2)
     
191.           st7=st7+(abs(cd11(i,j))-u9)*(abs(cd11(i,j))-u9);
     
192.        end
     
193.    end
     
194.    z9=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
     
195. %====================================================================
     
196.    Csize7=size(cd12);
     
197.    s7=0;
     
198.    for i=1:Csize7(1)
     
199.        for j=1:Csize7(2)
     
200.           s7=s7+abs(cd12(i,j));
     
201.        end
     
202.    end
     
203.    u10=(1/(Csize7(1)*Csize7(2)))*s7;
     
204. 
     
205.    st7=0;
     
206.    for i=1:Csize7(1)
     
207.       for j=1:Csize7(2)
     
208.           st7=st7+(abs(cd12(i,j))-u10)*(abs(cd12(i,j))-u10);
     
209.        end
     
210.    end
     
211.    z10=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
     
212. %=================================================
     
213.    Csize7=size(cd13);
     
214.    s7=0;
     
215.    for i=1:Csize7(1)
     
216.       for j=1:Csize7(2)
     
217.           s7=s7+abs(cd13(i,j));
     
218.        end
     
219.    end
     
220.    u11=(1/(Csize7(1)*Csize7(2)))*s7;
     
221. 
     
222.    st7=0;
     
223.    for i=1:Csize7(1)
     
224.       for j=1:Csize7(2)
     
225.          st7=st7+(abs(cd13(i,j))-u11)*(abs(cd13(i,j))-u11);
     
226.       end
     
227.    end
     
228.    z11=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
     
229. %=================================================
     
230.    Csize7=size(cd14);
     
231.    s7=0;
     
232.    for i=1:Csize7(1)
     
233.        for j=1:Csize7(2)
     
234.           s7=s7+abs(cd14(i,j));
     
235.        end
     
236.    end
     
237.    u12=(1/(Csize7(1)*Csize7(2)))*s7;
     
238. 
     
239.    st7=0;
     
240.    for i=1:Csize7(1)
     
241.        for j=1:Csize7(2)
     
242.           st7=st7+(abs(cd14(i,j))-u12)*(abs(cd14(i,j))-u12);
     
243.        end
     
244.    end
     
245.    z12=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
     
246. %=================================================
     
247.    Csize7=size(ca14);
     
248.    s7=0;
     
249.    for i=1:Csize7(1)
     
250.       for j=1:Csize7(2)
     
251.           s7=s7+abs(ca14(i,j));
     
252.       end
     
253.    end
     
254.    u13=(1/(Csize7(1)*Csize7(2)))*s7;
     
255. 
     
256.    st7=0;
     
257.    for i=1:Csize7(1)
     
258.       for j=1:Csize7(2)
     
259.          st7=st7+(abs(ca14(i,j))-u13)*(abs(ca14(i,j))-u13);
     
260.       end
     
261.    end
     
262.    z13=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
     
263. 
     
264.    %semble feature vector========================================
     
265.    % 应用准则函数把对分类影响不大的特征去除
     
266.    feature_wl=[u1,z1,u2,z2,u3,z3,u4,z4,u5,z5,u6,z6,u7,z7,u8,z8,u9,z9,u10,z10,u11,z11,u12,z12,u13,z13];
     
267.    %feature_wl=[u1,z1,u2,z2,u3,z3,u4,z4,u5,z5,u6,z6,u7,z7,u9,z9,u13,z13];

复制代码

orchis_2005

请问有没有关于时频分析
       的chirplet变换matlab程序呢？谢谢！

gaojizhen

请问楼主:
   我刚刚接触matlab和小波变换,   [c2,l2]=wavedec2(X1,4,'sym4');中对其参数有什么要求吗?为什么每次执行到这条语句的时候就出问题啊?

zmmwc

楼主真慷慨，谢谢啦！请问有没有用小波神经网络进行图象识别的程序？

千芮

有没有关于说话人识别中利用小波变换提取特征参数的程序，不胜感激，邮箱zbl@mail2.lut.cn

点点

LZ有关于小波奇异性检测方面的源程序么？
谢谢了:loveliness:
真是个慷慨的好LZ啊！