%212 Format:360Fre -> 540(1 Channal)|1080(2 Channal).
%%%%------------->心電信號的讀取!
%readheader('E:\MATLAB\develop\100.hea')
clear all;
NN={100:350 200:250 100:300 100:150};           %各層選取的只含噪聲點
fid1 = fopen('100.dat','rb');
[A] = fread(fid1,1080*4,'uchar');
j=1;
for i=1:360*4
   %212->12bits.
   B(i)=bitand(A(j+1),15)*2^8+A(j);             %取低4位，左移8位，成12位
   C(i)=bitand(A(j+1),240)*2^4+A(j+2);          %取高4位，左移4位，成12位
   j=j+3;
   %Signed.
   if B(i)>=2048
      B(i)=(4096-B(i))*(-1);
   end   
   
   if C(i)>=2048
      C(i)=(4096-C(i))*(-1);
   end      
end   
figure;
subplot(211);plot(B);
subplot(212);plot(C);

fid2 = fopen('100A.txt','wt');                  %文件寫入
fprintf(fid2,'%d\n',B);
fid3 = fopen('100B.txt','wt');
fprintf(fid3,'%d\n',C);

fclose(fid1);
fclose(fid2);
fclose(fid3);

%%
signal=B;
level=5;wf='bior1.5';                           %濾波選用的小波函數

%%%%%%%%%%%%%%%%%%%%%各種濾波方法
%-------->1.軟閾值濾波
S_thr1=wden(signal,'rigrsure','s','mln',5,'sym8'); %規則:rigrsure,尺度改變比例:mln
figure;
subplot(211);
plot(B);axis tight;grid on;axis([0 1500 800 1300]);
title('原始信號');
subplot(212);
plot(S_thr1);axis tight;grid on;axis([0 1500 800 1300]);
title('軟閾值濾波后波形');

%-------->2.硬閾值濾波
S_thr2=wden(signal,'rigrsure','h','mln',5,'sym8'); %規則:rigrsure,尺度改變比例:mln
figure;
subplot(211);
plot(B);axis tight;grid on;axis([0 1500 800 1300]);
title('原始信號');
subplot(212);
plot(S_thr2);axis tight;grid on;axis([0 1500 800 1300]);
title('硬閾值濾波后波形');

%-------->3.FIR濾波器進行濾波(切比雪夫Ⅱ型 低通濾波器)
in=B;
fs=360;                                         %采樣頻率
fn=fs/2;                                        %奈奎斯特采樣率
Wp=10/fn;Ws=100/fn;                             %通帶截止頻率及阻帶截止頻率
Rp=3;As=60;                                     %通帶最大衰減Rp=3 dB-----阻帶最小衰減As=60
[n,Wn]=cheb2ord(Wp,Ws, Rp,As)   
[b, a]=cheby2(n,Rp,Wn,'low');
%%繪制信號
figure;
subplot(3, 2, [1 2]);
plot(in); 
grid on;
ylabel('幅度');
title('輸入信號');
out1=filter(b, a, in);
subplot(3, 2, [5 6]);
plot(out1); 
grid on;axis([0 1500 800 1300]);
ylabel('幅度');
title('FIR低通濾波后輸出信號');

%-------->4.IIR濾波器進行濾波(巴特沃斯 低通數字濾波器)
in=B;
fp=40;fs=310;Fs=1000;Rp=3;Rs=60;T=1/Fs;         %指標設計
Wlp=2*tan(2*pi*fp*T/2)/pi;Wls=2*tan(2*pi*fs*T/2)/pi;    %求歸一化頻率
[N,Wn]=buttord(Wlp,Wls,Rp,Rs);                  %確定最小階數N和頻率參數Wn
[bp,ap]=butter(N,1,'s');                        %得歸一化低通原型
[bs,as]=lp2lp(bp,ap,Wn*pi*Fs);                  %轉換為模擬低通
[b,a]=bilinear(bs,as,Fs);                       %雙線性法進行模數轉換
out2=filter(b,a,in);                             %濾波
                %%繪制圖形
figure;
subplot(3, 2, [1 2]);
plot(in); 
grid on;
ylabel('幅度');
title('輸入信號');
subplot(3, 2, [5 6]);
plot(out2); 
grid on;axis([0 1500 800 1300]);
ylabel('幅度');
title('IIR低通濾波后輸出信號');

%-------->5.空域相關法濾波

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);             %設計指定濾波器
[swa,swd] = swt(signal,level,Lo_D,Hi_D);        %離散平穩小波變換
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(NN{j});
    Noise_var=var(Noise_d1);                    %以信號的前只含有噪聲的點估計噪聲在各層的方差，個數為NN數組中定義
    Pw_var=var(swd_org(j,:));
    Corr=swd_org(j,:).*swd_org(j+1,:);          %定義相關系數為相鄰兩層的乘積。
      
    cc=1.2;                                     %===》用以設定停止迭代的噪聲能量閾值，需要根據情況調節
    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);           %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);
xcrr=signal_n-B;                                %求濾波誤差信號。


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

figure;                                         %空域法濾波與其他濾波方法的比較。
subplot(6,1,1);
plot(signal); axis tight;grid on;axis([0 1500 900 1230]);
title('原始信號');
subplot(612);
plot(S_thr1);axis tight;grid on;axis([0 1500 900 1230]);
title('軟閾值濾波');
subplot(6,1,3);
plot(S_thr2); axis tight;grid on;axis([0 1500 900 1230]);
title('硬閾值濾波');
subplot(614);
plot(out1);axis tight;grid on;axis([0 1500 900 1230]);
title('FIR濾波');
subplot(615);
plot(out2);axis tight;grid on;axis([0 1500 900 1230]);
title('IIR濾波');
subplot(616);
plot(signal_n);axis tight;grid on;axis([0 1500 900 1230]);
title('空域相關濾波');

xzb1=10*log10(var(S_thr1)/var(signal-S_thr1));
xzb2=10*log10(var(S_thr2)/var(signal-S_thr2));
xzb3=10*log10(var(out1)/var(signal-out1));
xzb4=10*log10(var(out2)/var(signal-out2));
xzb5=10*log10(var(signal_n)/var(signal-signal_n));
fprintf('軟閾值濾波信噪比為：%f',xzb1)
fprintf('硬閾值濾波信噪比為：%f',xzb2)
fprintf('FIR濾波信噪比為：%f',xzb3)
fprintf('IIR濾波信噪比為：%f',xzb4)
fprintf('空域相關濾波信噪比為：%f',xzb5)