%===========================================================================================%
%            該程序完成8個脈沖信號的脈壓、MTI/MTD                                           %
%===========================================================================================%
close all; %關閉所有圖形
clear all; %清除所有變量
%===================================================================================%
%                                    雷達參數(shù)                                       %
%===================================================================================%
C=3.0e8;  %光速(m/s)
RF=3.140e9;  %雷達射頻
Lambda=C/RF;%雷達工作波長
PulseNumber=8;   %回波脈沖數(shù) 
BandWidth=2.0e6;  %發(fā)射信號帶寬
TimeWidth=42.0e-6; %發(fā)射信號時寬
PRT=120e-6;   % 雷達發(fā)射脈沖重復周期(s),120us對應1/2*120*300=18000米
PRF=1/PRT;
Fs=2.0e6;  %采樣頻率
NoisePower=-12;%(dB);%噪聲功率（目標為0dB）
%---------------------------------------------------------------%
SampleNumber=fix(Fs*PRT);%計算一個脈沖周期的采樣點數(shù)；
TotalNumber=SampleNumber*PulseNumber;%總的采樣點數(shù)；
BlindNumber=fix(Fs*TimeWidth);%計算一個脈沖周期的盲區(qū)-遮擋樣點數(shù)；
%===================================================================================%
%                                    目標參數(shù)                                       %
%===================================================================================%
TargetNumber=3;%目標個數(shù)
SigPower(1:TargetNumber)=[1 1 1];%目標功率,無量綱
TargetDistance(1:TargetNumber)=[3000 8025 10125];%目標距離,單位m
DelayNumber(1:TargetNumber)=fix(Fs*2*TargetDistance(1:TargetNumber)/C);% 把目標距離換算成采樣點（距離門）
TargetVelocity(1:TargetNumber)=[0 50 450];%目標徑向速度 單位m/s
TargetFd(1:TargetNumber)=2*TargetVelocity(1:TargetNumber)/Lambda; %計算目標多卜勒
%====================================================================================%
%                                   產生線性調頻信號                                     %
%====================================================================================%
number=fix(Fs*TimeWidth);%回波的采樣點數(shù)=脈壓系數(shù)長度=暫態(tài)點數(shù)目+1
if rem(number,2)~=0
   number=number+1;
end   
for i=-fix(number/2):fix(number/2)-1
   Chirp(i+fix(number/2)+1)=exp(j*(pi*(BandWidth/TimeWidth)*(i/Fs)^2));
end
coeff=conj(fliplr(Chirp));
%-------------------------產生目標回波串------------------------%
SignalAll=zeros(1,TotalNumber);%所有脈沖的信號,先填0
for k=1:TargetNumber% 依次產生各個目標
   SignalTemp=zeros(1,SampleNumber);% 一個脈沖
   SignalTemp(DelayNumber(k)+1:DelayNumber(k)+number)=sqrt(SigPower(k))*Chirp;%一個脈沖的1個目標（未加多普勒速度）
   Signal=zeros(1,TotalNumber);
   for i=1:PulseNumber
      Signal((i-1)*SampleNumber+1:i*SampleNumber)=SignalTemp;
   end
   FreqMove=exp(j*2*pi*TargetFd(k)*(0:TotalNumber-1)/Fs);%目標的多普勒速度*時間=目標的多普勒相移
   Signal=Signal.*FreqMove;
   SignalAll=SignalAll+Signal;
end
figure(1);
subplot(2,1,1);plot(real(SignalAll),'r-');title('目標信號的實部');
subplot(2,1,2);plot(imag(SignalAll));title('目標信號的虛部');
grid on;zoom on;

%====================================================================================%
%                                   產生系統(tǒng)噪聲信號                                  %
%====================================================================================%
SystemNoise=normrnd(0,10^(NoisePower/10),1,TotalNumber)+j*normrnd(0,10^(NoisePower/10),1,TotalNumber);
%====================================================================================%
%                                   總的回波信號                                     %
%====================================================================================%
Echo=SignalAll+SystemNoise;% +SeaClutter+TerraClutter;
for i=1:PulseNumber   %在接收機閉鎖期,接收的回波為0
      Echo((i-1)*SampleNumber+1:(i-1)*SampleNumber+number)=0;
