% mallet_wavelet.m

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

% 此函数用于消噪处理

%角度赋值

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

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

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

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

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

low_construct=l;

L_fre=20; %滤波器长度

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

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

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

coefficient=-1;

else

coefficient=1;

end

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

end

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

L_signal=100; %信号长度

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

f=10;

t=0.001;

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

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

zero2=zeros(1,30);

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

y_noise=y+noise;

figure(1);

subplot(2,1,1);

plot(y);

title('原信号');

subplot(2,1,2);

plot(y_noise);

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

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

check2=sum(low_decompose);

check3=norm(high_decompose);

check4=norm(low_decompose);

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

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

h_fre=conv(y_noise,high_decompose);

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

figure(2);

subplot(2,1,1)

plot(l_fre_down);

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

subplot(2,1,2);

plot(h_fre_down);

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

% 消噪处理

for i_decrease=31:44;

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

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

end

end

l_fre_pull=dyadup(l_fre_down); %0差值

h_fre_pull=dyadup(h_fre_down);

l_fre_denoise=conv(low_construct,l_fre_pull);

h_fre_denoise=conv(high_construct,h_fre_pull);

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

h_fre_keep=wkeep(h_fre_denoise,L_signal);

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

%平滑处理

for j=1:2

for i=60:70;

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

end;

end;

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

figure(3);

subplot(3,1,1)

plot(y);

ylabel('y'); %原信号

subplot(3,1,2);

plot(sig_denoise);

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

subplot(3,1,3);

plot(compare);

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

 

 

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

复制内容到剪贴板

代码:

