小波降噪-測試


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)