小波分解与重构

我理解的小波分解是将一个多频率组成的波通过小波分解将所有频率分解出来，重构就是将这些分频率加起来得到最后的重构结果，于是写了个这样的程序
clc
close all;
clear all;
clc;
fs=612;
[reg,sta,data]=readmydata('beijing08.dat');
data{1:end};
A=ans(2:end);
for i=1:609;
if A(i)>50.0;
A(i)=(A(i-12)+A(i+12))/2;
end
end
for i=609:612;
if A(i)>50.0;
A(i)=(A(i-12)+A(i-24))/2;
end
end
%信号时域波形
figure(1);
plot(1:612,A);
%使用db5小波进行尺度为7时的分解
[c,l]=wavedec(A,9,'db5');
%从小波分解结构[c,l]重构信号xdata
a0=waverec(c,l,'db5');
%检查重构效果
figure(2);
subplot(3,1,1);
plot(A);
title('原始信号')
subplot(3,1,2);
plot(a0);
title('重构信号')
subplot(3,1,3);
plot(A-a0);
title('误差信号')
err=max(abs(A-a0))
%重构第1~5层高频细节信号
d9=wrcoef('d',c,l,'db5',9);
d8=wrcoef('d',c,l,'db5',8);
d7=wrcoef('d',c,l,'db5',7);
d6=wrcoef('d',c,l,'db5',6);
d5=wrcoef('d',c,l,'db5',5);
d4=wrcoef('d',c,l,'db5',4);
d3=wrcoef('d',c,l,'db5',3);
d2=wrcoef('d',c,l,'db5',2);
d1=wrcoef('d',c,l,'db5',1);
%显示高频细节信号
figure(3);
subplot(9,1,1);
plot(d9,'LineWidth',2);
ylabel('d9');
subplot(9,1,2);
plot(d8,'LineWidth',2);
ylabel('d8');
subplot(9,1,3);
plot(d7,'LineWidth',2);
ylabel('d7');
subplot(9,1,4);
plot(d6,'LineWidth',2);
ylabel('d6');
subplot(9,1,5);
plot(d5,'LineWidth',2);
ylabel('d5');
subplot(9,1,6);
plot(d4,'LineWidth',2);
ylabel('d4');
subplot(9,1,7);
plot(d3,'LineWidth',2);
ylabel('d3');
subplot(9,1,8);
plot(d2,'LineWidth',2);
ylabel('d2');
xlabel('时间 t/s');
subplot(9,1,9);
plot(d1,'LineWidth',2);
ylabel('d1');
%第1层高频细节信号的包络谱
y=hilbert(d1);
ydata=abs(y);
y=y-mean(y);
nfft=1024;
p=abs(fft(ydata,nfft));
figure(4);
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
xlabel('频率 f/Hz');
ylabel('功率谱 P/W');

小波分解与重构程序
>> clear
I=imread('C:\Documents and Settings\Administrator\桌面\暑期/cidian.bmp');
I=rgb2gray(I);
[X,map]=gray2ind(I);
subplot(2,2,1);
imshow(X,map);
title('原始图像');
X=double(X);
sX=size(X);
[cA,cH,cV,cD]=dwt2(X,'db4');
A0=idwt2(cA,cH,cV,cD,' db4', sX);
subplot(2,2,2);
imshow(A0,map);
title('db4小波重构');
error1=max(max(abs(X-A0)))


程序很简单，也很基础。不要以为你见过，其实仔细看，有一些改进哦。

clc;
clear;
% 装载图像
load woman;
% X包含载入的图像
% 绘制原始图像
figure(1);
subplot(2,2,1);
image(X);
colormap(map);
title('原始图像');
% 使用sym5对X进行尺度为2的分解
[c,s] = wavedec2(X,2,'sym5');
% 从小波分解结构[c,s]进行尺度为1和2时的低频重构
a1 = wrcoef2('a',c,s,'sym5',1);
a2 = wrcoef2('a',c,s,'sym5',2);