clc;clear;
%% 1.正弦波定义
f1=50; % 频率1
f2=100; % 频率2
fs=2*(f1+f2); % 采样频率
Ts=1/fs; % 采样间隔
N=120; % 采样点数
n=1:N;
y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合
figure(1)
plot(y);
title('两个正弦信号')
figure(2)
stem(abs(fft(y)));
title('两信号频谱')
%% 2.小波滤波器谱分析
h=wfilters('db30','l'); % 低通
g=wfilters('db30','h'); % 高通
h=[h,zeros(1,N-length(h))]; % 补零（圆周卷积，且增大分辨率变于观察）
g=[g,zeros(1,N-length(g))]; % 补零（圆周卷积，且增大分辨率变于观察）
figure(3);
stem(abs(fft(h)));
title('低通滤波器图')
figure(4);
stem(abs(fft(g)));
title('高通滤波器图')
%% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)
sig1=ifft(fft(y).*fft(h)); % 低通(低频分量)
sig2=ifft(fft(y).*fft(g)); % 高通(高频分量)
figure(5); % 信号图
subplot(2,1,1)
plot(real(sig1));
title('分解信号1')
subplot(2,1,2)
plot(real(sig2));
title('分解信号2')
figure(6); % 频谱图
subplot(2,1,1)
stem(abs(fft(sig1)));
title('分解信号1频谱')
subplot(2,1,2)
stem(abs(fft(sig2)));
title('分解信号2频谱')
%% 4.MALLET重构算法
sig1=dyaddown(sig1); % 2抽取
sig2=dyaddown(sig2); % 2抽取
sig1=dyadup(sig1); % 2插值
sig2=dyadup(sig2); % 2插值
sig1=sig1(1,[1:N]); % 去掉最后一个零
sig2=sig2(1,[1:N]); % 去掉最后一个零
hr=h(end:-1:1); % 重构低通
gr=g(end:-1:1); % 重构高通
hr=circshift(hr',1)'; % 位置调整圆周右移一位
gr=circshift(gr',1)'; % 位置调整圆周右移一位
sig1=ifft(fft(hr).*fft(sig1)); % 低频
sig2=ifft(fft(gr).*fft(sig2)); % 高频
sig=sig1+sig2; % 源信号
%% 5.比较
figure(7);
subplot(2,1,1)
plot(real(sig1));
title('重构低频信号');
subplot(2,1,2)
plot(real(sig2));
title('重构高频信号');
figure(8);
subplot(2,1,1)
stem(abs(fft(sig1)));
title('重构低频信号频谱');
subplot(2,1,2)
stem(abs(fft(sig2)));
title('重构高频信号频谱');
figure(9)
plot(real(sig),'r','linewidth',2);
hold on;
plot(y);
legend('重构信号','原始信号')
title('重构信号与原始信号比较')

 

2维小波变换经典程序

复制内容到剪贴板

代码:

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.分解结果显示
figure(1);
colormap(map);
subplot(2,1,1);
image(f); % 原始图像
title('original pic');
subplot(2,1,2);
image(abs(decompose_pic)); % 分解后图像
title('decomposed pic');
figure(2);
colormap(map);
subplot(2,2,1);
image(abs(lt_pic)); % 左上方为低频分量--fi(x)*fi(y)
title('\Phi(x)*\Phi(y)');
subplot(2,2,2);
image(abs(rt_pic)); % 矩阵右上为--fi(x)*psi(y)
title('\Phi(x)*\Psi(y)');
subplot(2,2,3);
image(abs(lb_pic)); % 矩阵左下为--psi(x)*fi(y)
title('\Psi(x)*\Phi(y)');
subplot(2,2,4);
image(abs(rb_pic)); % 右下方为高频分量--psi(x)*psi(y)
title('\Psi(x)*\Psi(y)');


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

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

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

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

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

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

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


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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 结果显示
figure(3);
colormap(map);
subplot(2,1,1);
image(f); % 源图像显示
title('original pic');

subplot(2,1,2);
image(abs(construct_pic)); % 重构源图像显示
title('reconstructed pic');
error=abs(construct_pic-f); % 重构图形与原始图像误值
figure(4);
mesh(error); % 误差三维图像
title('absolute error display');

 

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

复制内容到剪贴板

代码:

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


end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cubic_average（立方b样条）
% 均值插值 % 初始化
s=[0 0 1 0 0] % 正弦波
% f=50;
% ts=1/200;
% n=0:16;
% s=sin(2*pi*f*n*ts); % 系数
se=[1/8,6/8,1/8];
so=[4/8,4/8] % 循环
for p=1:10;
t=length(s)-1;
o(1:t)=s(1:t)*so(1)+s(2:t+1)*so(2);
e(1)=s(t+1)*se(1)+s(1)*se(2)+s(2)*se(3);
e(2:t)=s(1:t-1)*se(1)+s(2:t)*se(2)+s(3:t+1)*se(3);
e(t+1)=s(t)*se(1)+s(t+1)*se(2)+s(1)*se(3);
s([1:2:2*t+1])=e([1:t+1]);
s([2:2:2*t])=o([1:t]);
plot(s);
drawnow;
end; % 抽取
t=length(s); % 总长度
p=128; % 需要点数 % 间隔
d=(t-1)/p; % 最终尺度函数
r=s(2:d:t-1); % 画图
figure(2);
plot(r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cubic_subdivision（立方插值）
% 细分插值 % %% 初始化（尺度函数）
% s=[0,0,1,0,0]; % 正弦函数
n=1:20;
f=50;
ts=1/200;
s=sin(2*pi*f*n*ts); % 指数函数
% n=0:16;
% s=exp(n); % % 系数
a=[-1/16,9/16,9/16,-1/16]; % 循环
for p=1:4;
t=length(s)-1;
o(1)=s(4)*a(1)+s(1)*a(2)+s(2)*a(3)+s(3)*a(4);
o(2:t-1)=s(1:t-2)*a(1)+s(2:t-1)*a(2)+s(3:t)*a(3)+s(4:t+1)*a(4);
o(t)=s(t-2)*a(4)+s(t+1)*a(3)+s(t)*a(2)+s(t-1)*a(1);
s([1:2:2*t+1])=s([1:t+1]);
s([2:2:2*t])=o([1:t]);
plot(s);
drawnow;
end; % % 抽取
% t=length(s); % 总长度
% p=128; % 需要点数
%
% % 间隔
% d=(t-1)/p;
%
% % 最终尺度函数
% r=s(2:d:t-1);
%
% % 画图
% figure(2);
% plot(r);

 

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

复制内容到剪贴板

代码:

clear;clc;
% 1.生成信号
f=50; % 频率
fs=800; % 采样率
T=128; % 信号长度
n=1:T;
y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 信号
% y=circshift(y.',3).'; %%
2.正变换 l1=wfilters('db4','l')*sqrt(2)/2;
% 参考低通滤波器
l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通滤波器1
h1=wfilters('db4','h')*sqrt(2)/2; % 参考高通滤波器
h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通滤波器1 low1=ifft(fft(y).*fft(l1_zeros)); % 低频分量1
high1=ifft(fft(y).*fft(h1_zeros)); % 高频分量1 l2=dyadup(l1); % 原滤波器插值
l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通滤波器2
h2=dyadup(h1); % 原滤波器插值
h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通滤波器2 low2=ifft(fft(low1).*fft(l2_zeros)); % 低频分量2
high2=ifft(fft(low1).*fft(h2_zeros)); % 高频分量2
%% 3.反变换 lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重构低通滤波器2
hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重构高通滤波器2
lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重构低通滤波器1
hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重构高通滤波器1 lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重构低频分量1(lowr=low1)
r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重构源信号(r_s=y)
%% 4.绘图 figure(1);
plot(y);
title('源信号'); figure(2);
plot(low1,'r');
hold on;
plot(low2,'b');
legend('第一层低频','第二层低频'); figure(3);
plot(high1,'r');
hold on;
plot(high2,'b');
legend('第一层高频','第二层高频'); figure(4);
plot(low1,'r');
hold on;
plot(lowr,'b.');
legend('第一层低频','重构第一层低频'); figure(5);
plot(y,'r');
hold on;
plot(r_s,'b.');
legend('源信号','重构信号'); disp(norm(low1-lowr))
disp(norm(y-r_s))

 

此程序实现构造小波基

复制内容到剪贴板

代码:

% 此程序实现构造小波基
% periodic_wavelet.m function ss=periodic_wavelet; clear;clc; % global MOMENT; % 消失矩阶数
% global LEFT_SCALET; % 尺度函数左支撑区间
% global RIGHT_SCALET; % 尺度函数右支撑区间
% global LEFT_BASIS; % 小波基函数左支撑区间
% global RIGHT_BASIS; % 小波基函数右支撑区间
% global MIN_STEP; % 最小离散步长
% global LEVEL; % 计算需要的层数(离散精度）
% global MAX_LEVEL; % 周期小波最大计算层数
[s2,h]=scale_integer;
[test,h]=scalet_stretch(s2,h);
wave_base=wavelet(test,h);
ss=periodic_waveletbasis(wave_base);
function [s2,h]=scale_integer; % 本函数实现求解小波尺度函数离散整数点的值
% sacle_integer.m MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
h=wfilters('db10','r'); % 滤波器系数 h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1; for i=LEFT_SCALET+1:RIGHT_SCALET-1
for j=LEFT_SCALET+1:RIGHT_SCALET-1
k=2*i-j+1;
if (k>=1&k<=RIGHT_SCALET+1)
a(i,j)=h(k); % 矩阵系数矩阵
else
a(i,j)=0;
end
end
end [s,w]=eig(a); % 求特征向量，解的基
s1=s(:,1);
s2=[0;s1/sum(s1);0]; % 根据条件SUM(FI(T))=1,求解; % 本函数实现尺度函数经伸缩后的离散值
% scalet_stretch.m function [s2,h]=scalet_stretch(s2,h); MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
for j=1:LEVEL % 需要计算到尺度函数的层数
t=0;
for i=1:2:2*length(s2)-3 % 需要计算的离散点取值（0，1，2，3 -> 1/2, 3/2, 5/2)
t=t+1;
fi(t)=0;
for n=LEFT_SCALET:RIGHT_SCALET; % 低通滤波器冲击响应紧支撑判断
if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) % 小波尺度函数紧支撑判断
fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1); % 反复应用双尺度方程求解
end
end
end
clear s
n1=length(s2);
n2=length(fi);
for i=1:length(s2)+length(fi) % 变换后的矩阵长度
if (mod(i,2)==1)
s(i)=s2((i+1)/2); % 矩阵奇数下标为小波上一层(0,1,2,3)离散值
else
s(i)=fi(i/2); % 矩阵偶数下标为小波下一层(1/2,3/2,5/2)(经过伸缩变换后)的离散值
end
end
s2=s;
end
% 采用双尺度方程求解小波基函数 PSI(T)
% wavelet.m function wave_base=wavelet(test,h); MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
i=0;
for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS; % 小波基支撑长度
s=0;
for n=1-RIGHT_SCALET:1-LEFT_SCALET % g(n)取值范围
if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET) % 尺度函数判断
s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1); % 计算任意精度的小波基函数值
end
end
i=i+1;
wave_base(i)=s;
end

 

