小波去噪软阈值和硬阈值的matlab仿真程序
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
%设置信噪比和随机种子值
snr=4;
init=2055615866;
%产生原始信号sref和高斯白噪声污染的信号s
[sref,s]=wnoise(1,11,snr,init);
%用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);
thr=1;
%进行硬阈值处理
ythard1=wthresh(d1,'h',thr);
ythard2=wthresh(d2,'h',thr);
ythard3=wthresh(d3,'h',thr);
c2=[a3 ythard3 ythard2 ythard1];
s3=waverec(c2,l,'db1');
%进行软阈值处理
ytsoftd1=wthresh(d1,'s',thr);
ytsoftd2=wthresh(d2,'s',thr);
ytsoftd3=wthresh(d3,'s',thr);
c3=[a3 ytsoftd3 ytsoftd2 ytsoftd1];
s4=waverec(c3,l,'db1');
%对上述信号进行图示
subplot(5,1,1);plot(sref);title('参考信号');
subplot(5,1,2);plot(s);title('染噪信号');
subplot(5,1,3);plot(s3);title('硬阈值处理');
subplot(5,1,4);plot(s4);title('软阈值处理');
%设置信噪比和随机种子值
snr=4;
init=2055615866;
%产生原始信号sref和高斯白噪声污染的信号s
[sref,s]=wnoise(1,11,snr,init);
%用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);
thr=1;
%进行硬阈值处理
ythard1=wthresh(d1,'h',thr);
ythard2=wthresh(d2,'h',thr);
ythard3=wthresh(d3,'h',thr);
c2=[a3 ythard3 ythard2 ythard1];
s3=waverec(c2,l,'db1');
%进行软阈值处理
ytsoftd1=wthresh(d1,'s',thr);
ytsoftd2=wthresh(d2,'s',thr);
ytsoftd3=wthresh(d3,'s',thr);
c3=[a3 ytsoftd3 ytsoftd2 ytsoftd1];
s4=waverec(c3,l,'db1');
%对上述信号进行图示
subplot(5,1,1);plot(sref);title('参考信号');
subplot(5,1,2);plot(s);title('染噪信号');
subplot(5,1,3);plot(s3);title('硬阈值处理');
subplot(5,1,4);plot(s4);title('软阈值处理');
Matlab小波消噪程序(原创)
[Yo,FS,NBITS,OPTS]=wavread('normal_s3.wav');
fs=FS
Y=Yo;
dt=1/FS;
YY=Yo(1800:9500)';
N=length(YY);
Y=YY+0.14*randn(1,N);
time=(0:N-1)*dt;
(N-1)*dt
fY=fs*(1:N/2)/N;
yYY=abs(fft(YY));
yY=abs(fft(Y));
[c,l]=wavedec(Y,4,'coif4');
a4=appcoef(c,l,'coif4',4);
a3=appcoef(c,l,'coif4',3);
a2=appcoef(c,l,'coif4',2);
a1=appcoef(c,l,'coif4',1);
d4=detcoef(c,l,4);
d3=detcoef(c,l,3);
d2=detcoef(c,l,2);
d1=detcoef(c,l,1);
figure(1)
subplot(511);plot(a4);title('a4');
subplot(512);plot(d4);title('d4');
subplot(513);plot(d3);title('d3');
subplot(514);plot(d2);title('d2');
subplot(5,1,5);plot(d1);title('d1');
SNRin=10*log(norm(YY)^2/(norm(Y-YY))^2) %去噪前信号的信噪比
figure(2)
subplot(211);plot(time,YY);axis([0 time(end) -1.5 1.5]);xlabel('time/s');title('原音信号');
subplot(212);plot(fY,yYY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('原心音信号频谱图');xlim([0 500])
figure(3)
subplot(211);plot(time,Y);axis([0 time(end) -1.5 1.5]);xlabel('time/s');title('加噪心音信号');
subplot(212);plot(fY,yY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('
加噪心音信号频谱');xlim([0 500])
%%%%%%%%%%%单级阈值去噪%%%%%%%%%%%%%%%
thr=yuzhi(Y,4,'coif4','sqtwolog','sln')
dd1=yuzhifunction(d1,thr(1),'soft');
dd2=yuzhifunction(d2,thr(2),'soft');
dd3=yuzhifunction(d3,thr(3),'soft');
dd4=yuzhifunction(d4,thr(4),'soft');
cl=[a4,dd4,dd3,dd2,dd1];
s=waverec(cl,l,'coif4');
N1=length(s);
times=(0:N1-1)*dt;
fs0=fs*(0:N1/2-1)/N1;
ys=abs(fft(s));
%%%%%%%%%%%多级阈值去噪%%%%%%%%%%%%%%
thr1=yuzhi(Y,4,'coif4','sqtwolog','mln')
[s1,cxc,clc,perf2,perf3]=wdencmp('lvd',Y,'coif4',4,thr1,'s');
N2=length(s1);
times1=(0:N2-1)*dt;
fs1=fs*(1:N2/2)/N2;
ys1=abs(fft(s1));
figure(4);
subplot(321);plot(time,YY);axis([0 time(end) -1.5 1.5]);xlabel('time/s');title('原音信号');
subplot(322);plot(fY,yYY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('原心音信号频谱图');xlim([0 400])
subplot(323);plot(times,s);axis([0 time(end) -1.5 1.5]);xlabel('time/s');title('单级阈值去噪心音信号');
subplot(324);plot(fs0,ys(1:floor(N1/2)));xlabel('f/hz');ylabel('频谱');title('单级阈值去噪心音信号频谱图');xlim([0 400])
subplot(325);plot(times1,s1);axis([0 time(end) -1.5 1.5]);xlabel('time/s');title('多级阈值去噪心音信号');
subplot(326);plot(fs1,ys1(1:floor(N2/2)));xlabel('f/hz');ylabel('频谱');title('多级阈值去噪心音信号频谱图');xlim([0 400])
想把这个程序。改为只用软硬阈值对比的心电信号去噪分析,该怎么办? 求高手。程序见下:
%%%%%%%%%%信号小波分解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%基于Haar小波%%%%%%%
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10000;
save ECGdata data;
fclose(ecg);
x=0:0.004:4.8;
figure(1);
subplot(321);
plot(x,data);
axis([0 4.8 -4 4]);
title('原始心电信号');
x=data;
wname='Haar';
level=5;
[c,l]=wavedec(x,level,wname);
a5=wrcoef('a',c,l,'Haar',5);
a4=wrcoef('a',c,l,'Haar',4);
a3=wrcoef('a',c,l,'Haar',3);
a2=wrcoef('a',c,l,'Haar',2);
a1=wrcoef('a',c,l,'Haar',1);
subplot(322);plot(a5);title('a5');axis([0 1201 -2 2]);
subplot(323);plot(a4);title('a4');axis([0 1201 -2 2]);
subplot(324);plot(a3);title('a3');axis([0 1201 -2 2]);
subplot(325);plot(a2);title('a2');axis([0 1201 -2 2]);
subplot(326);plot(a1);title('a1');axis([0 1201 -2 2]);
%%%%%%%%%%%%基于db6小波%%%%%
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10;
fclose(ecg);
x=0:0.004:4.8;
figure(2);
subplot(321);
plot(x,data);
axis([0 4.8 -4000 4000]);
title('原始心电信号');
x=data;
wname='db6';
level=5;
[c,l]=wavedec(x,level,wname);
a5=wrcoef('a',c,l,'db6',5);
a4=wrcoef('a',c,l,'db6',4);
a3=wrcoef('a',c,l,'db6',3);
a2=wrcoef('a',c,l,'db6',2);
a1=wrcoef('a',c,l,'db6',1);
subplot(322);plot(a5);title('a5');axis([0 1201 -2000 2000]);
subplot(323);plot(a4);title('a4');axis([0 1201 -2000 2000]);
subplot(324);plot(a3);title('a3');axis([0 1201 -2000 2000]);
subplot(325);plot(a2);title('a2');axis
([0 1201 -2000 2000]);
subplot(326);plot(a1);title('a1');axis([0 1201 -2000 2000]);
%%%%%%%%%%%%%%%%%基于sym3小波%%%%%%%%%%%%%%%%%%%%%%%
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10;
save ECGdata data;
fclose(ecg);
x=0:0.004:4.8;
figure(3);
subplot(321);
plot(x,data);
axis([0 4.8 -4000 4000]);
title('原始心电信号');
x=data;
wname='sym3';
level=5;
[c,l]=wavedec(x,level,wname);
a5=wrcoef('a',c,l,'sym3',5);
a4=wrcoef('a',c,l,'sym3',4);
a3=wrcoef('a',c,l,'sym3',3);
a2=wrcoef('a',c,l,'sym3',2);
a1=wrcoef('a',c,l,'sym3',1);
subplot(322);plot(a5);title('a5');axis([0 1201 -2000 2000]);
subplot(323);plot(a4);title('a4');axis([0 1201 -2000 2000]);
subplot(324);plot(a3);title('a3');axis([0 1201 -2000 2000]);
subplot(325);plot(a2);title('a2');axis([0 1201 -2000 2000]);
subplot(326);plot(a1);title('a1');axis([0 1201 -2000 2000]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%心电信号降噪
%%%%%%%%%%%%%%%Birge-Massart策略阈值降噪
%基于小波变换的心电信号的降噪
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10000;
fclose(ecg);
x=data;
wavename='db5';
level=4;
[c,l]=wavedec(x,level,wavename);
alpha=1.5;
sorh='h';
[thr,nkeep]=wdcbm(c,l,alpha);%使用硬阈值给信号降噪
[xc,cxc,lxc,perf0,perfl2]=wdencmp('lvd',c,l,wavename,level,thr,sorh);
t1=0:0.004:(length(x)-1)*0.004;
figure(4);
subplot(211);
plot(t1,x);
title('从人体采集的原始的ECG信号');
subplot(212);
plot(t1,xc);
title('Birge-Massart策略阈值降噪后的ECG信号(wname=db5 level=4)');
%%%%%%%%%%%%%%%%%%最优预测软阈值降噪
%基于小波变换的心电信号的压缩
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10000;
fclose(ecg);
x=data;
wavename='db3';
level=4;
[xd,cxd,lxd]=wden(x,'heursure','s','mln',level,wavename);
t=0:0.004:4.8;
figure(5);
subplot(311);
plot(t,x);
title('从人体采集的原始的ECG信号');
subplot(312);
plot(t,xd);
title('heursure规则阈值降噪后的ECG信号(db3 level=4)');
%%%%%%%%%%%%%%最小均方差软阈值降噪
%基于小波变换的心电信号的压缩
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10000;
fclose(ecg);
x=data;
wavename='db3';
level=4;
xd=wden(x,'minimaxi','s','mln',level,wavename);
t=0:0.004:4.8;
figure(6);
subplot(311);
plot(t,x);
title('从人体采集的原始的ECG信号');
subplot(313);
plot(t,xd);
title('Minimax规则阈值降噪后的ECG信号(wname=db3 level=4)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%心电信号压缩
%基于小波变换的心电信号的压缩
ecg=fopen('100.dat','r');
N=1201;
data=fread(ecg,N,'int16');
data=data/10000;
fclose(ecg);
x=data;
wavename='db5';
[c,l]=wavedec(x,4,wavename);
alpha=2;
sorh='h';
[thr,nkeep]=wdcbm(c,l,alpha);%使用硬阈值压缩信号
[xc,cxc,lxc,perf0,perfl2]=w
dencmp('lvd',c,l,wavename,4,thr,sorh);
t=0:0.004:4.8;
figure(7);
subplot(511);
plot(t,x);
title('原始的ECG信号');
axis([0 5 -3 3]);
subplot(512);
plot(t,xc);
axis([0 5 -3 3]);
title('在第4层对高频系数量化压缩后的ECG信号(alpha=2 db5)');
[c,l]=wavedec(x,4,wavename);
alpha=4;
sorh='h';
[thr,nkeep]=wdcbm(c,l,alpha);%使用硬阈值压缩信号
[xc1,cxc1,lxc1,perf01,perfl21]=wdencmp('lvd',c,l,wavename,4,thr,sorh);
subplot(513);
plot(t,xc1);
axis([0 5 -3 3]);
title('在第4层对高频系数量化压缩后的ECG信号alpha=4 db5');
[c,l]=wavedec(x,2,wavename);
alpha=4;
sorh='h';
[thr,nkeep]=wdcbm(c,l,alpha);%使用硬阈值压缩信号
[xc1,cxc1,lxc1,perf01,perfl21]=wdencmp('lvd',c,l,wavename,2,thr,sorh);
subplot(514);
plot(t,xc1);
axis([0 5 -3 3]);
title('在第2层对高频系数量化压缩后的ECG信号alpha=4 db5');
[c,l]=wavedec(x,2,'db3');
alpha=4;
sorh='h';
[thr,nkeep]=wdcbm(c,l,alpha);%使用硬阈值压缩信号
[xc1,cxc1,lxc1,perf01,perfl21]=wdencmp('lvd',c,l,'db3',2,thr,sorh);
subplot(515);
plot(t,xc1);
axis([0 5 -3 3]);
title('在第2层对高频系数量化压缩后的ECG信号alpha=4 db3');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Matlab小波去噪(默认,强制,给定三种情况)
%% 利用小波分析对监测采集的信号进行去噪处理,恢复原始信号
%小波分析进行去噪有3中方法:
%1、默认阈值去噪处理。该方法利用函数ddencmp( )生成信号的默认阈值,然后利用函数wdencmp( )进行去噪处理;
%2、给定阈值去噪处理。在实际的去噪处理过程中,阈值往往可通过经验公式获得,且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh( );
%3、强制去噪处理。该方法是将小波分解结构中的高频系数全部置0,即滤掉所有高频部分,然后对信号进行小波重构。这种方法比较简单,且去噪后的信号比较平滑,但是容易丢失信号中的有用成分。
%% 利用小波分析对监测采集的水轮机信号进行去噪处理,恢复原始信号
%Program Start
%% 载入监测所得信号
load default.txt; %装载采集的信号
x= default;
lx=length(x);
t=[0:1:length(x)-1]';
%% 绘制监测所得信号
subplot(2,2,1);
plot(t,x);
title('原始信号');
grid on
set(gcf,'color','w')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
%% 用db1小波对原始信号进行3层分解并提取小波系数
[c,l]=wavedec(x,3,'db1');%sym8
ca3=appcoef(c,l,'db1',3);%低频部分
cd3=detcoef(c,l,3);%高频部分
cd2=detcoef(c,l,2);%高频部分
cd1=detcoef(c,l,1);%高频部分
%% 对信号进行强制去噪处理并图示
cdd3=zeros(1,length(cd3));
cdd2=zeros(1,length(cd2));
cdd1=zeros(1,length(cd1));
c1=[ca3,cdd3,cdd2,cdd1];
x1=waverec(c1,1,'db1');
subplot(2,2,2);
plot(x1);
title('强制去噪后信号');
grid on
set(gcf,'color
','w')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
%% 默认阈值对信号去噪并图示
%用ddencmp( )函数获得信号的默认阈值,使用wdencmp( )函数实现去噪过程
[thr,sorh,keepapp]=ddencmp('den','wv',x);
x2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);
subplot(2,2,3);
plot(x2);
title('默认阈值去噪后信号');
grid on
set(gcf,'color','w')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
%% 给定的软阈值进行去噪处理并图示
cd1soft=wthresh(cd1,'x',1.465);%经验给出软阈值数
cd2soft=wthresh(cd2,'x',1.823); %经验给出软阈值数
cd3soft=wthresh(cd3,'x',2.768); %经验给出软阈值数
c2=[ca3,cd3soft,cd2soft,cd1soft];
x3=waverec(c2,1,'db1');
subplot(2,2,4);
plot(x3);
title('给定软阈值去噪后信号');
grid on
set(gcf,'color','w')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
以上就是三中小波去噪的原程序。
但红色标注的地方,我却运行不过去。
提示分别问:
1,??? Error using ==> horzcat
CAT arguments dimensions are not consistent.
Error in ==> tt at 33
c1=[ca3,cdd3,cdd2,cdd1];
2,??? Error using ==> wthresh at 30
Invalid argument value.
Error in ==> tt at 54
cd1soft=wthresh(cd1,'x',1.465);%经验给出软阈值数
请问为什么???
我解决这个问题了,c1=[a3' dd3 dd2 dd1]
c2=[ca3' cd3soft' cd2soft' cd1soft'];
写一个小波去噪的matlab程序
sound=wavread('000');sound=sound';
subplot(4,1,1);plot(sound);title('原始信号');grid;set(gca,'ytick',-1:1:1)
count1=length(sound);
noise=0.07*randn(1,count1);%count1=28000
for i=1:count1
signal(i)=sound(i)+noise(i);
end;
for i=1:count1
y(i)=signal(i);
end;
subplot(4,1,2); plot(y);
title('含噪信号');grid;set(gca,'ytick',-1:1:1)
%读出带噪语音信号,存为‘101.wav'
wavwrite(y,8000,16,'c101');
SNR1 = 10*log10(sum(sound.^2)/sum((y-sound).^2)); %计算输入信噪比
title(['Noise Added (SNR=',num2str(SNR1),'dB)']); % num2str把数字转换成字符串
[C,L]=wavedec(y,3,'db3'); %将信号y分解为3级
ca=appcoef(C,L,'db3',3); %提取一维近似系数,3表示级数,db3为字符串小波的名字
[cd1,cd2,cd3]=detcoef(C,L,[1,2,3]); %3层小波分解
one=mad(cd1); % mad计算样本数据平均绝对偏差
two=mad(cd2); %30---38为阈值函数设置
three=mad(cd3);
n1=one/0.6745;
n2=two/0.6745;
n3=three/0.6745;
men1=n1*sqrt(2*log(count1))/5; %固定阈值形式,阈值为与√2*信号的长度,其中count1为信号长度
men2=n2*sqrt(2*log(count1))/5;
men3=n3*sqrt(2*log(count1))/5;
E1=sum((abs(cd1)).^2)/length(cd1);
E2=sum((abs(cd2)).^2)/length(cd2);
E3=sum((abs(cd3)).^2)/length(cd3);
E=[1,E3/E1,E2/E1];
Em=max(E);
z=E3/E1;
if ((Em==1)&&(E3/E1<0.8))
f=1;
cd1soft=wthresh(cd1,'s',men1); %其中S表示软阈值,如果是H则表示硬阈值
cd2soft=wthresh(cd2,'s',men2);
cd3soft=wthresh(cd3,'s',men3);
c1=[ca,cd3s
oft,cd2soft,cd1soft];
else
f=0;
casoft=wthresh(ca,'s',0.0050);
cd1soft=wthresh(cd1,'s',men1);
cd2soft=wthresh(cd2,'s',men2);
cd3soft=wthresh(cd3,'s',men3);
c1=[casoft,cd3soft,cd2soft,cd1soft];
end;
t=f;
s3=waverec(c1,L,'db3');
s=s3;
subplot(4,1,3);plot(s);title('消噪信号');grid;set(gca,'ytick',-1:1:1)
SNR2 = 10*log10(sum(sound.^2)/sum((s-sound).^2));
wavwrite(s,8000,16,'c102');
求小波去噪软阈值和硬阈值的matlab仿真程序
硬阈值、软阈值
这里有一段不知道有用没
%设置信噪比和随机种子值
snr=4;
init=2055615866;
%产生原始信号sref和高斯白噪声污染的信号s
[sref,s]=wnoise(1,11,snr,init);
%用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);
thr=1;
%进行硬阈值处理
ythard1=wthresh(d1,'h',thr);
ythard2=wthresh(d2,'h',thr);
ythard3=wthresh(d3,'h',thr);
c2=[a3 ythard3 ythard2 ythard1];
s3=waverec(c2,l,'db1');
%进行软阈值处理
ytsoftd1=wthresh(d1,'s',thr);
ytsoftd2=wthresh(d2,'s',thr);
ytsoftd3=wthresh(d3,'s',thr);
c3=[a3 ytsoftd3 ytsoftd2 ytsoftd1];
s4=waverec(c3,l,'db1');
%对上述信号进行图示
subplot(5,1,1);plot(sref);title('参考信号');
subplot(5,1,2);plot(s);title('染噪信号');
subplot(5,1,3);plot(s3);title('硬阈值处理');
subplot(5,1,4);plot(s4);title('软阈值处理');
求助Matlab小波语音去噪程序修改
诸位我编了个小波语音去噪的程序,可是分帧清浊音判断后的效果不佳,不进行清浊音判断直接阈值效果反而好,能否修改一下,谢谢
sound=wavread('c12345');
count1=length(sound);
noise=0.05*randn(1,count1);
for i=1:count1
y(i)=sound(i)+noise(i);
end
figure(1)
plot(sound);
title('原始信号');
set(gco,'linewidth',1);
set(gca,'fontname','宋体','fontsize',10);
xlabel('采样点','fontname','宋体','fontsize',10);
ylabel('幅值','fontname','宋体','fontsize',10);
figure(2)
plot(y);
title('含噪信号');
set(gco,'linewidth',1);
set(gca,'fontname','宋体','fontsize',10);
xlabel('采样点','fontname','宋体','fontsize',10);
ylabel('幅值','fontname','宋体','fontsize',10);
%在小波基'db3'下进行一维离散小波变换
N=250;
n=fix(count1/N);
t=[];
s = [0 0];
for m=1:n
xd=y(N*(m-1)+1:N+N*(m-1));
[C,L]=wavedec(xd,4,'sym4');
ca=appcoef(C,L,'sym4');
[cd1,cd2,cd3,cd4]=detcoef(C,L,[1,2,3,4]);
one=median(abs(cd1));
two=median(abs(cd2));
three=median(abs(cd3));
four=median(abs(cd4));
five=median(abs(ca));
n1=one/0.6745;
n2=two/0.6745;
n3=three/0.6745;
n4=four/0.6745;
n5=five/0.6745;
men1=n1*sqrt(2*log(count1))/log(1+1);
men2=n2*sqrt(2*log(count1))/log(2+1);
men3=n3*sqrt(2*log(count1))/log(3+1);
men4=n4*sqrt(2*log(count1))/log(4+1);
men5=n5*sqrt(2*log(count1))/log(4+1);
E1=sum((abs(cd1)).^2)/length(cd1);
E2=sum((abs(cd2)).^2)/length(cd2);
E3=sum((abs(cd3)).^2)/length(cd3);
E4=sum((abs(cd4)).^2)/length(cd4);
E=[1,E4/E1,E3/E1,E2/E1];
Em=max(E);
z=E4/E1;
if ((Em==1)&&(E4/E1<0.9))
f=1;
cd1soft=wthresh(cd1,'h',men1);
cd2soft=wthresh(cd2,'h',men2);
cd3soft=wthresh(cd3,'h',men3);
cd4soft=wthresh(cd4,'h',men4);
c1=[ca,cd4soft,cd3soft,cd2soft,cd1soft];
else
f=0;
casoft=wthresh(ca,'h',men5);
cd1soft=wthresh(cd1,'h',men1);
cd2soft=wthresh(cd2,'h',men2);
cd3soft=wthresh(cd3,'h',men3);
cd4soft=wthresh(cd4,'h',men4);
c1=[casoft,cd4soft,cd3soft,cd2soft,cd1soft];
end;
t=f
s3=waverec(c1,L,'sym4');
s=[s s3];
figure(3)
plot(s);
title('去噪后信号');
set(gco,'linewidth',1);
set(gca,'fontname','宋体','fontsize',10);
xlabel('采样点','fontname','宋体','fontsize',10);
ylabel('幅值','fontname','宋体','fontsize',10);
wavwrite(y,8000,16,'c101');
wavwrite(s,8000,16,'c102');
end;
loadleleccum
%将信号装入MATLAB的工作环境
s=leleccum
subplot(221);plot(s);%画出原始信号波形
title('原始信号');
%下面利用默认阈值进行消噪处理
[thr,sorh,keepapp]=ddencmp('den','wv',s);%获得信号的默认阈值
s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);
subplot(223);plot(s2);
title('默认阈值消噪后的信号');grid;
%下面利用给定的软阈值进行消噪
cd1soft=wthresh(cd1,'s',t);
cd2soft=wthresh(cd2,'s',t);
cd3soft=wthresh(cd3,'s',t);
c2=[ca3,cd3soft,cd2soft,cd1soft];
s3=waverec(c2,l,'db1');%[C2,L]为新的分解结构
subplot(224);plot(s3);
title('给定软阈值消噪后的信号');grid
硬阈值:Y=wthresh(X,'h',T)
不知对否,大家给意见
给你一个小波应用软阈值和硬阈值处理,然后计算信噪比snr和均方差的mse的程序!
你可以参考一下!~
%--------------------------小波分析心电信号---------------------------------
E=M(:,2);
E=E';
n=size(E);
s=E(1:2000);
%使用db小波对s进行3层分解
for jj=1:8
%jj表示小波的类型
switch jj
case 1
db='db1';
case 2
db='db2';
case 3
db='db3';
case 4
db='db4';
case 5
db='db5';
case 6
db='db6';
case 7
db='db7';
case 8
db='db8';
otherwise
end;
[C L]=wavedec(E,3,db);
% 从c中提取尺度3下的近似小波系数
cA3=detcoef(C,L,3);
%从信号c中提取尺度3,2,1下的细节小波系数
[cD1,cD2,cD3,cD4,cD5]=detcoef(C,L,[1,2,3,4,5]);
%重构尺度5下的近似小波系数
A5=wrcoef('a',C,L,db,5);
% 重构尺度3,2,1下的细节小波系数.
D1=wrcoef('d',C,L,db,1);
D2=wrcoef('d',C,L,db,2);
D3=wrcoef('d',C,L,db,3);
D4=wrcoe
f('d',C,L,db,4);
D5=wrcoef('d',C,L,db,5);
%-------------------使用阈值法去噪-----------------------------
%利用'ddencmp'得到除噪的默认参数
%thr,sorh,keepapp分别表示是阈值,软阈值(或硬阈值)和允许用户保存的低频系数.
%'den','wv',s分别为降噪('cmp'表示压缩),小波('wp'小波包),和降噪信号.
[thr,sorh,keepapp]=ddencmp('den','wv',E);
%去噪
clean=wdencmp('gbl',C,L,db,5,thr,sorh,keepapp);
%---------去噪效果衡量(SNR越大效果越好,MSE越小越好)--------------------------------
%选取信号的长度。
N=n(2);
x=E;
y=clean;
F=0;
M=0;
for ii=1:N
m(ii)=(x(ii)-y(ii))*(x(ii)-y(ii));
t(ii)=y(ii)*y(ii);
f(ii)=t(ii)/m(ii);
F=F+f(ii);
M=M+m(ii);
end;
SNR(jj)=10*log10(F);
MSE(jj)=M/N;
end;
%对比原始信号和除噪后的信号
subplot(2,1,1);
plot(s(1:1000));
title('原始信号');
subplot(2,1,2);
plot(clean(1:1000));
title('除噪后的信号');
SNR,MSE
求助一个小波变换阈值设定的程序
sound=wavread('000');sound=sound';
subplot(4,1,1);plot(sound);title('原始信号');grid;set(gca,'ytick',-1:1:1)
count1=length(sound);
noise=0.07*randn(1,count1);%count1=28000
for i=1:count1
signal(i)=sound(i)+noise(i);
end;
for i=1:count1
y(i)=signal(i);
end;
subplot(4,1,2); plot(y);
title('含噪信号');grid;set(gca,'ytick',-1:1:1)
%读出带噪语音信号,存为‘101.wav'
wavwrite(y,8000,16,'c101');
SNR1 = 10*log10(sum(sound.^2)/sum((y-sound).^2)); %计算输入信噪比
title(['Noise Added (SNR=',num2str(SNR1),'dB)']); % num2str把数字转换成字符串
[C,L]=wavedec(y,3,'db3'); %将信号y分解为3级
ca=appcoef(C,L,'db3',3); %提取一维近似系数,3表示级数,db3为字符串小波的名字
[cd1,cd2,cd3]=detcoef(C,L,[1,2,3]); %3层小波分解
one=mad(cd1); % mad计算样本数据平均绝对偏差
two=mad(cd2); %30---38为阈值函数设置
three=mad(cd3);
n1=one/0.6745;
n2=two/0.6745;
n3=three/0.6745;
men1=n1*sqrt(2*log(count1))/5; %固定阈值形式,阈值为与√2*信号的长度,其中count1为信号长度
men2=n2*sqrt(2*log(count1))/5;
men3=n3*sqrt(2*log(count1))/5;
E1=sum((abs(cd1)).^2)/length(cd1);
E2=sum((abs(cd2)).^2)/length(cd2);
E3=sum((abs(cd3)).^2)/length(cd3);
E=[1,E3/E1,E2/E1];
Em=max(E);
z=E3/E1;
if ((Em==1)&&(E3/E1<0.8))
f=1;
cd1soft=wthresh(cd1,'s',men1); %其中S表示软阈值,如果是H则表示硬阈值
cd2soft=wthresh(cd2,'s',men2);
cd3soft=wthresh(cd3,'s',men3);
c1=[ca,cd3soft,cd2soft,cd1soft];
else
f=0;
casoft=wthresh(ca,'s',0.0050);
cd1soft=wthresh(cd1,'s',men1);
cd2soft=wthresh(cd2,'s',men2);
cd3soft=wthresh(cd3,'s',men3);
c1=[casoft,cd3soft,cd2soft,cd1soft];
end;
t=f;
s3=waverec(c1,L,'db3');
s=s3;
subplot(4,1,3);plot(s);title('消噪信号');grid;set(gca,'ytick',-1:1:1)
SNR2 = 10*log10(sum(sound.^2)/sum((s-sound).^2));
wavwrit
e(s,8000,16,'c102');
此程序是语音信号采用小波变换,改进阈值,的程序,LZ可参考改写
对一个给定信号进行小波变换,就是将该信号按某一小波函数簇展开,即将信号表示为一系列不同尺度和不同时移的小波函数的线性组合,其中每一项的系数称为小波系数,而同一尺度下所有不同时移的小波函数的线性组合称为信号在该尺度下的小波分量 。这是百度百科关于小波系数的原话。可以看出小波系数包含信号在某个尺度下(对应小波分解层数)的高频分量和低频分量。小波变换也就是把信号在小波域里分解成高频分量和低频分量。处理一般都是去除高频分量保留低频分量。