硬阈值和软阈值区别分析及选取技巧

　　阀值化处理就是当灰度值低于一个值（阀值，根据你的需要设定）时，让其成为一个小的值;当灰度值大于一个值（阀值）时，让其成为一个大的值。

　　比如：灰度在0-255， 当灰度小于128时赋值为0，大于128时赋值为255，即亮的地方更亮，暗的地方更暗。即实现了阀值分割。这样图像就黑白分明，对比度加大了。

　　硬阀值：是被检测数据所不能逾越的阀值。

　　软阀值：则规定被检测数据的变动范围。

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

　　小波阈值去噪，比较了软阈值，硬阈值及当今各种阈值计算方法和阈值函数处理方法的性能，通过信噪比及均方差的比较，得出各种的算法的优劣。

　　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