%MTD濾波器組設計
%
%tic;
clear all;
clc;
N=32;%階數可變
Tr=*****;Fr=1/Tr;%可修改的脈沖重復周期，%數值由自己填寫
f=0:2:Fr;  
delta2_clutter=16;  %雜波方差確定零陷凹口的大小
sigma2_clutter=10^6;  %雜波功率確定凹口深度
sigma2_noise=1;    % 噪聲功率
%for fi=200:65:4035
Df=40;   %主瓣頻率和副瓣衰減
interval=20;  %interval_former=20
f_interrupt=2:interval:Fr; %干擾頻率    
L=length(f_interrupt);  %干擾頻率個數
%%%%%%%%%%%%%%%%%%%%%%基本參數設定%%%%%%%%%%%%%%%%%%%%%%
a_fi=a_f(N,fi,Tr); %導頻矢量函數
%
for j=1:N          %雜波協方差矩陣Rc
    for i=j:N
     Rc(i,j)=sigma2_clutter*exp((-2)*(pi^2)...
         *delta2_clutter*(j-i)^2*Tr^2);
     Rc(j,i)=Rc(i,j);
    end
end
%%
for i=1:length(f)
    a_f_fd(:,:,i)=a_f(N,f(i),Tr);    %導頻矢量組成的矩陣并且降維過后的函數
end
[a_f_fd_row,a_f_fd_col,z]=size(a_f_fd);
for i=1:a_f_fd_row
    for j=1:z
        a_f_fd2(i,j)=a_f_fd(i,1,j);
    end
end
%初始化干擾功率
for i=1:L
    sigma2_interrupt(1,i)=0;
end
%%%%%%%%%%%%%%%%%以上不需改動%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=0;
%%%%%%%%%%%%%%%%%修改參數部分%%%%%%%%%%%%%%%%%%%%%%%%%%%
K=0.08;                               %迭代增益0.1
K1=4;                                 %K1=3
K2=8;                                 %K2=5
K3=1.05;
K4=1.03;
for k=0:20                            %迭代次數20
Rc_inter=zeros(size(Rc));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%加權和
for j=1:L
a_f2=a_f(N,f_interrupt(j),Tr)*ctranspose(a_f(N,f_interrupt(j),Tr));
Rc_interrupt(:,:,j)=sigma2_interrupt(k+1,j).*a_f2;
Rc_inter=Rc_inter+Rc_interrupt(:,:,j);
end
%構造權矢量函數
Ra(:,:,k+1)=Rc+sigma2_noise.*eye(size(Rc))+Rc_inter;   %協方差矩陣
Wi(:,:,k+1)=inv(Ra(:,:,k+1))*a_fi;%權矢量函數

%%%%%%%%%%頻響應函數%%%%%%%%%%%%%%%%%
Hi(k+1,:)=ctranspose(Wi(:,:,k+1))*a_f_fd2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h(k+1)=abs(max(Hi(k+1,:)));
for i=1:L
    d(i,k+1)=h(k+1)/(10^(Df/20));
end
for i=1:L
    eta(i)=sigma2_interrupt(k+1,i)+K*(abs(Hi(k+1,((i-1)*interval+2)/2+1))...
                                       -d(i,k+1));
    sigma2_interrupt(k+2,i)=max(0,eta(i));
    if(i<=1|i>=L-1)
        sigma2_interrupt(k+2,i)=0;
    end
  if((i>(fi-K1*interval)/interval)&(i<(fi+K2*interval)/interval))  
      sigma2_interrupt(k+2,i)=0;
  end
 if(i>=L-2)
       sigma2_interrupt(k+2,i)=K3*sigma2_interrupt(k+2,i);
   end
  if(i>=8&i<=10)
     sigma2_interrupt(k+2,i)=K4*sigma2_interrupt(k+2,i);  
   end
    %%%%%%%%%%
end
end
W(:,:,1)=Wi(:,:,k+1);
plot(f,20*log10(abs(Hi(k+1,:)/max(abs(Hi(k+1,:))))));
% hold on;
% plot(f,abs(Hi(k+1,:)));
% title('修改后第一個濾波器的增益函數');xlabel('頻率/Hz');ylabel('增益值');