%===========================================================================================%
%            該程序完成艦載雷達中頻接收信號的產生(為脈壓產生）                   %
%===========================================================================================%
close all; 
clear all;
%信號的產生方法按照雷達回波形式
%===================================================================================%
%                                    雷達參數                                       %
%===================================================================================%
C=3.0e8;  %光速(m/s)
RF=3.140e9;  %雷達射頻
Lambda=C/RF;%雷達工作波長
PulseNumber=2;   %回波脈沖數 
BandWidth=2.0e6;  %發射信號帶寬
TimeWidth=42.0e-6; %發射信號時寬
PRT=230e-6;   % 雷達發射脈沖重復周期(s)
PRF=1/PRT;   %雷達發射脈沖重復頻率
Fs=20.0e6;  %采樣頻率
NoisePower=-12;%(dB);
%---------------------------------------------------------------%
SampleNumber=fix(Fs*PRT);%計算一個脈沖周期的采樣點數；
TotalNumber=SampleNumber*PulseNumber;%總的采樣點數,920；
BlindNumbe=fix(Fs*TimeWidth);%計算一個脈沖周期的盲區-遮擋樣點數；
%===================================================================================%
%                                    目標參數                                       %
%===================================================================================%
TargetNumber=3;%目標個數
SigPower(1:TargetNumber)=[1 .5 .3];%目標功率,無量綱
TargetDistance(1:TargetNumber)=[7000 14000 23000];%目標距離,單位m
DelayNumber(1:TargetNumber)=fix(Fs*2*TargetDistance(1:TargetNumber)/C);
TargetVelocity(1:TargetNumber)=[0 30 100];%目標徑向速度 單位m/s
TargetFd(1:TargetNumber)=2*TargetVelocity(1:TargetNumber)/Lambda  %計算目標多卜勒
%====================================================================================%
%                                   產生線性調頻信號 即產生線性調頻方式的發射信號                                    %
%====================================================================================%
%線性調頻信號的突出優點是匹配濾波器對回波信號的多普勒頻移不敏感，即使回波信號有較大的多普勒頻移，原來的匹配濾波器仍能起到脈沖壓縮的使用。
number=fix(Fs*TimeWidth);%BlindNumbe
if rem(number,2)~=0%不等于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(Chirp);   %coeff為匹配濾波的系數，conj為取共軛函數
%Chirp_fft=fft(Chirp,128);
figure(20);
subplot(2,1,1);plot(real(Chirp));title('LFM信號實部');
 subplot(2,1,2);plot(imag(Chirp));title('LFM信號虛部');
figure(10);plot(abs(fft(Chirp)));title('LFM信號');

%-------------------------產生目標回波串------------------------%
% 目標回波串信號=發射信號經過一定的延時后在加上多普勒平移。即SignalAll=Chirp(t-t0)*exp(-jwt)  
SignalAdd=zeros(1,TotalNumber);%所有脈沖的信號
for k=1:TargetNumber% 依次產生各個目標
   SignalTemp=zeros(1,SampleNumber);% 一個脈沖
   SignalTemp(DelayNumber(k)+1:DelayNumber(k)+number)=sqrt(SigPower(k))*Chirp;%一個脈沖的某個目標
   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;
   SignalAdd=SignalAdd+Signal;
end
SignalAll=SignalAdd;   %SignalAll為模擬的目標回波信號，沒有加入噪聲
 figure(1);
  subplot(2,1,1);plot(real(SignalAll));title('目標信號的實部');
  subplot(2,1,2);plot(imag(SignalAll));title('目標信號的虛部');
 grid on;zoom on;

%====================================================================================%
%                                   產生系統噪聲信號                                  %
%====================================================================================%
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');
 
%================================時域脈壓=================================%
%===========時域脈沖壓縮直接對雷達回波信號進行卷積運算，其算式如下===========%
%   s(n)=s1(n)+jsQ(n);h(n)=hI(n)+jhQ(n)
%   y(n)=s(n)*h(n)     (1)
%   式中，s(n)為總的回波信號；h(n)為匹配濾波器的沖激響應信號；y(n)為經匹配濾波器壓縮后的輸出信號。
%   對線性調頻信號的脈沖壓縮處理，就是讓信號通過一個與其相匹配的濾波器實現的
pc_time0=conv(Echo,coeff);%conv為卷積函數
figure(3);plot(abs(pc_time0));title('時域脈壓結果的幅度,有暫態點');
pc_time1=pc_time0(number:TotalNumber+number-1);%去掉暫態點 number-1個

fo=fopen('td_pc.dat','wt');
for i=1:920
    fprintf(fo,'%f\r\n',abs(pc_time0(i)));
end
fclose(fo);


%================================頻域脈壓=================================%
%頻域脈沖壓縮先對輸入的總的回波序列進行FFT變換，然后點乘以匹配濾波器沖擊響應的離散譜，再用逆FFT還原成壓縮后的時間離散信號。其算式如下：
% %S（k）=FFT(s(n));H(k)=FFT(h(n))
%y(n)=IFFT(S(k)×H(k))=IFFT(FFT(s(n))×FFT(h(n)))   
 Echo_fft=fft(Echo,TotalNumber+number-1);
 coeff_fft=fft(coeff,TotalNumber+number-1);%匹配濾波器沖擊響應的離散譜
Echo_fft=fft(Echo,1024);%理應進行TotalNumber+number-1點FFT,但為了提高運算速度,進行了1024點的FFT
coeff_fft=fft(coeff,1024);
pc_fft=Echo_fft.*coeff_fft;
pc_freq0=ifft(pc_fft);
figure(4);plot(abs(pc_freq0));title('頻域脈壓結果的幅度,有暫態點和后填充點,紅色為時域頻域脈壓的差別');
hold on;plot(abs(pc_time0(1:TotalNumber+number-1)-pc_freq0(1:TotalNumber+number-1)),'r');
pc_freq1=pc_freq0(number:TotalNumber+number-1);%去掉暫態點 number-1個,后填充點若干(1024-number+1-TotalNumber=21個)

fo=fopen('fd_pc.dat','wt');%存放頻域脈壓數據
for i=1:920
    fprintf(fo,'%f\r\n',abs(pc_freq0(i)));
end
fclose(fo);

%================MTI,對消靜止目標和低速目標---相當于雜波=================================%
%目標信號=總的回波信號(有多普勒頻移)-發射信號(無多普勒頻移)，這樣進行MTI
mti=pc_time1(SampleNumber+1:2*SampleNumber)-pc_time1(1:SampleNumber);%
figure(5);plot(abs(mti));title('MTI 結果');

%load dsp_mti_result.dat;
%dsp_mti_result;
%error=dsp_mti_result'-mti;
%figure(6);plot(abs(error));title('MTI 結果誤差');

fo=fopen('matlab_mti_result.dat','wt');  %存取mti數據
for i=1:460   %   2個周期=920點
   fprintf(fo,'%f\r\n',abs(mti(i)));
end
fclose(fo);


%save matlab_mti_result.dat     abs(mti)  -ascii -double

fo=fopen('coeff_fft_re.dat','wt');
for i=1:1024
    fprintf(fo,'%f\r\n',real(coeff_fft(i)));%將Matlab產生的LFM系數（84點）變換到頻域（1024點），實部存入文件
end
fclose(fo);
fo=fopen('coeff_fft_im.dat','wt');
for i=1:1024
    fprintf(fo,'%f\r\n',imag(coeff_fft(i)));%虛部存入文件
end
fclose(fo);

fo=fopen('Echo_re.dat','wt');%Matlab產生的LFM信號（920點），實部存入文件
for i=1:920
    fprintf(fo,'%f\r\n',real(Echo(i)));
end
fclose(fo);
fo=fopen('Echo_im.dat','wt');%Matlab產生的LFM信號（920點），虛部存入文件
for i=1:920
    fprintf(fo,'%f\r\n',imag(Echo(i)));
end
fclose(fo);