gd=0;
end
function[snr,mse]=snrmse（I,In）
%计算信噪比函数
%I:
原始信号
%In:
去噪后信号
snr=0;
Ps=sum（sum（（I-mean（mean（I）））.^2））;%signalpower
Pn=sum（sum（（I-In）.^2））; %noisepower
snr=10*log10（Ps/Pn）;
mse=Pn/length（I）;
QQ截图20130516175535.png （11.18KB,下载次数:
0）
改进函数表达式
本帖最后由罗志雄于2013-5-1621:
58编辑
function[snr,mse]=snrmse（I,In）
%计算信噪比函数
%I:
原始信号
%In:
去噪后信号
snr=0;
Ps=sum（sum（（I-mean（mean（I）））.^2））;%signalpower
Pn=sum（sum（（I-In）.^2））; %noisepower
snr=10*log10（Ps/Pn）;
mse=Pn/length（I）;
修改后程序清单如下：
clearall;clc;closeall;
fs=8000; %语音信号采样频率为8000
xx=wavread（'lw1.wav'）;
x1=xx（:
1）;%取单声道
x1=x1-mean（x1）;
t=（0:
length（x1）-1）/8000;
y1=fft（x1,2048）; %对信号做2048点FFT变换
f=fs*（0:
1023）/2048;
figure
（1）
plot（t,x1） %做原始语音信号的时域图形
y=awgn（x1',10,'measured'）; %加10db的高斯白噪声
[snr,mse]=snrmsel（x1',y） %求得信噪比均方误差
snr1=SNR_singlech（x1',y）
figure
（2）
plot（t,y） %做加噪语音信号的时域图形
[c,l]=wavedec（y,3,'db1'）;%多尺度一维分解
%用db1小波对信号进行3层分解并提取系数
a3=appcoef（c,l,'db1',3）;
%a2=appcoef（c,l,'db1',2）;
%a1=appcoef（c,l,'db1',1）;
d3=detcoef（c,l,3）;
d2=detcoef（c,l,2）;
d1=detcoef（c,l,1）;
thr1=thselect（d1,'rigrsure'）;%阈值获取，使用Stein的无偏风险估计原理
thr2=thselect（d2,'rigrsure'）;
thr3=thselect（d3,'rigrsure'）;
%利用改进阈值函数进行去噪处理
gd1=Garrote_gg（d1,thr1）;
gd2=Garrote_gg（d2,thr2）;
gd3=Garrote_gg（d3,thr3）;
c1=[a3gd3gd2gd1];
functiongd=Garrote_gg（a,b）%a为信号分解后的小波系数，b为获得的阈值
m=0.2*（（a.*a）-（b*b））;
if（abs（a）>=b）
gd=sign（a）*（abs（a）-b/exp（m））;
else
gd=zeros（size（a））;
end
y1=waverec（c1,l,'db1'）;%多尺度重构
[snr,mse]=snrmsel（x1',y1）%求得信噪比均方误差
figure（3）;
plot（t,y1）;
小波去噪软阈值和硬阈值的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=[a3ythard3ythard2ythard1];
s3=waverec（c2,l,'db1'）;
%进行软阈值处理
ytsoftd1=wthresh（d1,'s',thr）;
ytsoftd2=wthresh（d2,'s',thr）;
ytsoftd3=wthresh（d3,'s',thr）;
c3=[a3ytsoftd3ytsoftd2ytsoftd1];
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小波除噪，为何硬阈值和软阈值除躁信噪比一样了？
matlab小波除噪，为何硬阈值和软阈值除躁信噪比一样了？
loadleleccum;
index=1:
1024;
f1=leleccum（index）;%产生含噪信号
init=2055615866;
randn（'seed',init）;
f2=f1+18*randn（size（x））;
snr=SNR_singlech（f1,f2）%信噪比
subplot（2,2,1）;plot（f1）;title（'含噪信号'）;%axis（[1,1024,-1,1]）;
subplot（2,2,2）;plot（f2）;title（'含噪信号'）;%axis（[1,1024,-1,1]）;
%用db5小波对原始信号进行3层分解并提取系数
[c,l]=wavedec（f2,3,'db6'）;
a3=appcoef（c,l,'db6',3）;
d3=detcoef（c,l,3）;
d2=detcoef（c,l,2）;
d1=detcoef（c,l,1）;
sigma=wnoisest（c,l,1）;
thr=wbmpen（c,l,sigma,2）;
%进行硬阈值处理
ythard1=wthresh（d1,'h',thr）;
ythard2=wthresh（d2,'h',thr）;
ythard3=wthresh（d3,'h',thr）;
c2=[a3ythard3ythard2ythard1];
f3=waverec（c2,l,'db6'）;
%进行软阈值处理
ytsoftd1=wthresh（d1,'s',thr）;
ytsoftd2=wthresh（d2,'s',thr）;
ytsoftd3=wthresh（d3,'s',thr）;
c3=[a3ytsoftd3ytsoftd2ytsoftd1];
f4=waverec（c3,l,'db6'）;
%对上述信号进行图示
subplot（2,2,3）;plot（f3）;title（'硬阈值处理'）;%axis（[1,1024,-1,1]）;
subplot（2,2,4）;plot（f4）;title（'软阈值处理'）;%axis（[1,1024,-1,1]）;
snr=SNR_singlech（f1,f3）
snr=SNR_singlech（f1,f4）
信噪比函数SNR_singlech（I,In）
functionsnr=SNR_singlech（I,In）
%计算信噪比函数
%I：
riginalsignal
%In:
noisysignal（ie.originalsignal+noisesignal）
snr=0;
Ps=sum（sum（（I-mean（mean（I）））.^2））;%signalpower
Pn=sum（sum（（I-In）.^2））; %noisepower
snr=10*log10（Ps/Pn）;
小波去噪程序
Matlab小波去噪（默认，强制，给定三种情况）
%%利用小波分析对监测采集的信号进行去噪处理，恢复原始信号
%小波分析进行去噪有3中方法：
%1、默认阈值去噪处理。
该方法利用函数ddencmp（）生成信号的默认阈值，然后利用函数wdencmp（）进行去噪处理；
%2、给定阈值去噪处理。
在实际的去噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。
在进行阈值量化处理时可利用函数wthresh（）;
%3、强制去噪处理。
该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。
这种方法比较简单，且去噪后的信号比较平滑，但是容易丢失信号中的有用成分。
%%利用小波分析对监测采集的水轮机信号进行去噪处理，恢复原始信号
%ProgramStart
%%载入监测所得信号
loaddefault.txt; %装载采集的信号
x=default;
lx=length（x）;
t=[0:
1:
length（x）-1]';
%%绘制监测所得信号
subplot（2,2,1）;
plot（t,x）;
title（'原始信号'）;
gridon
set（gcf,'color','w'）
set（gca,'fontname','timesNewRoman'）
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（'强制去噪后信号'）;
gridon
set（gcf,'color','w'）
set（gca,'fontname','timesNewRoman'）
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（'默认阈值去噪后信号'）;
gridon
set（gcf,'color','w'）
set（gca,'fontname','timesNewRoman'）
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（'给定软阈值去噪后信号'）;
gridon
set（gcf,'color','w'）
set（gca,'fontname','timesNewRoman'）
set（gca,'fontsize',14.0）
以上就是三种小波去噪的原程序。
但红色标注的地方，我却运行不过去。
提示分别问：
1，?
?
?
Errorusing==>horzcat
CATargumentsdimensionsarenotconsistent.
Errorin==>ttat33
c1=[ca3,cdd3,cdd2,cdd1];
2，?
?
?
Errorusing==>wthreshat30
Invalidargumentvalue.
Errorin==>ttat54
cd1soft=wthresh（cd1,'x',1.465）;%经验给出软阈值数
请问为什么？
？
？
关于小波新阈值去噪改进方法里，matlab程序编写的讨论
新阈值能有效地克服软阈值去噪方法中由于估计值与真实值之间的恒定偏差而带来的去噪误差，也能有效地抑制硬阈值去噪方法中易产生的信号振荡现象。
但不懂它程序的编写。
请教高手帮忙解决一下！
谢谢了！
functiony=wthresh（x,sorh,t）
%WTHRESHPerformsoftorhardthresholding.
% Y=WTHRESH（X,SORH,T）returnssoft（ifSORH='s'）
% orhard（ifSORH='h'）T-thresholding oftheinput
% vectorormatrixX.Tisthethresholdvalue.
%
% Y=WTHRESH（X,'s',T）returnsY=SIGN（X）.（|X|-T）+,soft
% thresholdingisshrinkage.
%
% Y=WTHRESH（X,'h',T）returnsY=X.1_（|X|>T）,hard
% thresholdingiscruder.
%
% SeealsoWDEN,WDENCMP,WPDENCMP.
% M.Misiti,Y.Misiti,G.Oppenheim,J.M.Poggi12-Mar-96.
% LastRevision:
24-Jul-2007.
% Copyright1995-2007TheMathWorks,Inc.
% $Revision:
1.11.4.3$
switchsorh
case's'
tmp=（abs（x）-t）;
tmp=（tmp+abs（tmp））/2;
y =sign（x）.*tmp;
case'h'
y =x.*（abs（x）>t）;
otherwise
error（'Wavelet:
FunctionArgVal:
Invalid_ArgVal',...
'Invalidargumentvalue.'）
end
这是传统的阈值去噪法，只要你有新的算法，直接改就行了，但是我对该法有些疑惑，wthresh（x,sorh,thr）,若我用赵瑞珍法算出了每个尺度（假设J=3）的阈值，那么还用wthresh（x,sorh,thr）去噪时，thr是不是应该在之前写成thr=[thr1,thr2,thr3]呢？
求助高手：
Matlab小波去噪中遇见的问题
小波程序如下：
function[y]=f1（x）;
X=xlsread（'shoulian.xls'）;%读入表格数据
x=X（:
1）;%将表格数据的第一列赋值给x轴
y=X（:
2）;%将表格数据的第二列赋值给y轴
subplot（2,2,1）;
plot（x,y,'-m*'）%画出图形
holdon;
title（'原始图形'）;
grid;%画出网格
%用db4小波对原始信号进行3层分解并提取系数 还没有写好:
要增加高频系数图
[c,l]=wavedec（y,3,'db4'）;%对y进行3层分解
ca3=appcoef（c,l,'db4',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'cdd3cdd2cdd1];
s1=waverec（c1,l,'db4'）;
subplot（2,2,2）;
plot（x,s1,'-m*'）;
holdon;
title（'强制去噪后的信号'）;grid;
%用默认阈值对信号进行去噪处理并图示结果
%用ddencmp（）函数获得信号的默认阈值,使用wdencmp（）命令函数实现去噪过程
[thr,sorh,keepapp]=ddencmp（'den','wv',y）;
s2=wdencmp（'gbl',c,l,'db4',3,thr,sorh,keepapp）;
subplot（2,2,3）;
plot（x,s2,'-m*'）;
holdon;
title（'默认阈值去噪后的信号'）;grid;
%用给定的软阈值进行去噪处理
cdd1soft=wthresh（cd1,'s',1.465）;
cdd2soft=wthresh（cd2,'s',1.823）;
cdd3soft=wthresh（cd3,'s',2.768）;
c2=[ca3cdd3softcdd2softcdd1soft];
s3=waverec（c2,l,'db4'）;
subplot（2,2,4）;
plot（x,s3,'-m*'）;
title（'给定软阈值去噪后的信号'）;grid
clear;
错误显示：
?
?
?
Errorusing==>horzcat
Allmatricesonarowinthebracketedexpressionmusthavethe
samenumberofrows.
Errorin==>xnnat41
c2=[ca3cdd3softcdd2softcdd1soft];
i以至于最后一个图形没法显示
同时处理后的图形数据怎么调出来，以及第一个数据应该是【0，0】，但是处理后的不是。
请问是我选择的分解重构函数部队还是其他原因
请问如何修改
需要说明的是：
收敛表格的两组数据为：
0
0.00
2
0.49
4
0.94
5
1.03
6
1.03
7
1.48
10
1.49
11
1.86
14
1.96
19
1.78
26
2.75
28
2.78
31
2.71
34
2.54
39
2.85
47
2.71
49
2.71
54
5.43
61
5.47
68
5.47
75
5.20
请教程序问题，心电信号去噪的matlab程序
想把这个程序。
改为只用软硬阈值对比的心电信号去噪分析，该怎么办？
求高手。
程序见下：
%%%%%%%%%%信号小波分解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%基于Haar小波%%%%%%%
ecg=fopen（'100.dat','r'）;
N=1201;
升级会员 