% 绘制尺度为1时的低频图像
subplot(2,2,3);
image(a1);colormap(map);
title('尺度为1时的低频图像');
% 绘制尺度为2时

的低频图像
subplot(2,2,4);
image(a2);colormap(map);
title('尺度为2时的低频图像');
% 从小波分解结构[c,s]在尺度为2时重构高频
% 'h' 是水平方向
% 'v' 是垂直方向
% 'd' 是对角方向
hd2 = wrcoef2('h',c,s,'sym5',2);
vd2 = wrcoef2('v',c,s,'sym5',2);
dd2 = wrcoef2('d',c,s,'sym5',2);
% 绘制高频图像
figure(2);
subplot(2,2,1);
image(hd2);colormap(map);
title('尺度为2时的水平高频图像');
subplot(2,2,2);
image(vd2);colormap(map);
title('尺度为2时的垂直高频图像');
subplot(2,2,3);
image(dd2);colormap(map);
title('尺度为2时的对角高频图像');
subplot(2,2,4);
image(hd2+dd2+vd2+a1);colormap(map);
% 验证这些图像的长度都是sX
sX = size(X)
sa1 = size(a1)
shd2 = size(hd2)
norm(hd2+dd2+vd2+a1-X)



结果：
sX =
256 256

sa1 =
256 256

shd2 =
256 256

ans =
8.106902321771278e-010
分析：
由ans的值可以看到，重构的效果几乎完美！


小波变换的模极大值重建信号源程序
这是用小波变换模极大值重建信号的源程序，数据是一心电信号，在matlab6.5下实现，来源于胡广书的《现代信号处理教程》附属光盘程序包含四个部份，其中wave_peak实现信号的分解并求出模极大值%--------------------------------------------------------------------------
% exa130301.m 例13.3.1: 利用小波变换模极大重建原信号
%--------------------------------------------------------------------------
close all;
points=1024; level=6; sr=360; num_inter=6; wf='db3';
%所处理数据的长度 分解的级数 抽样率 迭代次数 小波名称

offset=0;
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);

%计算小波分解系数和模极大序列
[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset);
% signal: 原始信号; swa:小波概貌; swd:小波细节;
% ddw: 局部极大位置; wpeak:小波变换的局部极大序列]

pswa=swa(level,: ); % pswa: 为待重建的信号
wframe=(wpeak~=0);
%迭代初始化
w0=zeros(1,points);
[a,d]=swt(w0,level,Lo_D,Hi_D);
w2=d; % w2为待重建小波
for j=1:num_inter
w2=Py_Pgama(d,wpeak,wframe,1,sr); % 先进行Py投影和 Pgama投影
w0=iswt(pswa,w2,Lo_R,Hi_R); % 再进行Pv投影
[a,d]=swt(w0,level,Lo_D,Hi_D); % Pv
end
pswa=iswt(swa(level,: ),w2,Lo_R,Hi_R); % 计算重建信号

% 原信号和由模极大重建信号的比较
figure,
subplot(211)
plot(pswa(1:points));
subplot(212)
plot(signal(1:points),'r');

%分别计算重建小波以及原信号的信噪比
werr=w2-swd;
% 原信号的小波变换（swd)和重建后的小波变换（w2）的比较
figure,
for m=1:level
wsnr(m)=20*log10(norm(swd(m,: ))/norm(werr(m,: )))
subplot(level+1,1,m);
plot(swd(m,: )),hold on,
plot(w2(m,: ),'r');grid on;ylabel(strcat('j=',num2str(m))),axis tight;
end
err=ps

wa(1:points)-signal(1:points);
snr=20*log10(norm(signal)/norm(err));