end
figure(2);
subplot(2,1,1);plot(real(Echo));title('總回波信號的實部,閉鎖期為0');
subplot(2,1,2);plot(imag(Echo));title('總回波信號的虛部,閉鎖期為0');
%================================時域脈壓=================================%
pc_time0=conv(Echo,coeff);
figure(3);plot(abs(pc_time0));title('時域脈壓結果的幅度,有暫態(tài)點');
pc_time1=pc_time0(number:TotalNumber+number-1);%去掉暫態(tài)點 number-1個
%================================頻域脈壓=================================%
% Echo_fft=fft(Echo,TotalNumber+number-1);
% coeff_fft=fft(coeff,TotalNumber+number-1);
Echo_fft=fft(Echo,2048);%理應進行TotalNumber+number-1點FFT,但為了提高運算速度,進行了2048點的FFT
coeff_fft=fft(coeff,2048);
pc_fft=Echo_fft.*coeff_fft;
pc_freq0=ifft(pc_fft);
figure(4);subplot(2,1,1);plot(abs(pc_freq0));title('頻域脈壓結果的幅度,有前暫態(tài)點和后暫態(tài)點');
hold on;
subplot(2,1,2);plot(abs(pc_time0(1:TotalNumber+number-1)-pc_freq0(1:TotalNumber+number-1)),'r');title('紅色為時域頻域脈壓的差別');
pc_freq1=pc_freq0(number:TotalNumber+number-1);%去掉暫態(tài)點 number-1個,后填充點若干(2048-number+1-TotalNumber=45個)
%================按照脈沖號、距離門號重排數(shù)據(jù)=================================%
for i=1:PulseNumber
      pc(i,1:SampleNumber)=pc_freq1((i-1)*SampleNumber+1:i*SampleNumber);
end
figure(5);
plot(abs(pc(1,:)));
%mesh(abs(pc));title('脈壓 結果');
%================MTI（動目標顯示）,對消靜止目標和低速目標---可抑制雜波=================================%
for i=1:PulseNumber-1  %滑動對消，少了一個脈沖
   mti(i,:)=pc(i+1,:)-pc(i,:);
end
figure(6);
mesh(abs(mti));title('MTI 結果');
%plot(abs(mti(1,:)));
%================MTD（動目標檢測）,區(qū)分不同速度的目標，有測速作用=================================%
mtd=zeros(PulseNumber,SampleNumber);
for i=1:SampleNumber
   buff(1:PulseNumber)=pc(1:PulseNumber,i);
   buff_fft=fft(buff);
   mtd(1:PulseNumber,i)=buff_fft(1:PulseNumber)';
end
 %figure(7);plot(abs(mtd));title('MTD  結果');
 figure(7);mesh(abs(mtd));title('MTD  結果');
%==以下是為DSP程序提供回波數(shù)據(jù)、脈壓系數(shù)=====================================================%
fo=fopen('e:\coeff_fft_re.dat','wt');%頻域脈壓系數(shù)的實部
for i=1:2048
    fprintf(fo,'%f\r\n',real(coeff_fft(i)));
end
fclose(fo);
fo=fopen('e:\coeff_fft_im.dat','wt');%頻域脈壓系數(shù)的虛部
for i=1:2048
    fprintf(fo,'%f\r\n',imag(coeff_fft(i)));
end
fclose(fo);
fo=fopen('e:\echo_re.dat','wt');%8次回波的實部
for i=1:TotalNumber
    fprintf(fo,'%f\r\n',real(Echo(i)));
end
fclose(fo);
fo=fopen('e:\echo_im.dat','wt');%8次回波的虛部
for i=1:TotalNumber
    fprintf(fo,'%f\r\n',imag(Echo(i)));
end
fclose(fo);