clear;
clc;
close all;
fr=1e3;%脈沖重復頻率
tr=1/fr;%脈沖重復周期
fs=1e6;%采樣頻率
ts=1/fs;%采樣周期
%%%%%%%%%%%%%%%%%%%%%%%%%生成測試模塊%%%%%%%%%%%%%%%%%%%%%%
Point=4000;
doppler_vector=2*pi*linspace(0,fr,Point);
t_doppler=(1:16)*tr;
test_matrix=t_doppler'*doppler_vector;
test_matrix_i=cos(test_matrix);
test_matrix_q=sin(test_matrix);
test_matrix=test_matrix_i+test_matrix_q*j;


%%%%%%%%%%%%%%%%%%%%  產生16個FIR濾波器的濾波器組  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a0=[ 5.93141234807e-008,4.347755861717e-007,1.303747145161e-006,2.772196904973e-006,...
  4.74857800029e-006,6.910166596115e-006,8.774077218793e-006,9.856915099804e-006,...
  9.856915099804e-006,8.774077218793e-006,6.910166596115e-006,4.74857800029e-006,...
  2.772196904973e-006,1.303747145161e-006,4.347755861717e-007,5.93141234807e-008]*1e8;
n=1:16;
f=linspace(0.5,1.5,16)*pi;
fn=f'*n;
fnexp=exp(j*fn);
a016=repmat(a0,1,16);
a016=reshape(a016,16,16)';
b=(a016.*fnexp);
%%%%%%%% FIR濾波器的濾波器組的頻率響應曲線  %%%%%%%%%%%%
figure('Name','FIR濾波器的濾波器組的頻率響應曲線');
for n=1:16
    freqz(b(n,:),1,1024,'whole',fr)
    hold on;
end
axis([0 fr -50 100]);
%%%%%%%% FIR濾波器的濾波器組的頻率響應曲線  %%%%%%%%%%%%

yt=test_matrix;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  采用FIR進行MTD處理  %%%%%%%%%%%%%%%%%%%
M=16;
L=length(yt);
for i=1:L
    s_temp=yt(:,i); 
    s_mtd(:,i)=b*s_temp;    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  采用FIR進行MTD處理%%%%%%%%%%%%%%%%%%%%%

mtd_abs=abs(s_mtd);
mtd_abs=mtd_abs/max(max(mtd_abs));
mtd_abs=20*log10(mtd_abs+eps);
t=linspace(0,fr,Point);


figure('Name','實測FIR濾波器組的頻率響應曲線 1');

plot(t,abs(s_mtd(1,:)),'-r'),hold on
plot(t,abs(s_mtd(2,:)),'-g'),hold on
plot(t,abs(s_mtd(3,:)),'-b'),hold on
plot(t,abs(s_mtd(4,:)),'-c'),hold on
plot(t,abs(s_mtd(5,:)),'-m'),hold on;
plot(t,abs(s_mtd(6,:)),'-y'),hold on;
plot(t,abs(s_mtd(7,:)),'-k'),hold on
plot(t,abs(s_mtd(8,:)),'-c'),hold on
plot(t,abs(s_mtd(9,:)),'-r'),hold on
plot(t,abs(s_mtd(10,:)),'-g'),hold on
plot(t,abs(s_mtd(11,:)),'-b'),hold on
plot(t,abs(s_mtd(12,:)),'-c'),hold on
plot(t,abs(s_mtd(13,:)),'-m'),hold on
plot(t,abs(s_mtd(14,:)),'-y'),hold on
plot(t,abs(s_mtd(15,:)),'-k'),hold on
plot(t,abs(s_mtd(16,:)),'-c');grid on;


figure('Name','實測FIR濾波器組的頻率響應曲線 2');

plot(t,mtd_abs(1,:),'-r'),hold on
plot(t,mtd_abs(2,:),'-g'),hold on
plot(t,mtd_abs(3,:),'-b'),hold on
plot(t,mtd_abs(4,:),'-c'),hold on
plot(t,mtd_abs(5,:),'-m'),hold on;
plot(t,mtd_abs(6,:),'-y'),hold on;
plot(t,mtd_abs(7,:),'-k'),hold on
plot(t,mtd_abs(8,:),'-c'),hold on
plot(t,mtd_abs(9,:),'-r'),hold on
plot(t,mtd_abs(10,:),'-g'),hold on
plot(t,mtd_abs(11,:),'-b'),hold on
plot(t,mtd_abs(12,:),'-c'),hold on
plot(t,mtd_abs(13,:),'-m'),hold on
plot(t,mtd_abs(14,:),'-y'),hold on
plot(t,mtd_abs(15,:),'-k'),hold on
plot(t,mtd_abs(16,:),'-c');grid on;
axis([0,fr,-70,0])

figure('Name','實測FIR濾波器組的頻率響應曲線包絡');
plot(t,20*log10(max(abs(s_mtd))/max(max(abs(s_mtd)))),'-r');grid on;