小波变换中能量泄漏问题的研究
程序代码：
clear all;clc;
f1=16.7;% 频率1
f2=33.4;
f3=50.1;
Fs=1000; % 采样频率
Ts=1/Fs; % 采样间隔
N=1024; % 采样点数
t=[0:1/Fs:(N-1)/Fs]; %采样时刻
lev=5;
x=sin(2*pi*f1*t);
figure(1);
plot(t,x);
%%小波分解
[C,L]=wavedec(x,lev,'db10'); %C由[cAj,cDj,cDj-1,...,cD1]
%提取多尺度小波变换的高频系数
figure(2);
lev_1=lev+1;
for i=1:lev
cD=detcoef(C,L,i);
subplot(lev_1,1,i)
plot(cD);
ylabel(['cD',num2str(i)]);
title(['Detail cD',num2str(i)]);
end
%提取多尺度小波变换的低频系数
figure(2);

cA=appcoef(C,L,'db10',lev);
subplot(lev_1,1,lev_1)
plot( cA);
ylabel(['cA',num2str(lev)]);
title(['Approximation cA',num2str(lev)])

%%小波重构
figure(3);
for i=1:lev
D=wrcoef('d',C,L,'db10',i);
subplot(lev_1,1,i)
plot(D);
ylabel(['D',num2str(i)]);
title(['Detail D',num2str(i)]);
end

figure(3);

A=wrcoef('a',C,L,'db10',lev);
subplot(lev_1,1,lev_1)
plot( A);
ylabel(['A',num2str(lev)]);
title(['Approximation A',num2str(lev)]);

%%计算能量谱
figure(4);
[Ea,Ed] = wenergy(C,L);
EE1=sum(Ed)+Ea;
for i=1:lev
Ed1(lev_1-i)= Ed(i);
end

E1=[Ea Ed1]/EE1;
x2=1:lev_1; %p=rand(6,1);
bar(x2,E1);

str = cell(1,lev_1); %建立单元数组1*10的空数组
for i = 0:lev
str{i+1} = ['E',num2str(i)]; %num2str(i)数字转为字符串输出。
end
set(gca,'xticklabel',str);


求Matlab小波能量谱的程序
给你这个程序参考下：
figure(1)
fs=2000;
t=0:1/fs:1;
for i=1:length(t)
if t(i)<=0.2
x(i)=8*sin(60*pi*t(i)+pi/12)+6*sin(240*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6);
else
x(i)=8*sin(60*pi*t(i)+pi/12)+6*sin(240*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6)+20*sin(500*pi*t(i)+pi/3).*exp(-20*(t(i)-0.2));
end
end
plot(t,x)
grid on
figure(2)
wpt=wpdec(x,4,'db3');
plot(wpt);
for i=0:15
figure(3);
rfs=wprcoef(wpt,[4 i]);
subplot(4,4,i+1);
plot(rfs);
grid on;
title(['[4 ',num2str(i),']']);
axis([0 2000 min(rfs) max(rfs)]);
end
for i=0:3
figure(4);
rfs=wprcoef(wpt,[4 i]);
subplot(4,1,i+1);
plot(rfs);
grid on;
title(['[4 ',num2str(i),']']);
axis([0 2000 min(rfs) max(rfs)]);
end
for i=1:4
wpt=wpdec(x,i,'db3');
e=wenergy(wpt);
E=zeros(1,length(e));
for j=1:2^i
E(j)=sum(abs(wprcoef(wpt,[i,j-1])).^2);
end
figure(5)
subplot(4,1,i);
bar(e);
axis([0 length(e) 0 130]);
title(['第 ',num2str(i), ' 层']);
for j=1:length(e)
text(j-0.2,e(j)+20,num2str(e(j),'%2.2f'));
end
end



0.0000 10.8567 26.5764
0.0005 9.4834 21.235

