%%%%%%%%%%%%%%%%%%%%%%%%%空域相關濾波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%參考文獻：
%Wavelet transform domain filters:a spatially selective noise filtration technique
%不能選正交小波，文獻中提出用二次樣條小波（quadratic-spline）.此處用雙正交小波：bior 1.5
%優點：1.對于信噪比高的信號濾波效果好；
%      2.對于邊沿的保護強過閾值濾波，不會產生閾值濾波情況下的過于平滑與Gibbs現象。
%缺點：1.由于對邊沿信號沒做任何處理，所以邊沿可能會有脈沖噪聲保留下來；
%      2.計算相關系數中，如果計算出來的小波系數點位置偏差大，則相關系數計算受影響；
%      3.需要迭代運算，迭代的噪聲能量閾值選取很重要，這里以開始段無信號處估計噪聲；
%      4.需要迭代運算，所以運算量比閾值法大；
%      5.受分解層次影響，在大尺度上小波系數點位置偏差更大，相關系數計算不準確。
%需要具體調整的地方：1.分解的尺度；
%                   2.選定用什么信號作為噪聲的估計；
%                   3.設定停止迭代的噪聲能量閾值參數cc。
%%%%%%%%%%%%%%%%%%%%%%%%%空域相關濾波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clc;
snr=5;
init=2055615866;
x=noisdopp
refx=x;
signal=x;
points=1024;  level=5;  wf='bior 1.5'; %sym8,bior 1.5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %閾值消噪：
S_thr=wden(signal,'rigrsure','s','mln',5,'sym8'); 
% % rigrsure;heursure;sqtwolog; minimaxi
% % one;sln;mln
% subplot(211);plot(signal);
% title('閾值濾波');
% subplot(212);plot(S_thr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%進行二進制小波變換（離散平穩小波變換），并給出各級波形：
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);
[swa,swd] = swt(signal,level,Lo_D,Hi_D);
% figure;
% subplot(level,1,1); plot(real(signal)); grid on;axis tight;
% for i=1:level
%     subplot(level+1,2,2*(i)+1);
%     plot(swa(i,:)); axis tight;grid on;xlabel('time');
%     ylabel(strcat('a   ',num2str(i)));
%     subplot(level+1,2,2*(i)+2);
%     plot(swd(i,:)); axis tight;grid on;
% ylabel(strcat('d   ',num2str(i)));
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%小波系數的處理：
Swd_n=swd; 
swd_org=swd;
mask_n=zeros(size(Swd_n)); %先把系數處理矩陣設置為全0。

for j=1:(level-1) 
    %在1:(level-1)分解層次上對高頻系數處理，最后一層無法求相關系數，所以不作處理。
    Noise_d1=swd_org(j,:);
    Noise_d1=Noise_d1(1:80);
    Noise_var=var(Noise_d1); %以信號的前80個只含有噪聲的點估計噪聲在各層的方差。
    Pw_var=var(swd_org(j,:));
    Corr=swd_org(j,:).*swd_org(j+1,:); %定義相關系數為相鄰兩層的乘積。
      
    cc=1.7; %_______用以設定停止迭代的噪聲能量閾值，需要根據情況調節。________%
    while Pw_var>cc*Noise_var
    Pw=sum(abs(swd(j,:)).^2); %計算小波能量
    Pcorr=sum(abs(Corr).^2); %計算相關系數能量
    Corr_new=Corr.*((Pw/Pcorr)^0.5); %歸一化
     
    corr_mod=abs(Corr_new);
    w_mod=abs(swd(j,:));
    swd_n=swd(j,:).*(corr_mod>w_mod);
    swd_n1=(swd_n~=0);
    mask_n(j,:)=mask_n(j,:)+swd_n1; %將選出的點賦給系數處理矩陣相應位置。
    swd_n0=ones(size(swd_n1));
    swd_n0=swd_n0-swd_n1;
    
    swd(j,:)=swd(j,:).*swd_n0; %將高頻系數選出大值后的地方置0。
    Pw_var=var(swd(j,:));
    Corr_new=Corr_new.*swd_n0; %將相關系數選出大值后的地方置0。
    Corr=Corr_new;
    end
end  
    
mask_max=ones(1,length(mask_n));
mask_n=[mask_n((1:(level-1)),:);mask_max]; %最后一層系數處理矩陣全置1。
Swd_reg=swd_org.*mask_n;

signal_n=iswt(swa,Swd_reg,wf);
% S_mix=wden(signal_n,'sqtwolog','s','sln',5,'sym8'); %rigrsure;heursure;sqtwolog;minimaxi
xcrr=signal_n-xref; %求濾波誤差信號。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%畫圖：
figure; %空域法處理后的高頻系數。
subplot(level+1,1,1); plot(real(signal)); grid on;axis tight;
title('空域法處理后的高頻系數');
for i=1:level
    subplot(level+1,1,i+1);
    plot(Swd_reg(i,:)); axis tight;grid on;
ylabel(strcat('j=   ',num2str(i)));
end

figure; %高頻系數處理前后的比較。
for i=1:level
    subplot(level,2,2*(i)-1);
    plot(swd_org(i,:)); axis tight;grid on;
    ylabel(strcat('d   ',num2str(i)));
    subplot(level,2,2*(i));
    plot(Swd_reg(i,:)); axis tight;grid on;
    ylabel(strcat('j=   ',num2str(i)));
end

figure; %信號濾波前后比較。
subplot(3,1,1);
plot(signal); axis tight;grid on;
axis([0 1000 -5 22]);
title('原始信號');
subplot(3,1,2);
plot(signal_n); axis tight;grid on;
axis([0 1000 -5 22]);
title('空域法濾波后信號');
subplot(3,1,3);
plot(xcrr); axis tight;grid on;
axis([0 1000 -5 22]);
title('濾波誤差信號');

figure; %空域法濾波與閾值濾波的比較。
subplot(3,1,1);
plot(signal); axis tight;grid on;
title('原始信號');
subplot(312);
plot(S_thr);axis tight;grid on;
title('閾值濾波');
subplot(3,1,3);
plot(signal_n); axis tight;grid on;
title('空域法濾波');