小波降噪-测试


function [Y S]=xiaobojiangzao(input,low_f,high_f,delt,w_name)
%input_x:输入信号
% low_f:带通滤波下截止频率
% high_f:带通滤波上截止频率
% delt:采样率的倒数
% w_name:小波族名称


ls=length(input);
S=zeros(ls, 16);
Y=zeros(ls, 16);


for i=1:16
x=input(:,i); %提取要处理的一路时域信号
y=fft_filter(x,low_f,high_f,delt); %带通滤波


%对肌电信号进行小波分解
[c,l]=wavedec(y,5,w_name);%5层小波分解(层:尺度)
ca5=appcoef(c,l,w_name,5);% 提取一维小波变换低频系数(第五层的逼近系数)
cd5=detcoef(c,l,5);%提取一维小波变换高频系数 (细节系数)
cd4=detcoef(c,l,4);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);


figure;
subplot(611);plot(ca5);title(‘低频系数’);
subplot(612);plot(cd5);title(‘第五层高频系数’);
subplot(613);plot(cd4);title(‘第四层高频系数’);
subplot(614);plot(cd3);title(‘第三层高频系数’);
subplot(615);plot(cd2);title(‘第二层高频系数’);
subplot(616);xlabel(‘小波的各层频率系数’);plot(cd1);title(‘第一层高频系数’);


%小波去噪关键步骤(阈值选取,可改进之处)
%消除噪声处理
[thr,sorh,keepapp]=ddencmp(‘den’,’wv’,y); %函数自动生成小波消噪或压缩的阈值选取方案
s1=wdencmp(‘gbl’,c,l,w_name,5,thr,sorh,keepapp);
S(:,i)=s1;


%’rigrsure’为无偏似然估计阈值类型
thr1=thselect(cd1,’rigrsure’);
thr2=thselect(cd2,’rigrsure’);
thr3=thselect(cd3,’rigrsure’);
thr4=thselect(cd4,’rigrsure’);
thr5=thselect(cd5,’rigrsure’);
%各层的阈值
TR=[thr1,thr2,thr3,thr4,thr5];
%’s’为软阈值;’h’硬阈值。
SORH=’s’;
[s1,CXC,LXC,PERF0,PERF2]=wdencmp(‘lvd’,y,w_name,5,TR,SORH);
% XC为去噪后信号
% [CXC,LXC]为的小波分解结构
% PERF0和PERF2是恢复和压缩的范数百分比。
% ‘lvd’为允许设置各层的阈值,
% ‘gbl’为固定阈值。
% 5为阈值的长度


 


figure;
subplot(411);plot(x);
title(‘初始EMG信号’);
subplot(412);plot(y);
title(‘带通滤波后的信号’);
subplot(413);plot(s1);
title(‘小波去噪后的信号’);


%去除噪声,重构信号,对一维小波系数进行单支重构(重构之后得到的是信号,单支重构的目的在于获得原信号中某频率段的分量信号。
%分量信号的长度(点数)、采样频率均与原信号相同。将单支重构得到的各个分量信号直接按对应点求和,可恢复原信号。)


%重构第5层逼近信号(重构各分解尺度(层)上的信号)
a5=wrcoef(‘a’,c,l,w_name,5);%是对尺度5上的低频部分进行重构,a是低频,d是高频,c,l是wavedec所得
%重构第1~5层细节信号
d5=wrcoef(‘d’,c,l,w_name,5);
d4=wrcoef(‘d’,c,l,w_name,4);
d3=wrcoef(‘d’,c,l,w_name,3);
d2=wrcoef(‘d’,c,l,w_name,2);
d1=wrcoef(‘d’,c,l,w_name,1);


% subplot(414);
y=a5+d4;
Y(:,i)=y;
% plot(y);
% title(‘小波分解重构去噪后的信号’);


end


 


%计算信噪比
F=0;
M=0;
T=0;
% for i=1:ls
% m(i)=(x(i)-y(i))^2;
% t(i)=y(i)^2;
% f(i)=t(i)/m(i);
% F=F+f(i);
% M=M+m(i);
% end;


for i=1:ls
M=M+x(i)^2;
T=T+y(i)^2;
end;
F=M/T;


SNR=10*log10(F)