消失矩作用的程序

复制内容到剪贴板

代码:

clear;clc;
f=50;
T=0.001;
n=1:50;
y=sin(2*pi*f*n*T);
noise=[zeros(1,17),0.2*randn(1,15),zeros(1,18)];
y_noise=y+noise;
figure(1);
subplot(2,1,1);
plot(y);
title('signal');
subplot(2,1,2);
plot(y_noise);
title('noise & signal');
[yl2,yh2]=dwt(y,'db2');
[yl10,yh10]=dwt(y,'db10');
figure(2);
subplot(2,1,1);
plot(yl2);
title('db2 low frequency signal');
subplot(2,1,2);
plot(yh2);
title('db2 high frequency signal');
figure(3);
subplot(2,1,1);
plot(yl10);
title('db10 low frequency signal');
subplot(2,1,2);
plot(yh10);
title('db10 high frequency signal');
[yl2,yh2]=dwt(y_noise,'db2');
[yl10,yh10]=dwt(y_noise,'db10');
figure(4);
subplot(2,1,1);
plot(yl2);
title('db2 low frequency noise & signal');
subplot(2,1,2);
plot(yh2);
title('db2 high frequency noise & signal');
figure(5);
subplot(2,1,1);
plot(yl10);
title('db10 low frequency noise & signal');
subplot(2,1,2);
plot(yh10);
title('db10 high frequency noise & signal');

 