8
0.0010 8.4153 18.0315
0.0015 6.5842 15.2849
0.0020 5.6687 13.4539
0.0025 5.3635 10.4021
0.0030 5.6687 8.8762
0.0035 4.4480 8.2659
0.0040 2.9221 6.1296
0.0045 3.9902 6.2822
0.0050 2.9221 5.2141
0.0055 3.0747 5.0615
0.0060 2.0066 3.9934
0.0065 2.0066 4.2986
以上是数据形式，是一个txt文件，要求第一列和第二列做一个图，第一列和第三列一个图，然后分别进行小波分解与重构，下面是程序但总有错，望各位高手帮帮忙！
[c11 c22 c33] =textread('2ms.txt','%f%f%f');
s=[c11, c22];
subplot(2,2,1);
plot(s);
title('原始信号');
%=============================
%用db1小波对原始信号进行3层分解并提取系数
[c,l]=wavedec(s,3,'db1');
a3=appcoef(c,l,'db1',3);
d3=detcoef(c,l,3);
d2=detcoef(c,l,2);
d1=detcoef(c,l,1);
%=============================
%对信号进行强制性消噪处理并图示结果
dd3=zeros(1,length(d3));
dd2=zeros(1,length(d2));
dd1=zeros(1,length(d1));
c1=[a3 dd3 dd2 dd1];
s1=waverec(c1,l,'db1');
subplot(2,2,2);
plot(s1);grid;
title('强制消噪后的信号');
%=============================
%用默认阈值对信号进行消噪处理并图示结果
%用ddencmp函数获得信号的默认阈值
[thr,sorh,keepapp]=ddencmp('den','wv',s);
s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);
subplot(2,2,3);
plot(s2);grid;
title('默认阈值消噪后的信号');
%=============================
%用给定的软阈值进行消噪处理
softd1=wthresh(d1,'s',1.465);
softd2=wthresh(d2,'s',1.823);
softd3=wthresh(d3,'s',2.768);
c2=[a3 softd3 softd2 softd1];
s3=waverec(c2,l,'db1');
subplot(2,2,4);
plot(s3);grid;
title('给定软阈值消噪后的信号');
问题补充：c1=[a3 dd3 dd2 dd1]; 这里提示有错

c1=[a3‘ dd3’ dd2‘ dd1’];


Matlab语音信号小波分解重构

softd1=wthresh(d1,'s',1.465);
softd2=wthresh(d2,'s',1.823);
softd3=wthresh(d3,'s',2.768);
c2=[a3,softd3,softd2,softd1];
s3=waverec(c2,l,'db1');
语音信号小波分解重构我将信号分解后 用软阈值进行降噪 得到降噪后的低频系数
但是在重构时 c2=[a3,softd3,softd2,softd1];
这句出现错误（我将逗号换成空格也是如此,我的版本是2009a)
希望高手指点

用db1小波对原始信号进行3层分解并提取系数[c,l]=wavedec(s,3,'db1');a3=appcoef(c,l,'db1',3);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);%=============================%对信号进行强制性消噪处理并图示结果dd3=zeros(1,length(d3));dd2=zeros(1,length(d2));dd1=zeros(1,length(d1));c1=[a3 dd3 dd2 dd1];s1=waverec(c1,l,'db1');subplot(2,2,2);plot(s1);grid;title('强制消噪后的信号');%=============================%用默认阈值对信号进行消噪处理并图示结果 %用ddencmp函数获得信号的默认阈值[thr,sorh,keepapp]=ddencmp('den','wv',s);s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);subplot(2,2,3);plot(s2);grid;title('默认阈值消噪后

的信号');%=============================%用给定的软阈值进行消噪处理softd1=wthresh(d1,'s',1.465);softd2=wthresh(d2,'s',1.823);softd3=wthresh(d3,'s',2.768);c2=[a3 softd3 softd2 softd1];s3=waverec(c2,l,'db1');subplot(2,2,4);plot(s3);grid;title('给定软阈值消噪后的信号')

用*.mat装在信号可以实现
但是我装在*.wav后 就无法实现 说c2=[a3 softd3 softd2 softd1]句中有错误