二维小波变换（正和逆变换）
二维小波变换（正和逆变换），用二维卷积实现的，其中有个子函数用于产生小波系数的（小波分解和小波重构的系数都有）！你自己也可以往里面加系数！function =wavefilter(wname,type)
error(nargchk(1,2,nargin));
if(nargin==1 & nargout ~=4) | (nargin==2 & nargout~=2)
error('Invalid number of output arguments.');
end
if nargin==1 & ~ischar(wname)
error('WNAME must be a string.');
end
if nargin==2 & ~ischar(type)
error('TYPE must be a string.');
end
switch lower(wname)
case {'haar','db1'}
ld=/sqrt(2); hd=[-1 1]/sqrt(2);
lr=ld; hr=-hd;
case 'db4'
ld=[-1.059740178499728e-002 3.288301166698295e-002 3.084138183598697e-002 -1.870348117188811e-001 ...
-2.798376941698385e-002 6.308807679295904e-001 7.148465705525415e-001 2.303778133088552e-001];
t=;
hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
lr=ld; lr(end:-1:1)=ld;
hr=cos(pi*t).*ld;
case 'sym4'
ld=[-7.576571478927333e-002 -2.963552764599851e-002 4.976186676320155e-001 8.037387518059161e-001 ...
2.978577956052774e-001-9.921954357684722e-002 -1.260396726203783e-002 3.222310060404270e-002];
t=;
hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
lr=ld; lr(end:-1:1)=ld;
hr=cos(pi*t).*ld;
case 'bior6.8'
ld=[0 1.908831736481291e-003 -1.9142861290088767e-003 -1.699063986760234e-002 1.1934565279772926e-002 ...
4.973290349094079e-002 -7.726317316720414e-002 -9.405920349573646e-002 4.207962846098268e-001 ...
8.259229974584023e-001 4.207962846098268e-001 -9.405920349573646e-002 -7.726317316720414e-002 ...
4.973290349094079e-002 1.193456527972926e-002 -1.699063986760234e-002 -1.914286129088767e-003 ...
1.908831736481291e-003];
hd=[0 0 0 1.442628250562444e-002 -1.446750489679015e-002 -7.872200106262882e-002 4.036797903033992e-002 ...
4.178491091502746e-001 -7.589077294536542e-001 4.178491091502746e-001 4.036797903033992e-002 ...
-7.872200106262882e-002 -1.446750489679015e-002 1.442628250562444e-002 0 0 0 0];
t=;
lr=cos(pi*(t+1)).*hd;
hr=cos(pi*t).*ld;
case 'jpeg9.7'
ld=[0 0.02674875741080976 -0.01686411844287495 -0.07822326652898785 0.2668641184428723 ...
0.6029490182363579 0.2668641184428723 -0.07822326652898785 -0.01686411844287495...
0.02674875741080976];
hd=[0 -0.09127176311424948 0.05754352622849957 0.5912717631142470 -1.115087052456994 ...
0.5912717631142470 0.05754352622849957 -0.09127176311424948 0 0];
t=;
lr=cos(pi*(t+1)).*hd;
hr=cos(pi*t).*ld;
otherwise
error('Unrecongizable wavelet name (WNAME).');
end
if(nargin==1)
varargout(1:4)={ld,hd,lr,hr};
else
switch lower(type(1))
case 'd'
varargout={ld,hd};
case 'r'
varargout={lr,hr};
otherwise
error('Unrecongizable filter TYPE.');
end
end
function =wavefast(x,n,varargin)
error(nargchk(3,4,nargin));
if nargin==3
if ischar(varargin{1})
=wavefilter(varargin{1},'d');
else
error('Missing wavelet name.');
end
else
lp=varargin{1}; hp=varargin{2};
end
fl=length(lp);sx=size(x);
if (ndims(x)~=2) | (min(sx)<2) | ~isreal(x) | ~isnumeric(x)
error('X must be a real ,numeric matric.');
end
if(ndims(lp)~=2) | ~ isreal(lp) | ~isnumeric(lp) | (ndims(hp) ~=2) | ~ isreal(hp) | ~isnumeric(hp) ...
| (fl~=length(hp)) | rem(fl,2)~=0
error(['LP and HP must be ever and equal length real numeric filter vector.']);
end
if ~isreal(n) | ~isnumeric(n) |(n<1) |(n>log2(max(sx)))
error(['N must be a real scalar between 1 and log2(max(size(X))).']);
end
c=[];s=sx;app=double(x);
for i=1:n
=symextend(app,fl);
rows=symconv(app,hp,'row',fl,keep);
coefs=symconv(rows,hp,'col',fl,keep);
c=; s=;
coefs=symconv(rows,lp,'col',fl,keep);
c=;
rows=symconv(app,lp,'row',fl,keep);
coefs=symconv(rows,hp,'col',fl,keep);
c=;
app=symconv(rows,lp,'col',fl,keep);
end
c=; s=;
function =symextend(x,fl)
keep=floor((fl+size(x)-1)/2);
y=padarray(x,[(fl-1) (fl-1)],'symmetric','both');
function y=symconv(x,h,type,fl,keep)
if strcmp(type, 'row')
y=conv2(x,h);
y=y(:,1:2:end);
y=y(:,fl/2+1:fl/2+keep(2));
else
y=conv2(x,h');
y=y(1:2:end,:);
y=y(fl/2+1:fl/2+keep(2),:);
end
function =waveback(c,s,varargin)
error(nargchk(3,5,nargin));
error(nargchk(1,2,nargout));
if(ndims(c) ~=2) | (size(c,1) ~=1)
error('C must be a row vector .');
end
if (ndims(s) ~=2) | ~isreal(s) | ~isnumeric(s) | (size(s,2) ~=2)
error('S must be a real,numeric two-column array.');
end
elements=prod(s,2);
if (length(c) <elements(end)) | ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))
error([' must be a standard wavelet decomposition structure.']);
end
nmax=size(s,1)-2;
wname=varargin{1};filterchk=0;nchk=0;
switch nargin
case 3
if ischar(wname)
=wavefilter(wname,'r');
n=nmax;
else
error('Undefined filter.');
end
if nargout~=1
error('Wrong number of output arguments.');
end
case 4
if ischar(wname)
=wavefilter(wname,'r');
n=varargin{2};nchk=1;
else
lp=varargin{1};
hp=varargin{2};
filterchk=1;n=nmax;
if nargout ~=1
error('Wrong number of output arguments.');
end
end
case 5
lp=varargin{1};hp=varargin{2};filterchk=1;
n=varargin{3};nchk=1;
otherwise
error('Improper number of input arguments.');
end
fl=length(lp);
if filterchk
if (ndims(lp) ~=2) | ~isreal(lp) | ~isnueric(lp) |(ndims(hp) ~=2) | ~isreal(hp) ...
| ~isnumeric(hp) |(fl ~=length(hp)) | rem(fl,2) ~=0
error(['LP and HP must be even and equal length real,numeric filter vectors.']);
end
end
if nchk & (~isnumeric(n) | ~isreal(n))
error('N must be a real numeric.');
end
if(n~=nmax) & (nargout ~=2)
error('Not enough output arguments.');
end
nc=c;ns=s;nnmax=nmax;
for i=1:n
a=symconvup(wavecopy('a',nc,ns),lp,lp,fl,ns(3,:))+ ...
symconvup(wavecopy('h',nc,ns,nnmax),hp,lp,fl,ns(3,:))+ ...
symconvup(wavecopy('v',nc,ns,nnmax),lp,hp,fl,ns(3,:))+ ...
symconvup(wavecopy('d',nc,ns,nnmax),hp,hp,fl,ns(3,:));
nc=nc(4*prod(ns(1,:))+1:end);
nc=;
ns=ns(3:end,:);
ns=;
nnmax=size(ns,1)-2;
end
if nargout ==1
a=nc; nc=repmat(0,ns(1,:)); nc(:)=a;
end
varargout{1}=nc;
if nargout==2
varargout{2}=ns;
end
function z=symconvup(x,f1,f2,fln,keep)
y=zeros(.*size(x));
y(1:2:end,:)=x;
y=conv2(y,f1');
z=zeros(.*size(y));
z(:,1:2:end)=y;
z=conv2(z,f2);
z=z(fln-1:fln+keep(1)-2,fln-1:fln+keep(2)-2);
第二代小波变换源码
%%第二代小波变换源码%%此程序用提升法实现第二代小波变换
%%我用的是非整数阶小波变换
%%采用时域实现,步骤先列后行
%%正变换：分裂，预测，更新；
%%反变换：更新，预测，合并
%%只做一层（可以多层，而且每层的预测和更新方程不同）
clear;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%1.调原始图像矩阵
load wbarb;%下载图像
f=X; %原始图像
% f=[0 0 0 0 0 0 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 0 0 0 0 0 0 ;];%原始图像矩阵
N=length(f); %图像维数
T=N/2; %子图像维数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%正变换%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 1.列变换
%A.分裂（奇偶分开）
f1=f(,:);%奇数
f2=f(,:); %偶数
% f1(:,T+1)=f1(:,1);%补列
% f2(T+1,:)=f2(1,:);%补行
%B.预测
for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
end;
% high_frequency_column(T+1,:)=high_frequency_column(1,:);%补行
%C.更新
for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
end;
%D.合并
f_column(,:)=low_frequency_column(,:);
f_column(,:)=high_frequency_column(,:);
figure(1)
colormap(map);
image(f);
figure(2)
colormap(map);
image(f_column);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.行变换
%A.分裂（奇偶分开）
f1=f_column(:,);%奇数
f2=f_column(:,); %偶数
% f2(:,T+1)=f2(:,1); %补行
%B.预测
for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
end;
% high_frequency_row(:,T+1)=high_frequency_row(:,1);%补行
%C.更新
for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
end;
%D.合并
f_row(:,)=low_frequency_row(:,);
f_row(:,)=high_frequency_row(:,);
figure(3)
colormap(map);
image(f_row);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 1.行变换
% A.提取（低频高频分开）
f1=f_row(:,);%奇数
f2=f_row(:,); %偶数
% f2(:,T+1)=f2(:,1); %补行
%B.更新
for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
end;
%C.预测
for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
end;
% high_frequency_row(:,T+1)=high_frequency_row(:,1);%补行
%D.合并(奇偶分开合并)
f_row(:,)=low_frequency_row(:,);
f_row(:,)=high_frequency_row(:,);
figure(4)
colormap(map);
image(f_row);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.列变换
%A.提取（低频高频分开）
f1=f_row(,:);%奇数
f2=f_row(,:); %偶数
% f1(:,T+1)=f1(:,1);%补列
% f2(T+1,:)=f2(1,:);%补行
%B.更新
for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
end;
%C.预测
for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
end;
% high_frequency_column(T+1,:)=high_frequency_column(1,:);%补行
%D.合并(奇偶分开合并）
f_column(,:)=low_frequency_column(,:);
f_column(,:)=high_frequency_column(,:);
figure(5)
colormap(map);
image(f_column);
利用小波变换实现对电能质量检测的算法实现
N=10000;s=zeros(1,N);
for n=1:N
if n<0.4*N||n>0.8*N
s(n)=31.1*sin(2*pi*50/10000*n);
else
s(n)=22.5*sin(2*pi*50/10000*n);
end
end
l=length(s);
=wavedec(s,6,'db5'); %用db5小波分解信号到第六层
subplot(8,1,1);
plot(s);
title('用db5小波分解六层：s=a6+d6+d5+d4+d3+d2+d1');
Ylabel('s');%对分解结构【c，l】中第六层低频部分进行重构
a6=wrcoef('a',c,l,'db5',6);
subplot(8,1,2);
plot(a6);
Ylabel('a6');%对分解结构【c，l】中各层高频部分进行重构
for i=1:6
decmp=wrcoef('d',c,l,'db5',7-i);
subplot(8,1,i+2);
plot(decmp);
Ylabel(['d',num2str(7-i)]);
end
%-----------------------------------------------------------
rec=zeros(1,300);
rect=zeros(1,300);
ke=1;
u=0;
d1=wrcoef('d',c,l,'db5',1);
figure(2);
plot(d1);
si=0;
N1=0;
N0=0;
sce=0;
for n=20:N-30
rect(ke)=s(n);
ke=ke+1;
if(ke>=301)
if(si==2)
rec=rect;
u=2;
end;
si=0;
ke=1;
end;
if(d1(n)>0.01) % the condition of abnormal append.
N1=n;
if(N0==0)
N0=n;
si=si+1;
end;
if(N1>N0+30)
Nlen=N1-N0;
Tab=Nlen/10000;
end;
end;
if(si==1)
for k=N0:N0+99 %testing of 1/4 period signals to
sce=sce+s(k)*s(k)/10000;
end;
re=sqrt(sce*200) %re indicate the pike value of .
sce=0;
si=si+1;
end;
end;
Nlen
N0
n=1:300;
figure(3)
plot(n,rec);
基于小波变换的图象去噪 Normalshrink算法
function =threshold_2_N(img,levels)% reference :image denoising using wavelet thresholding
=size(img);
HH=img((xx/2+1):xx,(yy/2+1):yy);
delt_2=(std(HH(:)))^2;%(median(abs(HH(:)))/0.6745)^2;%
T_img=img;
for i=1:levels
temp_x=xx/2^i;
temp_y=yy/2^i;
% belt=1.0*(log(temp_x/(2*levels)))^0.5;
belt=1.0*(log(temp_x/(2*levels)))^0.5;%2.5 0.8
%HL
HL=img(1:temp_x,(temp_y+1):2*temp_y);
delt_y=std(HL(:));
T_1=belt*delt_2/delt_y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_HL=sign(HL).*max(0,abs(HL)-T_1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_img(1:temp_x,(temp_y+1):2*temp_y)=T_HL;
Sub_T(3*(i-1)+1)=T_1;
%LH
LH=img((temp_x+1):2*temp_x,1:temp_y);
delt_y=std(LH(:));
T_2=belt*delt_2/delt_y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_LH=sign(LH).*max(0,abs(LH)-T_2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_img((temp_x+1):2*temp_x,1:temp_y)=T_LH;
Sub_T(3*(i-1)+2)=T_2;
%HH
HH=img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y);
delt_y=std(HH(:));
T_3=belt*delt_2/delt_y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_HH=sign(HH).*max(0,abs(HH)-T_3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y)=T_HH;
Sub_T(3*(i-1)+3)=T_3;
end
基于小波变换模极大的多尺度图像边缘检测
clear all;load wbarb;
I = ind2gray(X,map);imshow(I);
I1 = imadjust(I,stretchlim(I),);figure;imshow(I1);
= size(I);
h = ;
g = ;
delta = ;
J = 3;
a(1:N,1:M,1,1:J+1) = 0;
dx(1:N,1:M,1,1:J+1) = 0;
dy(1:N,1:M,1,1:J+1) = 0;
d(1:N,1:M,1,1:J+1) = 0;
a(:,:,1,1) = conv2(h,h,I,'same');
dx(:,:,1,1) = conv2(delta,g,I,'same');
dy(:,:,1,1) = conv2(g,delta,I,'same');
x = dx(:,:,1,1);
y = dy(:,:,1,1);
d(:,:,1,1) = sqrt(x.^2+y.^2);
I1 = imadjust(d(:,:,1,1),stretchlim(d(:,:,1,1)),);figure;imshow(I1);
lh = length(h);
lg = length(g);
for j = 1:J+1
lhj = 2^j*(lh-1)+1;
lgj = 2^j*(lg-1)+1;
hj(1:lhj)=0;
gj(1:lgj)=0;
for n = 1:lh
hj(2^j*(n-1)+1)=h(n);
end
for n = 1:lg
gj(2^j*(n-1)+1)=g(n);
end
a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same');
dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same');
x = dx(:,:,1,j+1);
y = dy(:,:,1,j+1);
dj(:,:,1,j+1) = sqrt(x.^2+y.^2);
I1 = imadjust(dj(:,:,1,j+1),stretchlim(dj(:,:,1,j+1)),);figure;imshow(I1);
end
利用小波变根据二进制数（水印）来改变图片，提取其中一个子带的直方图
close allclear all
clc;
Original_image = imread('lena512_marked.bmp');
Original_image = double(Original_image);
=ADWT2(Original_image,1);
HL1 = v1{1};
HL1_1D = sort(HL1(:)).';
bin_size = 2;
range_value = 81;
step =[-range_value:bin_size:range_value];
Extract_histogram1 = histc(HL1_1D,step);% 提取的 histogram 步长为1
% Extract_histogram1 = histc(HL1_1D,);
= size(Extract_histogram1);
relation1 = zeros(1,n1);
for i = 2:3:(n1/2)
relation1(1,i) = 2*Extract_histogram1(1,i)/(Extract_histogram1(1,i+1)+Extract_histogram1(1,i-1));
end
for i = ((n1/2)+4):3:(n1-2)
relation1(1,i) = 2*Extract_histogram1(1,i)/(Extract_histogram1(1,i+1)+Extract_histogram1(1,i-1));
end
%%% Extract watermark W
W_number1 = 24; % 水印的个数 W_number
H1 = zeros(1,3*W_number1);
H1(1:(3*W_number1/2)) = Extract_histogram1(1:(3*W_number1/2));% H的第一个直方图的起点是-range_value
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
a = zeros(1,W_number1);
b = zeros(1,W_number1);
c = zeros(1,W_number1);
H_matrix1 = reshape(H1,3,W_number1);
a = H_matrix1(1,:);
b = H_matrix1(2,:);
c = H_matrix1(3,:); % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
extract_W1=zeros(1,W_number1);
for i=1:W_number1
if (2*b(1,i) >= a(1,i)+c(1,i))
extract_W1(1,i)=1;
else if (2*b(1,i) < (a(1,i)+c(1,i)))
extract_W1(1,i)=-1;
end
end
end
extract_watermark1=extract_W1
%%% 缩小
Original_image = imread('lena512_marked.bmp');
% scaled_image = imresize(Original_image,);
scaled_image = imresize(Original_image,.9);
scaled_image = double(scaled_image);
=ADWT2(scaled_image,1);
HL2 = v2{1};
HL2_1D = sort(HL2(:)).'; % 将HL子带变为一维向量
Extract_histogram2 = histc(HL2_1D,step); % 提取的 histogram 步长为1
= size(Extract_histogram2);
% hist(sort(HL(:)),n);
relation2 = zeros(1,n1);
for i = 2:3:(n1/2)
relation2(1,i) = 2*Extract_histogram2(1,i)/(Extract_histogram2(1,i+1)+Extract_histogram2(1,i-1));
end
for i = ((n1/2)+4):3:(n1-2)
relation2(1,i) = 2*Extract_histogram2(1,i)/(Extract_histogram2(1,i+1)+Extract_histogram2(1,i-1));
end
%%% Extract watermark W
W_number2 = 24; % 水印的个数 W_number
H2 = zeros(1,3*W_number2);
H2(1:(3*W_number2/2)) = Extract_histogram2(1:(3*W_number2/2));% H的第一个直方图的起点是-range_value
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
a = zeros(1,W_number2);
b = zeros(1,W_number2);
c = zeros(1,W_number2);
H_matrix2 = reshape(H2,3,W_number2);
a = H_matrix2(1,:);
b = H_matrix2(2,:);
c = H_matrix2(3,:); % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
extract_W2=zeros(1,W_number2);
for i=1:W_number2
if (2*b(1,i) >= a(1,i)+c(1,i))
extract_W2(1,i)=1;
else if (2*b(1,i) < (a(1,i)+c(1,i)))
extract_W2(1,i)=-1;
end
end
end
extract_watermark2=extract_W2
%%% 放大
Original_image = imread('lena512_marked.bmp');
scaled_image = imresize(Original_image,1.3);
scaled_image = double(scaled_image);
=ADWT2(scaled_image,1);
HL3 = v3{1};
HL3_1D = sort(HL3(:)).'; % 将HL子带变为一维向量
Extract_histogram3 = histc(HL3_1D,step); % 提取的 histogram 步长为1
= size(Extract_histogram3);
% hist(sort(HL(:)),n);
relation3 = zeros(1,n1);
for i = 2:3:(n1/2)
relation3(1,i) = 2*Extract_histogram3(1,i)/(Extract_histogram3(1,i+1)+Extract_histogram3(1,i-1));
end
for i = ((n1/2)+4):3:(n1-2)
relation3(1,i) = 2*Extract_histogram3(1,i)/(Extract_histogram3(1,i+1)+Extract_histogram3(1,i-1));
end
%%% Extract watermark W
W_number3 = 24; % 水印的个数 W_number
H3 = zeros(1,3*W_number3);
H3(1:(3*W_number3/2)) = Extract_histogram3(1:(3*W_number3/2));% H的第一个直方图的起点是-range_value
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
a = zeros(1,W_number3);
b = zeros(1,W_number3);
c = zeros(1,W_number3);
H_matrix3 = reshape(H3,3,W_number3);
a = H_matrix3(1,:);
b = H_matrix3(2,:);
c = H_matrix3(3,:); % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
extract_W3=zeros(1,W_number3);
for i=1:W_number3
if (2*b(1,i) >= a(1,i)+c(1,i))
extract_W3(1,i)=1;
else if (2*b(1,i) < (a(1,i)+c(1,i)))
extract_W3(1,i)=-1;
end
end
end
extract_watermark3=extract_W3
%%% 旋转
Original_image = imread('lena512_marked.bmp');
Rescale_image = imrotate(Original_image,10);
Rescale_image = double(Rescale_image);
=ADWT2(Rescale_image,1);
HL4 = v4{1};
HL4_1D = sort(HL4(:)).'; % 将HL子带变为一维向量
Extract_histogram4 = histc(HL4_1D,step); % 提取的 histogram 步长为1
= size(Extract_histogram4);
% hist(sort(HL(:)),n);
relation4 = zeros(1,n1);
for i = 2:3:(n1/2)
relation4(1,i) = 2*Extract_histogram4(1,i)/(Extract_histogram4(1,i+1)+Extract_histogram4(1,i-1));
end
for i = ((n1/2)+4):3:(n1-2)
relation4(1,i) = 2*Extract_histogram4(1,i)/(Extract_histogram4(1,i+1)+Extract_histogram4(1,i-1));
end
%%% Extract watermark W
W_number4 = 24; % 水印的个数 W_number4
H4 = zeros(1,3*W_number4);
H4(1:(3*W_number4/2)) = Extract_histogram4(1:(3*W_number4/2));% H的第一个直方图的起点是-range_value
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
a = zeros(1,W_number4);
b = zeros(1,W_number4);
c = zeros(1,W_number4);
H_matrix4 = reshape(H4,3,W_number4);
a = H_matrix4(1,:);
b = H_matrix4(2,:);
c = H_matrix4(3,:); % a(i),b(i) and c(i) 是嵌入第i个水印的三个bin的大小
extract_W4=zeros(1,W_number4);
for i=1:W_number4
if (2*b(1,i) >= a(1,i)+c(1,i))
extract_W4(1,i)=1;
else if (2*b(1,i) < (a(1,i)+c(1,i)))
extract_W4(1,i)=-1;
end
end
end
extract_watermark4=extract_W4
X1=2:1:(n1-2);
Y1=relation1(1,X1);
X2=2:1:(n1-2);
Y2=relation2(1,X2);
X3=2:1:(n1-2);
Y3=relation3(1,X3);
X4=2:1:(n1-2);
Y4=relation4(1,X4);
figure;plot(X1,Y1,'r*',X2,Y2,'b+',X3,Y3,'kd',X4,Y4,'gd');
title('The estimation of bin relation');
用小波函数构建神经网络的源程序
该程序是用小波函数构建神经网络的源程序。用以分析心电信号、脑电信号等等。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
load data
M1=20;
epo=15;
A=4;
B=18;
B2=B/2+1
N=500;
M=(A+1)*(B+1);
fora0=1:A+1;
forb0=1:B+1;
i=(B+1)*(a0-1)+b0;
b_init(i)=((b0-B2)/10)/(2^(-A)); a_init(i)=1/(2^(-A));
c_init(i)=(20-A)/2;
end
end
w0=ones(1,M);
for i=1:N
for j=1:M
t=x(i);
t= a_init(j)*t-b_init(j);
%P0(i,j)= (cos(1.75*t)*exp(-t*t/2))/2^c_init(j);
P0(i,j)= ((1-t*t)*exp(-t*t/2))/2^c_init(j);
end
end
%calculation of output of network
for i=1:N
u=0;
for j=1:M
u=u+w0(j)*P0(i,j);%w0?aè¨?μ
end
y0(i)= u;% y(p)= u=??W(j)*phi(p,j)= ??W(j)* |μj(t)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:M
W(:,k)=P0(:,k);
end
for k=2:M
u=0;
for i=1:k-1
aik(i)=(P0(:,k)'*W(:,i))/(W(:,i)'*W(:,i));
u=u+aik(i) *W(:,i);
end
W(:,k)=P0(:,k)-u;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:M
g(i)= (W(:,i)'*d')/( W(:,i)'* W(:,i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=0;
for i=1:M
u=u+g(i)*(W(:,i)'*W(:,i));
end
DD=u;
for i=1:M
Erro(i)=(g(i)^2)*(W(:,i)'*W(:,i))/DD;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=1;
while(k<=M1)
u=1;
fori=2:M
if abs(Erro(u))<abs(Erro(i));
u=i
else u=u
end
end
I(k)=u;
Erro(u)=0
k=k+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:M1
u=I(k);
a(k)=a_init(u);
b(k)=b_init(u);
c(k)=c_init(u);
w1(k)=w0(u);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epoch=1;
error=0.1;
err=0.0001;
lin=0.5;
while (error>=err & epoch<=epo)
for i=1:N
for j=1:M1
t=x(i);
t= a(j)*t-b(j);
%P1(i,j)= (cos(1.75*t)*exp(-t*t/2))/2^c(j);
P1(i,j)= ((1-t*t)*exp(-t*t/2))/2^c(j);
end
end
%calculation of output of network
for i=1:N
u=0;
for j=1:M1
u=u+w1(j)*P1(i,j);
end
y1(i)= u;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=0;
for i=1:N
u=u+(d(i)-y1(i))^2;
end
u=u/2;%u=1/2??(d-p)^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if u>error
lin=lin*0.8;
else
lin=lin*1.2;
end
error=u; %error=u=1/2??(d-p)^2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:M1
u=0;
for i=1:N
u=u+(d(i)-y1(i))*P1(i,j);
end
EW(j)=-u;
end
if epoch==1
SW=-EW;
w1_=w1;
else
SW=-EW+((EW*EW')*SW_)/(EW_*EW_');
end
EW_=EW;
SW_=SW;
w1=w1_+SW*lin;
w1_=w1;
%number of epoch increase by 1
epoch=epoch+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1);plot(x,d);
subplot(2,1,2);plot(x,y1);
利用小波和霍夫曼对单声道文件进行压缩编码 并解码输出
x=wavread('vqdeout.wav');subplot(3,2,1);
plot(x);
=wavedec(x,3,'db1');
cA3=appcoef(c,l,'db1',3);
= detcoef(c,l,);
subplot(3,2,2);
plot(cA3);
subplot(3,2,3);
plot(cD1);
subplot(3,2,4);
plot(cD2);
subplot(3,2,5);
plot(cD3);
%quantize
=py_quantize(cA3);
=quantiz(cA3,partition,codebook);
=py_quantize(cD1);
=quantiz(cD1,partition1,codebook1);
=py_quantize(cD2);
=quantiz(cD2,partition2,codebook2);
=py_quantize(cD3);
=quantiz(cD3,partition3,codebook3);
%probability
g=py_probability(quants,l(1,1),codebook);
g1=py_probability(quants1,l(4,1),codebook1);
g2=py_probability(quants2,l(3,1),codebook2);
g3=py_probability(quants3,l(2,1),codebook3);
%huffman编解码
= huffmandict(codebook,g);
comp = huffmanenco(quants,dict); % Encode the data.
dsig = huffmandeco(comp,dict);
isequal(quants,dsig)
%huffman编解码
= huffmandict(codebook1,g1);
comp1 = huffmanenco(quants1,dict1); % Encode the data.
dsig1 = huffmandeco(comp1,dict1);
isequal(quants1,dsig1)
%huffman编解码
= huffmandict(codebook2,g2);
comp2 = huffmanenco(quants2,dict2); % Encode the data.
dsig2 = huffmandeco(comp2,dict2);
isequal(quants2,dsig2)
%huffman编解码
= huffmandict(codebook3,g3);
comp3 = huffmanenco(quants3,dict3); % Encode the data.
dsig3 = huffmandeco(comp3,dict3);
isequal(quants3,dsig3)
%dsig3=dsig3';
%dsig2=dsig2';
%dsig1=dsig1';
%dsig=dsig';
%重构
cA2=idwt(dsig,dsig3,'db1');
icA2(1:13950)=cA2(1:13950);
cA1=idwt(icA2,dsig2,'db1');
icA1(1:27899)=cA1(1:27899);
XX=idwt(icA1,dsig1,'db1');
XX=XX';
subplot(3,2,6);
plot(XX);
MSE=sum((x-XX).^2)/length(XX);
PSNR=10*log10((max(max(XX)))^2/MSE);
用小波神经网络来对时间序列进行预测
%File name : nprogram.m
%Description : This file reads thedata from its source into their respectivematrices prior to
% performing wavelet decomposition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and variables
clc;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% user desired resolution level (Tested: resolution = 2 is best)
level = menu('Enter desired resolution level: ', '1',...
'2 (Select this for testing)', '3', '4');
switch level
case 1, resolution = 1;
case 2, resolution = 2;
case 3, resolution = 3;
case 4, resolution = 4;
end
msg = ['Resolution level to be used is ', num2str(resolution)];
disp(msg);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% user desired amount of data to use
data = menu('Choose amount of data to use: ', '1 day', '2 days', '3 days', '4 days',...
'5 days', '6 days', '1 week (Select this for testing)');
switch data
case 1, dataPoints = 48; %1 day = 48 points
case 2, dataPoints = 96; %2 days = 96 points
case 3, dataPoints = 144; %3 days = 144 points
case 4, dataPoints = 192; %4 days = 192 points
case 5, dataPoints = 240; %5 days = 240 points
case 6, dataPoints = 288; %6 days = 288 points
case 7, dataPoints = 336; %1 weeks = 336 points
end
msg = ['No. of data points to be used is ', num2str(dataPoints)];
disp(msg);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Menu for data set selection
select = menu('Use QLD data of: ', 'Jan02',...
'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02');
switch select
case 1, demandFile = 'DATA200601_QLD1';
case 2, demandFile = 'DATA200602_QLD1';
case 3, demandFile = 'DATA200603_QLD1';
case 4, demandFile = 'DATA200604_QLD1';
case 5, demandFile = 'DATA200605_QLD1';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reading the historical DEMAND data into tDemandArray
selectedDemandFile=;
...
= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
%Display no. of points in the selected time series demand data
= size(tDemandArray);
msg = ['The no. of points in the selected Demand data is ', num2str(demandDataPoints)];
disp(msg);
%Decompose historical demand data signal
= swtmat(tDemandArray, resolution, 'db2');
approx = dD(resolution, :);
%Plot the original demand data signal
figure (1);
subplot(resolution + 2, 1, 1); plot(tDemandArray(1: dataPoints))
legend('Demand Original');
title('QLD Demand Data Signal');
%Plot the approximation demand data signal
for i = 1 : resolution
subplot(resolution + 2, 1, i + 1); plot(approx(1: dataPoints))
legend('Demand Approximation');
end
%After displaying approximation signal, display detail x
for i = 1: resolution
if( i > 1 )
detail(i, :) = dD(i-1, :)- dD(i, :);
else
detail(i, :) = tDemandArray' - dD(1, :);
end
if i == 1
subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
legendName = ['Demand Detail ', num2str(i)];
legend(legendName);
else
subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))
legendName = ['Demand Detail ', num2str(i)];
legend(legendName);
end
i = i + 1;
end
%Normalising approximation demand data
maxDemand = max(approx'); %Find largest component
normDemand = approx ./ maxDemand; %Right divison
maxDemandDetail = [ ];
normDemandDetail = [, ];
detail = detail + 4000;
for i = 1: resolution
maxDemandDetail(i) = max(detail(i, :));
normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reading the historical historical PRICE data into rrpArray
selectedPriceFile = ;
...
= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');
%Display no. of points in Price data
= size(rrpArray);
msg = ['The no. of points in Price data data is ', num2str(noOfDataPoints)];
disp(msg);
%Decompose historical Price data signal
= swtmat(rrpArray, resolution, 'db2');
approxP = dP(resolution, :);
%Plot the original Price data signal
figure (2);
subplot(resolution + 3, 1, 1); plot(rrpArray(2: dataPoints))
legend('Price Original');
title('Price Data Signal');
%Plot the approximation Price data signal
for i = 1 : resolution
subplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints))
legend('Price Approximation');
end
%After displaying approximation signal, display detail x
for i = 1: resolution
if( i > 1 )
detailP(i, :) = dP(i-1, :)- dP(i, :);
else
detailP(i, :) = rrpArray' - dP(1, :);
end
if i == 1
=butter(1,0.65,'low');
result =filter(B,A, detailP(i, 1: dataPoints));
subplot(resolution + 3, 1, resolution - i + 4);plot(result(i, 2: dataPoints))
legendName = ['low pass filter', num2str(i)];
legend(legendName);
subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
legendName = ['Price Detail ', num2str(i)];
legend(legendName);
else
subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints))
legendName = ['Price Detail ', num2str(i)];
legend(legendName);
end
i = i + 1;
end
%Normalising approximation Price data
maxPrice = max(approxP'); %Find largest component
normPrice = approxP ./ maxPrice; %Right divison
maxPriceDetail = [ ];
normPriceDetail = [, ];
detailP = detailP + 40;
for i = 1: resolution
maxPriceDetail(i) = max(detailP(i, :));
normPriceDetail(i, :) = detailP(i, :) ./maxPriceDetail(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function here allows repetitive options to,
% 1) Create new NNs, 2) Retrain the existing NNs,
% 3) Perform load demand forecasting and 4) Quit session
while (1)
choice = menu('Please select one of the following options: ',...
'CREATE new Neural Networks',...
'Perform FORECASTING of load demand', 'QUIT session...');
switch choice
case 1, scheme = '1';
case 2, scheme = '2';
case 3, scheme = '3';
case 4, scheme = '4';
end
%If scheme is 'c', call <nCreate> to create new NNs, train them then perform forecast
if(scheme == '1')
nCreate;
end
%If scheme is 'r', call <nRetrain> to retrain the existing NNs
if(scheme == '2')
nRetrain;
end
%If scheme is 'f', call <nForecast> to load the existing NN model
if(scheme == '3')
nForecast;
end
%If scheme is 'e', verifies and quit session if 'yes' is selected else continue
if(scheme == '4')
button = questdlg('Quit session?', 'Exit Dialog','Yes','No','No');
switch button
case 'Yes', disp(' ');
disp('Session has ended!!');
disp(' ');
break;
case 'No', quit cancel;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : ncreate.m
%Description : This file prepares the input &output data for the NNs. It creates then trains the
% NNs to mature them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and set target demand ouput to start at point 2
clc;
targetStartAt = 2;
disp('Program will now CREATE a Neural Network for training and forecasting...');
disp(' ');
disp('To capture the pattern of the signal, the model is ')
disp('set to accept dataPoints x 2 sets of training examples; ');
disp('1 set of demand + 1 sets of price. ');
disp(' ');
disp('The normalised demand data <point 2>, is to be taken as the ')
disp('output value for the first iteration of training examples. ');
disp(' ');
disp('Press ENTER key to continue...');
pause;
%Clear command screen then prompt for no. of training cycles
%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
clc;
cycle = input('Input the number of training cycles: ');
numOfTimes = resolution + 1;
%Creating and training the NNs for the respective
%demand and price coefficient signals
for x = 1: numOfTimes
%Clearing variables
clear targetDemand;
clear inputs;
clear output;
clc;
if(x == 1)
neuralNetwork = ['Neural network settings for approximation level ',
num2str(resolution)];
else
neuralNetwork = ['Neural network settings for detail level ', num2str(x - 1)];
end
disp(neuralNetwork);
disp(' ');
%Set no. of input nodes and hidden neurons for the
%respective demand and price coefficient signal
numOfInputs = 2;
inputValue = ['Number of neural network INPUT units is set at ', num2str(numOfInputs)];
disp(inputValue);
disp(' ');
numOfOutput = 1;
outValue = ['Output is set to ', num2str(numOfOutput)];
disp(outValue);
disp(' ');
numOfHiddens = input('Enter the no. of HIDDEN units for the NN hidden : ');
hiddenValue = ['Number of neural network HIDDEN units is set at ',
num2str(numOfHiddens)];
disp(hiddenValue);
disp(' ');
%Setting no. of training examples
trainingLength = dataPoints;
%Set target outputs of the training examples
if(x == 1)
targetDemand = normDemand(targetStartAt: 1 + trainingLength);
else
targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + trainingLength);
end
%Preparing training examples
%Setting training i/ps to be 2 with user defined no. of iterations (dataPoints)
y = 0;
while y < trainingLength
if(x == 1)
inputs(1, y + 1) = normDemand(y + 1);
inputs(2, y + 1) = normPrice(y + 1);
else
inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
output(y + 1, :) = targetDemand(y + 1);
y = y + 1;
end
inputs = (inputs');
%Setting no. of training cycles
= size(targetDemand); % <== tells the NN how long is 1 cycle;
size(targetDemand)
clear nn;
%Create new neural network for respective coefficient signal
%NET = MLP(NIN, NHIDDEN, NOUT, FUNC)
nn = mlp(numOfInputs, numOfHiddens, numOfOutput, 'linear');
%NN options
options = zeros(1, 18);
options(1) = 1; %Provides display of error values
options(14) = cycle * ni * np;
%Training the neural network
%netopt(net, options, x, t, alg);
nn = netopt(nn, options, inputs, output, 'scg');
%Save the neural network
filename = ['nn', num2str(x)];
save(filename, 'nn');
disp(' ');
msg = ['Neural network successfully CREATED and saved as => ', filename];
disp(msg);
if(x < 3)
disp(' ');
disp('Press ENTER key to continue training the next NN...');
else
disp(' ');
disp('Model is now ready to forecast!!');
disp(' ');
disp('Press ENTER key to continue...');
end
pause;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : nforecast.m
%Description : This file loads the trained NNs for load forcasting and %recombines the predicted
% outputs from the NNs to form the final predicted load series.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Clear command screen and variables
clc;
clear forecastResult;
clear actualDemand;
clear predicted;
clear actualWithPredicted;
disp('Neural networks loaded, performing electricity demand forecast...');
disp(' ');
pointsAhead = input('Enter no. of points to predict (should be < 336): ');
numOfTimes = resolution + 1;
%Predict coefficient signals using respective NNs
for x = 1 : numOfTimes
%Loading NN
filename = ['nn', num2str(x)];
clear nn
load(filename);
clear in;
numOfInputs = nn.nin;
y = 0;
%Preparing details to forecast
while y < pointsAhead
if(x == 1)
propData(1, y + 1) = normDemand(y + 1);
propData(2, y + 1) = normPrice(y + 1);
else
propData(1, y + 1) = normDemandDetail(x - 1, y + 1);
propData(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
y = y + 1;
end
if(x == 1)
forecast = mlpfwd(nn, propData') * maxDemand;
else
forecast = mlpfwd(nn, propData') * maxDemandDetail(x - 1) - 4000;
end
forecastResult(x, :) = forecast';
end
%For debugging
% forecastResult
actualDemand = tDemandArray(2: 1 + pointsAhead);
predicted = sum(forecastResult, 1)';
%Calculating Mean Absolute Error
AbsError = abs(actualDemand - predicted(1: pointsAhead)) ./ actualDemand;
msg = ['Mean Absolute Error = ', num2str(mean(AbsError(1: pointsAhead))), ' !!'];
disp(' ');
disp(msg);
%Plot actual time series against predicted result
figure(3)
actualWithPredicted(:, 1) = actualDemand;
actualWithPredicted(:, 2) = predicted(1: pointsAhead);
plot(actualWithPredicted);
graph = ['Mean Absolute Error = ', num2str(mean(AbsError))];
title(graph);
legend('Actual', 'Forecasted');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File name : nretrain.m
%Description : This file loads the existing NNs and trains them again.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
%Prompt for the starting point for training
disp('Program will now RETRAIN the Neural Networks ')
disp('with the SAME intial data series again...');
disp(' ');
disp('To capture the pattern of the signal, the model is ')
disp('set to accept dataPoints x 2 sets of training examples; ');
disp('1 set of demand + 1 sets of price. ');
disp(' ');
disp('The normalised demand data <point 2>, is to be taken as the ')
disp('output value for the first iteration of training examples. ');
disp(' ');
msg = ['Data points to be used for reTraining the NNs is from 1 to ',...
num2str(dataPoints)];
disp(msg);
disp(' ');
disp('Press ENTER key to continue...');
pause;
%Clear command screen
clc;
%Prompt for no. of training cycles
%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)
cycle = input('Input number of cycles to retrain the NNs: ');
numOfTimes = resolution + 1;
%Loading existing NNs for training
for x = 1: numOfTimes
%Re-initialising variables
clear targetDemand;
clear inputs;
clear output;
clc;
%Loading NN for the respective demand and temperature coefficient signals
filename = ['nn', num2str(x)];
clear nn
load(filename);
%Getting the size of NN
numOfInputs = nn.nin;
numOfHiddens = nn.nhidden;
numOfOutput = 1;
%Setting length of reTraining examples and target outputs
reTrainLength = dataPoints;
targetLength = reTrainLength;
targetStartAt = 2;
%Set target outputs of the training examples
if(x == 1)
targetDemand = normDemand(targetStartAt: 1 + targetLength);
else
targetDemand = normDemandDetail(x - 1, targetStartAt: 1 + targetLength);
end
%Preparing training examples
%Setting training i/ps to be 2 with user set no. of iterations (dataPoints)
y = 0;
while y < reTrainLength
if(x == 1)
inputs(1, y + 1) = normDemand(y + 1);
inputs(2, y + 1) = normPrice(y + 1);
else
inputs(1, y + 1) = normDemandDetail(x - 1, y + 1);
inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);
end
output(y + 1, :) = targetDemand(y + 1);
y = y + 1;
end
inputs = (inputs');
%Setting no. of training cycles
= size(targetDemand); % <== tells the NN how long is 1 cycle;
size(targetDemand) %With reference to line 106
%NN options
options = zeros(1, 18);
options(1) = 1; %Provides display of error values
options(14) = cycle * ni * np;
%Training the neural network
%netopt(net, options, x, t, alg);
nn = netopt(nn, options, inputs, output, 'scg');
%Save the neural network
filename = ['nn', num2str(x)];
save(filename, 'nn');
disp(' ');
msg = ['Neural network => ', filename, ' <= successfully RETRAINED and saved!! '];
disp(msg);
if(x < 3)
disp(' ');
disp('Press ENTER key to continue training the next NN...');
else
disp(' ');
disp('Model is now ready to forecast again!!');
disp(' ');
disp('Press ENTER key to continue...');
end
pause;
end
基于小波特征提取方法的图象匹配算法
function feature_wl=getf_wavelet(Input)A=double(Input);
Xrgb=0.2990*A(:,:,1)+0.5870*A(:,:,2)+0.1140*A(:,:,3);
NbColor=255;
%WCODEMAT Extended pseudocolor matrix scaling. Y = WCODEMAT(X,NBCODES,OPT,ABSOL) returns a coded version of input
% matrix X if ABSOL=0, or ABS(X) if ABSOL isnonzero, using the first NBCODES integers.
%Coding can be done row-wise (OPT='row' or 'r'), columnwise (OPT='col' or 'c'), or globally (OPT='mat' or 'm').
%Coding uses a regular grid between the minimum and the maximum values of each row (column or matrix,respectively).
%Y = WCODEMAT(X,NBCODES) is equivalent toY = WCODEMAT(X,NBCODES,'mat',1).
X1=wcodemat(Xrgb,NbColor);%伪彩矩阵压缩
% Function 'wavedec2' Multilevel 2-D wavelet decomposition. = WAVEDEC2(X,N,'wname') returns the wavelet
% decomposition of the matrix X at level N, using the wavelet named in string 'wname' (see WFILTERS).
% Outputs are the decomposition vector C and the corresponding bookkeeping matrix S.
% N must be a strictly positive integer (see WMAXLEV).
=wavedec2(X1,4,'sym4');%二维小波分解
% Function 'detcoef2' Extract 2-D detail coefficients.
% D = detcoef2(O,C,S,N) extracts from the wavelet decomposition structure , the horizontal, vertical
% or diagonal detail coefficients for O = 'h' (or 'v' or 'd',respectively), at level N. N must
% be an integer such that 1 <= N <= size(S,1)-2.
ch11=detcoef2('h',c2,l2,1);%提取小波分解细节系数
ch12=detcoef2('h',c2,l2,2);
ch13=detcoef2('h',c2,l2,3);
ch14=detcoef2('h',c2,l2,4);
cv11=detcoef2('v',c2,l2,1);
cv12=detcoef2('v',c2,l2,2);
cv13=detcoef2('v',c2,l2,3);
cv14=detcoef2('v',c2,l2,4);
cd11=detcoef2('d',c2,l2,1);
cd12=detcoef2('d',c2,l2,2);
cd13=detcoef2('d',c2,l2,3);
cd14=detcoef2('d',c2,l2,4);
% Function 'appcoef2' Extract 2-D approximation coefficients.APPCOEF2 computes the approximation coefficients of a
% two-dimensional signal.A = APPCOEF2(C,S,'wname',N) computes the approximation coefficients at level N using the wavelet decomposition
% structure (see WAVEDEC2). 'wname' is a string containing the wavelet name.Level N must be an integer such that 0 <= N <= size(S,1)-2.
ca14=appcoef2(c2,l2,'sym4',4);%提取小波分解概貌系数
%==========compute mean & var(求分解后的每个图象的均值和方差)==================
%image x1's u & std
Csize1=size(ch11);
s1=0;
for i=1:Csize1(1)
for j=1:Csize1(2)
s1=s1+abs(ch11(i,j));%abs求绝对值
end
end
u1=(1/(Csize1(1)*Csize1(2)))*s1;
st1=0;
for i=1:Csize1(1)
for j=1:Csize1(2)
st1=st1+(abs(ch11(i,j))-u1)*(abs(ch11(i,j))-u1);
end
end
z1=sqrt((1/(Csize1(1)*Csize1(2)-1))*st1);
%======================================================
Csize2=size(ch12);
s2=0;
for i=1:Csize2(1)
for j=1:Csize2(2)
s2=s2+abs(ch12(i,j));
end
end
u2=(1/(Csize2(1)*Csize2(2)))*s2;
st2=0;
for i=1:Csize2(1)
for j=1:Csize2(2)
st2=st2+(abs(ch12(i,j))-u2)*(abs(ch12(i,j))-u2);
end
end
z2=sqrt((1/(Csize2(1)*Csize2(2)-1))*st2);
%===========================================================
Csize3=size(ch13);
s3=0;
for i=1:Csize3(1)
for j=1:Csize3(2)
s3=s3+abs(ch13(i,j));
end
end
u3=(1/(Csize3(1)*Csize3(2)))*s3;
st3=0;
for i=1:Csize3(1)
for j=1:Csize3(2)
st3=st3+(abs(ch13(i,j))-u3)*(abs(ch13(i,j))-u3);
end
end
z3=sqrt((1/(Csize3(1)*Csize3(2)-1))*st3);
%================================================================
Csize3=size(ch14);
s3=0;
for i=1:Csize3(1)
for j=1:Csize3(2)
s3=s3+abs(ch14(i,j));
end
end
u4=(1/(Csize3(1)*Csize3(2)))*s3;
st3=0;
for i=1:Csize3(1)
for j=1:Csize3(2)
st3=st3+(abs(ch14(i,j))-u4)*(abs(ch14(i,j))-u4);
end
end
z4=sqrt((1/(Csize3(1)*Csize3(2)-1))*st3);
%=================================================================
Csize4=size(cv11);
s4=0;
for i=1:Csize4(1)
for j=1:Csize4(2)
s4=s4+abs(cv11(i,j));
end
end
u5=(1/(Csize4(1)*Csize4(2)))*s4;
st4=0;
for i=1:Csize4(1)
for j=1:Csize4(2)
st4=st4+(abs(cv11(i,j))-u5)*(abs(cv11(i,j))-u5);
end
end
z5=sqrt((1/(Csize4(1)*Csize4(2)-1))*st4);
%===============================================================
Csize5=size(cv12);
s5=0;
for i=1:Csize5(1)
for j=1:Csize5(2)
s5=s5+abs(cv12(i,j));
end
end
u6=(1/(Csize5(1)*Csize5(2)))*s5;
st5=0;
for i=1:Csize5(1)
for j=1:Csize5(2)
st5=st5+(abs(cv12(i,j))-u6)*(abs(cv12(i,j))-u6);
end
end
z6=sqrt((1/(Csize5(1)*Csize5(2)-1))*st5);
%================================================================
Csize6=size(cv13);
s6=0;
for i=1:Csize6(1)
for j=1:Csize6(2)
s6=s6+abs(cv13(i,j));
end
end
u7=(1/(Csize6(1)*Csize6(2)))*s6;
st6=0;
for i=1:Csize6(1)
for j=1:Csize6(2)
st6=st6+(abs(cv13(i,j))-u7)*(abs(cv13(i,j))-u7);
end
end
z7=sqrt((1/(Csize6(1)*Csize6(2)-1))*st6);
%==================================================================
Csize6=size(cv14);
s6=0;
for i=1:Csize6(1)
for j=1:Csize6(2)
s6=s6+abs(cv14(i,j));
end
end
u8=(1/(Csize6(1)*Csize6(2)))*s6;
st6=0;
for i=1:Csize6(1)
for j=1:Csize6(2)
st6=st6+(abs(cv14(i,j))-u8)*(abs(cv14(i,j))-u8);
end
end
z8=sqrt((1/(Csize6(1)*Csize6(2)-1))*st6);
%==============================================================
Csize7=size(cd11);
s7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
s7=s7+abs(cd11(i,j));
end
end
u9=(1/(Csize7(1)*Csize7(2)))*s7;
st7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
st7=st7+(abs(cd11(i,j))-u9)*(abs(cd11(i,j))-u9);
end
end
z9=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
%====================================================================
Csize7=size(cd12);
s7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
s7=s7+abs(cd12(i,j));
end
end
u10=(1/(Csize7(1)*Csize7(2)))*s7;
st7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
st7=st7+(abs(cd12(i,j))-u10)*(abs(cd12(i,j))-u10);
end
end
z10=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
%=================================================
Csize7=size(cd13);
s7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
s7=s7+abs(cd13(i,j));
end
end
u11=(1/(Csize7(1)*Csize7(2)))*s7;
st7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
st7=st7+(abs(cd13(i,j))-u11)*(abs(cd13(i,j))-u11);
end
end
z11=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
%=================================================
Csize7=size(cd14);
s7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
s7=s7+abs(cd14(i,j));
end
end
u12=(1/(Csize7(1)*Csize7(2)))*s7;
st7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
st7=st7+(abs(cd14(i,j))-u12)*(abs(cd14(i,j))-u12);
end
end
z12=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
%=================================================
Csize7=size(ca14);
s7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
s7=s7+abs(ca14(i,j));
end
end
u13=(1/(Csize7(1)*Csize7(2)))*s7;
st7=0;
for i=1:Csize7(1)
for j=1:Csize7(2)
st7=st7+(abs(ca14(i,j))-u13)*(abs(ca14(i,j))-u13);
end
end
z13=sqrt((1/(Csize7(1)*Csize7(2)-1))*st7);
%semble feature vector========================================
% 应用准则函数把对分类影响不大的特征去除
feature_wl=;
%feature_wl=; 请问有没有关于时频分析
的chirplet变换matlab程序呢？谢谢！ 请问楼主:
我刚刚接触matlab和小波变换, =wavedec2(X1,4,'sym4');中对其参数有什么要求吗?为什么每次执行到这条语句的时候就出问题啊? 楼主真慷慨，谢谢啦！请问有没有用小波神经网络进行图象识别的程序？ 有没有关于说话人识别中利用小波变换提取特征参数的程序，不胜感激，邮箱zbl@mail2.lut.cn LZ有关于小波奇异性检测方面的源程序么？
谢谢了:loveliness:
真是个慷慨的好LZ啊！