平移变换平移法(cycle_spinning)消除gibbs效应

复制内容到剪贴板

代码:

clear;
clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1.原始信号 f=50; % 信号频率
fs=800; % 采样频率
N=128; % 采样点 % 信号赋值
n=1:N;
y=sin(2*pi*f*n/fs); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.噪声
noise=0.4*rand(1,128); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3.染噪信号
y_noise=y+noise; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4.硬消噪采用cycle_spinning技术 % 累加量
z5=zeros(1,N); % 平移变换频移法
for i=1:N;
z=circshift(y_noise.',i-1).'; % 源信号右平移
[z1,z2]=lwt(z,'db3'); % 小波正变换
z2=zeros(1,N/2); % 高频分量全部为零（主要噪声，硬消噪）
z3=ilwt(z1,z2,'db3'); % 小波反变换
z4=circshift(z3.',-(i-1)).'; % 变换后信号左平移
z5=z5+z4/N; % 平均
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 5.显示 error=norm(y-z5)/norm(y); % 相对误差 figure(1);
subplot(2,1,1)
plot(y,'r');
legend('源信号');
subplot(2,1,2);
plot(y_noise);
legend('染噪信号');
figure(2);
subplot(2,1,1)
plot(y,'r');
legend('源信号');
title(error);
subplot(2,1,2);
plot(z5);
legend('消噪后信号');

 

基于小波消噪的雷达回波检测

基于小波消噪的雷达回波检测，可以检测雷达回波的有无及其准确的位置

复制内容到剪贴板 代码:close all;clear;clc;
signalw=1200;%信号宽度
m=120000; %数据采样点数
pfa=0.001;%虚警概率
dectectw=signalw/4;%滑窗宽度
a=wgn(1,m,13,50,'dbm','real');%噪声
b=gate(a,m,dectectw,pfa)%门限
signal=real(fmlin(signalw,0,0.25));%信号

%加高斯白噪声信号
mix=zeros(1,m);
signalpos=40000;
for l=signalpos:signalpos+signalw-1
    mix(1,l)=signal(l-(signalpos-1),1);
end
mix=mix+a;%加高斯白噪声信号

%检测信号有无及位置
g=0;
for l=1:dectectw
    g=g+mix(l)*(mix(l))';
   
end
   g0=g;
   for k=(dectectw+1):m
    g=g+mix(k)*(mix(k))'-mix(k-dectectw)*(mix(k-dectectw))';
    if (g>b)
        over=1
        index=k
        break
    end
end%检测信号有无及位置

pnoise=mean(abs(a).^2);
psignal=mean(abs(signal).^2);
snr=10*log10(psignal/pnoise) %信噪比

[mix1,c,l]=wden(mix,'rigrsure','s','sln',3,'db2');
[a1,c,l]=wden(a,'rigrsure','s','sln',3,'db2');
% figure(2)
% subplot(2,1,1)
% plot(a)
% subplot(2,1,2)
% plot(a1)
[signal1,c,l]=wden(signal,'rigrsure','s','sln',3,'db2');
% figure(3)
% subplot(2,1,1)
% plot(signal)
% subplot(2,1,2)
% plot(signal1)

%检测信号有无及位置
dectectw1=signalw*1/4;
b1=gate(a1,m,dectectw1,pfa)
g=0;
for l=1:dectectw1
    g=g+mix1(l)*(mix1(l))';
   
end
   g0=g;
   for k=(dectectw1+1):m
    g=g+mix1(k)*(mix1(k))'-mix1(k-dectectw1)*(mix1(k-dectectw1))';
    if (g>b1)
        over=1
        index=k
        break
    end
end%检测信号有无及位置

pnoise1=(sum(abs(mix1).^2)-sum(abs(signal1).^2))/m;
psignal1=mean(abs(signal1).^2);
snr1=10*log10(psignal1/pnoise1) %信噪比

subplot(2,1,1)
plot(mix)
axis([0 60000 -4 4])
TITLE('消噪前信号')
str={['SNR=',num2str(snr),'dB']};
text(45000,3,str)
subplot(2,1,2)
plot(mix1)
axis([0 60000 -4 4])
TITLE('消噪后信号')
str={['SNR=',num2str(snr1),'dB']};
text(45000,2.5,str)

二维小波变换（正和逆变换）

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

复制内容到剪贴板 代码:function [varargout]=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=[1 1]/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=[0:7];
        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=[0:7];
        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=[0:17];
        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=[0:9];
        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 [c,s]=wavefast(x,n,varargin)
error(nargchk(3,4,nargin));
if nargin==3
    if ischar(varargin{1})
        [lp,hp]=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
    [app,keep]=symextend(app,fl);
    rows=symconv(app,hp,'row',fl,keep);
    coefs=symconv(rows,hp,'col',fl,keep);
    c=[coefs(:)' c]; s=[size(coefs);s];
    coefs=symconv(rows,lp,'col',fl,keep);
    c=[coefs(:)' c];
    rows=symconv(app,lp,'row',fl,keep);
    coefs=symconv(rows,hp,'col',fl,keep);
    c=[coefs(:)' c];
    app=symconv(rows,lp,'col',fl,keep);
end
c=[app(:)' c]; s=[size(app);s];
function [y, keep]=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 [varargout]=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(['[C S] 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)
            [lp,hp]=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)
            [lp,hp]=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=[a(:)' nc];
    ns=ns(3:end,:);
    ns=[ns(1,:);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([2 1].*size(x));
y(1:2:end,:)=x;
y=conv2(y,f1');
z=zeros([1 2].*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);

二维小波变换（正和逆变换）

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

复制内容到剪贴板 代码:function [varargout]=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=[1 1]/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=[0:7];
        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=[0:7];
        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=[0:17];
        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=[0:9];
        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 [c,s]=wavefast(x,n,varargin)
error(nargchk(3,4,nargin));
if nargin==3
    if ischar(varargin{1})
        [lp,hp]=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
    [app,keep]=symextend(app,fl);
    rows=symconv(app,hp,'row',fl,keep);
    coefs=symconv(rows,hp,'col',fl,keep);
    c=[coefs(:)' c]; s=[size(coefs);s];
    coefs=symconv(rows,lp,'col',fl,keep);
    c=[coefs(:)' c];
    rows=symconv(app,lp,'row',fl,keep);
    coefs=symconv(rows,hp,'col',fl,keep);
    c=[coefs(:)' c];
    app=symconv(rows,lp,'col',fl,keep);
end
c=[app(:)' c]; s=[size(app);s];
function [y, keep]=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 [varargout]=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(['[C S] 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)
            [lp,hp]=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)
            [lp,hp]=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=[a(:)' nc];
    ns=ns(3:end,:);
    ns=[ns(1,:);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([2 1].*size(x));
y(1:2:end,:)=x;
y=conv2(y,f1');
z=zeros([1 2].*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);

使用小波包变换分析信号的MATLAB程序 复制内容到剪贴板 代码:%t=0.001:0.001:1; t=1:1000; s1=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t)); for t=1:500; s2(t)=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t)); end for t=501:1000; s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t)); end subplot(9,2,1) plot(s1) title('原始信号') ylabel('S1') subplot(9,2,2) plot(s2) title('故障信号') ylabel('S2') wpt=wpdec(s1,3,'db1','shannon'); %plot(wpt); s130=wprcoef(wpt,[3,0]); s131=wprcoef(wpt,[3,1]); s132=wprcoef(wpt,[3,2]); s133=wprcoef(wpt,[3,3]); s134=wprcoef(wpt,[3,4]); s135=wprcoef(wpt,[3,5]); s136=wprcoef(wpt,[3,6]); s137=wprcoef(wpt,[3,7]); s10=norm(s130); s11=norm(s131); s12=norm(s132); s13=norm(s133); s14=norm(s134); s15=norm(s135); s16=norm(s136); s17=norm(s137); st10=std(s130); st11=std(s131); st12=std(s132); st13=std(s133); st14=std(s134); st15=std(s135); st16=std(s136); st17=std(s137); disp('正常信号的特征向量'); snorm1=[s10,s11,s12,s13,s14,s15,s16,s17] std1=[st10,st11,st12,st13,st14,st15,st16,st17] subplot(9,2,3);plot(s130); ylabel('S130'); subplot(9,2,5);plot(s131); ylabel('S131'); subplot(9,2,7);plot(s132); ylabel('S132'); subplot(9,2,9);plot(s133); ylabel('S133'); subplot(9,2,11);plot(s134); ylabel('S134'); subplot(9,2,13);plot(s135); ylabel('S135'); subplot(9,2,15);plot(s136); ylabel('S136'); subplot(9,2,17);plot(s137); ylabel('S137'); wpt=wpdec(s2,3,'db1','shannon'); %plot(wpt); s230=wprcoef(wpt,[3,0]); s231=wprcoef(wpt,[3,1]); s232=wprcoef(wpt,[3,2]); s233=wprcoef(wpt,[3,3]); s234=wprcoef(wpt,[3,4]); s235=wprcoef(wpt,[3,5]); s236=wprcoef(wpt,[3,6]); s237=wprcoef(wpt,[3,7]); s20=norm(s230); s21=norm(s231); s22=norm(s232); s23=norm(s233); s24=norm(s234); s25=norm(s235); s26=norm(s236); s27=norm(s237); st20=std(s230); st21=std(s231); st22=std(s232); st23=std(s233); st24=std(s234); st25=std(s235); st26=std(s236); st27=std(s237); disp('故障信号的特征向量'); snorm2=[s20,s21,s22,s23,s24,s25,s26,s27] std2=[st20,st21,st22,st23,st24,st25,st26,st27] subplot(9,2,4);plot(s230); ylabel('S230'); subplot(9,2,6);plot(s231); ylabel('S231'); subplot(9,2,8);plot(s232); ylabel('S232'); subplot(9,2,10);plot(s233); ylabel('S233'); subplot(9,2,12);plot(s234); ylabel('S234'); subplot(9,2,14);plot(s235); ylabel('S235'); subplot(9,2,16);plot(s236); ylabel('S236'); subplot(9,2,18);plot(s237); ylabel('S237'); %fft figure y1=fft(s1,1024); py1=y1.*conj(y1)/1024; y2=fft(s2,1024); py2=y2.*conj(y2)/1024; y130=fft(s130,1024); py130=y130.*conj(y130)/1024; y131=fft(s131,1024); py131=y131.*conj(y131)/1024; y132=fft(s132,1024); py132=y132.*conj(y132)/1024; y133=fft(s133,1024); py133=y133.*conj(y133)/1024; y134=fft(s134,1024); py134=y134.*conj(y134)/1024; y135=fft(s135,1024); py135=y135.*conj(y135)/1024; y136=fft(s136,1024); py136=y136.*conj(y136)/1024; y137=fft(s137,1024); py137=y137.*conj(y137)/1024; y230=fft(s230,1024); py230=y230.*conj(y230)/1024; y231=fft(s231,1024); py231=y231.*conj(y231)/1024; y232=fft(s232,1024); py232=y232.*conj(y232)/1024; y233=fft(s233,1024); py233=y233.*conj(y233)/1024; y234=fft(s234,1024); py234=y234.*conj(y234)/1024; y235=fft(s235,1024); py235=y235.*conj(y235)/1024; y236=fft(s236,1024); py236=y236.*conj(y236)/1024; y237=fft(s237,1024); py237=y237.*conj(y237)/1024; f=1000*(0:511)/1024; subplot(1,2,1); plot(f,py1(1:512)); ylabel('P1'); title('原始信号的功率谱') subplot(1,2,2); plot(f,py2(1:512)); ylabel('P2'); title('故障信号的功率谱') figure subplot(4,2,1); plot(f,py130(1:512)); ylabel('P130'); title('S130的功率谱') subplot(4,2,2); plot(f,py131(1:512)); ylabel('P131'); title('S131的功率谱') subplot(4,2,3); plot(f,py132(1:512)); ylabel('P132'); subplot(4,2,4); plot(f,py133(1:512)); ylabel('P133'); subplot(4,2,5); plot(f,py134(1:512)); ylabel('P134'); subplot(4,2,6); plot(f,py135(1:512)); ylabel('P135'); subplot(4,2,7); plot(f,py136(1:512)); ylabel('P136'); subplot(4,2,8); plot(f,py137(1:512)); ylabel('P137'); figure subplot(4,2,1); plot(f,py230(1:512)); ylabel('P230'); title('S230的功率谱') subplot(4,2,2); plot(f,py231(1:512)); ylabel('P231'); title('S231的功率谱') subplot(4,2,3); plot(f,py232(1:512)); ylabel('P232'); subplot(4,2,4); plot(f,py233(1:512)); ylabel('P233'); subplot(4,2,5); plot(f,py234(1:512)); ylabel('P234'); subplot(4,2,6); plot(f,py235(1:512)); ylabel('P235'); subplot(4,2,7); plot(f,py236(1:512)); ylabel('P236'); subplot(4,2,8); plot(f,py237(1:512)); ylabel('P237'); figure %plottree(wpt)

二维小波变换（正和逆变换） 二维小波变换（正和逆变换），用二维卷积实现的，其中有个子函数用于产生小波系数的（小波分解和小波重构的系数都有）！你自己也可以往里面加系数！ 复制内容到剪贴板 代码:function [varargout]=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=[1 1]/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=[0:7];         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=[0:7];         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=[0:17];         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=[0:9];         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 [c,s]=wavefast(x,n,varargin) error(nargchk(3,4,nargin)); if nargin==3     if ischar(varargin{1})         [lp,hp]=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     [app,keep]=symextend(app,fl);     rows=symconv(app,hp,'row',fl,keep);     coefs=symconv(rows,hp,'col',fl,keep);     c=[coefs(:)' c]; s=[size(coefs);s];     coefs=symconv(rows,lp,'col',fl,keep);     c=[coefs(:)' c];     rows=symconv(app,lp,'row',fl,keep);     coefs=symconv(rows,hp,'col',fl,keep);     c=[coefs(:)' c];     app=symconv(rows,lp,'col',fl,keep); end c=[app(:)' c]; s=[size(app);s]; function [y, keep]=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 [varargout]=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(['[C S] 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)             [lp,hp]=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)             [lp,hp]=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=[a(:)' nc];     ns=ns(3:end,:);     ns=[ns(1,:);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([2 1].*size(x)); y(1:2:end,:)=x; y=conv2(y,f1'); z=zeros([1 2].*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([1:2:N-1],:);  %  奇数 f2=f([2:2:N],:);    %  偶数 % 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([1:1:T],:)=low_frequency_column([1:T],:); f_column([T+1:1:N],:)=high_frequency_column([1:T],:);         figure(1) colormap(map); image(f); figure(2) colormap(map); image(f_column); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%   2.行变换 %  A.分裂（奇偶分开） f1=f_column(:,[1:2:N-1]);  %  奇数 f2=f_column(:,[2:2:N]);    %  偶数 % 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(:,[1:1:T])=low_frequency_row(:,[1:T]); f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);     figure(3) colormap(map); image(f_row); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%   1.行变换 %   A.提取（低频高频分开） f1=f_row(:,[T+1:1:N]);  %  奇数 f2=f_row(:,[1:1:T]);    %  偶数 % 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(:,[2:2:N])=low_frequency_row(:,[1:T]); f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);     figure(4) colormap(map); image(f_row); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%   2.列变换 %  A.提取（低频高频分开） f1=f_row([T+1:1:N],:);  %  奇数 f2=f_row([1:1:T],:);    %  偶数 % 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([2:2:N],:)=low_frequency_column([1:T],:); f_column([1:2:N-1],:)=high_frequency_column([1:T],:);         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); [c,l]=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 [T_img,Sub_T]=threshold_2_N(img,levels) % reference :image denoising using wavelet thresholding [xx,yy]=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),[0,1]);figure;imshow(I1); [N,M] = size(I); h = [0.125,0.375,0.375,0.125]; g = [0.5,-0.5]; delta = [1,0,0]; 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)),[0 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)),[0 1]);figure;imshow(I1); end