比较了软阈值,硬阈值及当今各种阈值计算方法和阈值函数处理方法的性能

小波阈值去噪，比较了软阈值，硬阈值及当今各种阈值计算方法和阈值函数处理方法的性能，通过信噪比及均方差的比较，得出各种的算法的优劣
clear;clc;
%-----------------------------------------------------------------
%从程序的运行结果来看，文献1，3，4的去噪效果比较好
%其中文献4对高信噪比的的情况不是很好，在高信噪比时，软硬阈值的效果最好
%-----------------------------------------------------------------------
%测试数据的选取
% fun='blocks';
% snr=5;
% jN=5; %分解的层数
% N=13; %数据长度为2^N
% [x,s]=wnoise(4,N,sqrt(snr));
% ps=sum(x.^2)/length(x);
% sigma_noise=abs(sqrt(ps/(10^(snr/10))));
% noise=sigma_noise*randn(1,length(x));%noise噪声的方差是(sigma_noise.^2)
% s=x+noise;
% figure,
% subplot(211);plot(x);
% subplot(212);plot(s);
% subplot(211);plot(x);title('纯净信号x');
% subplot(212);plot(s);title('混合信号x');

%调幅信号,纯净信号
tic;
% fs=5e+6; %采样率50M
% ts=1/fs;
% fc=10.7e+6; %载频10.7M
% t0=2; %数据长度N=t0*fs;
% t=[0:ts:t0]; %模拟信号的数字化
% m=sinc(100*t); %消息信号
% c=cos(2*pi*fs.*t);%载波信号
% x=m.*c;
% N=t0*fs;
% toc;
%t0=.2; % signal duration
%ts=0.001; % sampling interval
%fc=250; % carrier frequency
%snr=20; % SNR in dB (logarithmic)
%fs=1/ts; % sampling frequency
%t=[-t0/2:ts:t0/2]; % time vector
%m=sinc(100*t); % the message signal
%c=cos(2*pi*fc.*t); % the carrier signal
%x=m.*c; % the DSB-AM modulated signal
%N=t0*fs;
%---------------------------------------------------------------------
%-------------------------BPSK信号,数字调制-----------------------------------
codes=6; %码元个数，即输入调制信号的长度
sigma=1; %调制信号的幅度
fs=600e3; %采样率600KHz
fb=1e3; %波特率1KHz，fb表示对输入调制信号的采样率
fc=100e3; %载频100KHz
Modulate=2; %为选择调制方式
N0=fs/fb; %一个码元周期内的采样点数，一个输入信号长度内的采样点数
N=N0*codes; %总的采样点数（已调信号的长度）
p0=pi*rand(1,1); %产生初始相位
symbols=randint(1,codes,[0,1]); %产生基带码元
x_B = ones(N0,1)*symbols;
x_BB = x_B(:)'; %根据波特率要求产生码元
signal_base = x_BB; %产生基带信号
signal=sigma*dmod(symbols,fc,fb,[fs p0],'psk',Modulate);
%产生psk调制信号，p0是载频的初始相位
x=si

gnal;
%-------------------------加入指定强度的噪声---------------------------------
snr=5;
ps=sum(x.^2)/N;
sigma_noise=abs(sqrt(ps/(10^(snr/10))));
nn=randn(1,N);
enn=sum(nn)/N; %随机数nn的均值
nn=nn-enn; %使nn均值为0
noise=sigma_noise*nn;
s=x+noise;
wname='db7';
jN=6; %分解的层数
[c,l]=wavedec(s,jN,wname);
snrs=20*log10(norm(x)/norm(s-x));
mmses=mmse(s-x);

%高频分量的索引
first = cumsum(l)+1;
first1=first;
first = first(end-2:-1:1);
ld = l(end-1:-1:2);
last = first+ld-1;
%--------------------------------------------------------------------------
%软阈值与硬阈值
cxdsoft=c;
cxdhard=c;
for j=1:jN %j是分解尺度
flk = first(j):last(j);%flk是di在c中的索引
thr(j)=sqrt(log(length(flk)))/log(j+1);
for k=0:(length(flk)-1)%k是位移尺度
djk=c(first(j)+k);%为了简化程序
djk2=c(first(j)+k);
absdjk=abs(djk);
thr1=thr(j);
%阈值处理
if absdjkthr1
djk=djk-thr1/exp((djk-thr1)/n);
else
djk=djk+thr1/exp((-djk-thr1)/n);
end
cxd51(first(j)+k)=djk;
end
end
%新方法二
cxd52=c;
for j=1:jN
flk = first(j):last(j);
for k=0:(length(flk)-1)
djk=c(first(j)+k);%为了简化程序
absdjk=abs(djk);
thr1=thr(j);
if absdjkthr1
djk=djk-thr1+thr1/(2*b+1);
else
djk=djk+thr1-thr1/(2*b+1);
end
cxd53(first(j)+k)=djk;
end
end
%新方法四
cxd54=c;
for j=1:jN
flk = first(j):last(j);
for k=0:(length(flk)-1)
djk=c(first(j)+k);%为了简化程序
absdjk=abs(djk);
thr1=thr(j);
bb=(2*aa+1+m*sqrt((2*aa+1).^2-pi.^2))/(pi*thr1.^(2*n+1));
kk=(1+(bb*thr1.^(2*n+1)).^2)/((2*a+1)*bb*thr1.^(2*n));
if absdjkthr1
djk=djk-thr1+kk*atan(bb*djk.^(2*aa+1));
else
djk=djk+thr1-kk*atan(bb*djk.^(2*aa+1));
end
cxd54(first(j)+k)=djk;
end
end
%信号重构
snew51=waverec(cxd51,l,wname);
snew52=waverec(cxd52,l,wname);
snew53=waverec(cxd53,l,wname);
snew54=waverec(cxd54,l,wname);
%信噪比及均方差分析
snry51=20*log10(norm(x)/norm(snew51-x));
snry52=20*log10(norm(x)/norm(snew52-x));
snry53=20*log10(norm(x)/norm(snew53-x));
snry54=20*log10(norm(x)/norm(snew54-x));
mmsey51=mmse(snew51-x);
mmsey52=mmse(snew52-x);
mmsey53=mmse(snew53-x);
mmsey54=mmse(snew54-x);
%--------------------------------------------------------------------------
%阈值处方法六
%文献【一种改进的小波域阈值去噪算法】作者付炜
m=40;
cxd6=c;
for j=1:jN
flk = first(j):last(j);
thr1=thselect(s,'heursure');
for k=0:(length(flk)-1

)
djk=c(first(j)+k);%为了简化程序
absdjk=abs(djk);
thr2=thr1/3;
if absdjk>thr1
djk=djk'*(absdjk.^m-thr1.^m)/absdjk.^m;
elseif djk