%==產生LFM信號-分析其頻譜====% 
clc;clear all;close all;warning off;              
snr=0;                           %信噪比
fc=0e10;                         % 載波頻率;
c=3e8;                           % 真空光速;
Tr=750e-6;                       % 脈沖重復周期s;
fd=1.6e6;                          % DOPPLOR 頻率
td=2e-6;                         % 時間延遲s
Tp=1e-5;                         % LFM調制脈沖寬度s；
B=8e6;                           % LFM調制脈沖帶寬hz；
fs=2*B;                          % LFM采樣頻率；
t=-Tp/2:1/fs:Tp/2;N=length(t);
lfm=10^(snr/10)*exp(j*2*pi*(fc*t+0.5*B/Tp*t.^2));
fft_lfm=fftshift(fft(lfm,N));
LFM_Back=10^(snr/10)*exp(j*2*pi*(fc*(t-td)+0.5*B/Tp*(t-td).^2)).*exp(-j*2*pi*fd*t);;
figure(1)
subplot(211),plot(-Tp/2+Tp/N:Tp/N:Tp/2,real(lfm));xlabel('LFM脈沖寬度(s)'); ylabel('幅度');title('實LFM信號波形');axis([-Tp,Tp,-2,2])
subplot(212),plot(-Tp/2+Tp/N:Tp/N:Tp/2,real(LFM_Back));xlabel('LFM脈沖寬度(s)'); ylabel('幅度');title('實LFM信號回波波形');axis([-Tp,Tp,-2,2])

figure(2),
     subplot(211),plot(-0.5+1/N:1/N:0.5,abs(fft_lfm(1:N)));xlabel('歸一化頻率'); ylabel('幅度');title('LFM信號頻域幅度');
     subplot(212),plot(-0.5+1/N:1/N:0.5,angle(fft_lfm(1:N)));xlabel('歸一化頻率'); ylabel('相角');title('LFM信號頻域相角');
%%%%%%%%%%%%%%%%%%%%%%%%%%% 頻域匹配濾波 % 放大可查看半功率寬度=1/B
out_mat=fftshift(ifft(fft((LFM_Back)).*conj(fft((lfm))))) ;
figure(3)
subplot(2,1,1),plot(-Tp/2:1/fs:Tp/2,20*log10(abs(out_mat)/max(abs(out_mat))));title('匹配輸出LFM');axis([-Tp/2,Tp/2,-200,0]);
%%%%%%%%%%%%%%%%%%%%%%%%%%% 頻域匹配濾波 % 放大可查看半功率寬度=1/B  , 加窗
out_mat=fftshift(ifft(fft((LFM_Back)).*conj(fft((lfm.*hamming(N)'))))) ;
subplot(2,1,2),plot(-Tp/2:1/fs:Tp/2,20*log10(abs(out_mat)/max(abs(out_mat))));title('匹配輸出LFM');axis([-Tp/2,Tp/2,-200,0]);

%%%%%%%%%%%%%%%%%%%%%%%%%%% 時域匹配濾波 
matcher=conj(fliplr(lfm));                % 時域匹配濾波 
weight=conj(fliplr(lfm.*hamming(N)'));    % 時域匹配濾波 加窗
out_matcher=conv(LFM_Back,matcher);
out_weight=conv(LFM_Back,weight);
figure(4)
subplot(2,1,1),plot(-Tp:1/fs:Tp,20*log10(abs(out_matcher)/max(abs(out_matcher))));title('時間匹配輸出LFM');axis([-Tp/2,Tp/2,-200,0]);
subplot(2,1,2),plot(-Tp:1/fs:Tp,20*log10(abs(out_weight)/max(abs(out_weight))));title('時間加窗匹配輸出LFM');axis([-Tp/2,Tp/2,-200,0]);



%%%%%%%%%%%%%%%%%%%%%%%%%%% 頻率時間聯合估計
fdd=(B/Tp)*Tp;
ii=0;
NN=20;% 頻域采樣數目
out_mat=zeros(1,2*N-1);
for fdx=-fdd:fdd/NN:fdd;
    ii=ii+1;
    matcher=conj(fliplr(lfm)).*exp(-j*2*pi*fdx*t); 
    out_mat(ii,:)=conv(LFM_Back,matcher);
end
figure(5);[F,T]=meshgrid(-Tp:1/fs:Tp,-fdd:fdd/NN:fdd);mesh(F,T,abs(out_mat));title('頻率時間聯合估計');axis([-Tp,Tp,-fdd,fdd,min(min(abs(out_mat))),max(max(abs(out_mat)))]);grid on;
figure(6);contour(F,T,abs(out_mat),4);title('頻率時間聯合估計-等高線圖');xlabel('時延td (s)');ylabel('多普勒頻移fd (Hz)');grid on;
figure(7)
subplot(2,1,1);plot(-fdd:fdd/NN:fdd,abs(out_mat(:,N)));title('|X(0,fd)| 多普勒頻移估計');xlabel('多普勒頻移fd (Hz)');ylabel('|x(0,fd)|');grid on;axis([-fdd,fdd,min(abs(out_mat(:,N))),max(abs(out_mat(:,N)))]);
subplot(2,1,2);plot(-Tp:1/fs:Tp,abs(out_mat(NN+1,:)));title('|X(td,0)| 延遲估計 ');xlabel('時延td (s)');ylabel('|x(td,0)|');grid on;axis([-Tp/2,Tp/2,-200,0]);axis([-Tp,Tp,min(abs(out_mat(NN+1,:))),max(abs(out_mat(NN+1,:)))]);




%%%%%%%%%%%%%%%%%%%%%%%%%%% 模糊函數,等高線圖(1)－基本原理
fdd=(B/Tp)*Tp;
ii=0;
NN=20;% 頻域采樣數目
out_mat=zeros(1,2*N-1);
for fdx=-fdd:fdd/NN:fdd;
    ii=ii+1;
    matcher=conj((lfm)).*exp(-j*2*pi*fdx*t); 
    out_mat(ii,:)=conv(lfm,matcher);
end
figure(8);[F,T]=meshgrid(-Tp:1/fs:Tp,-fdd:fdd/NN:fdd);mesh(F,T,abs(out_mat));axis([-Tp,Tp,-fdd,fdd,min(min(abs(out_mat))),max(max(abs(out_mat)))]);grid on;
figure(9);contour(F,T,abs(out_mat),16);title('等高線圖');xlabel('時延td (s)');ylabel('多普勒頻移fd (Hz)');grid on;
figure(10)
subplot(2,1,1);plot(-fdd:fdd/NN:fdd,abs(out_mat(:,N)));title('|X(0,fd)| 多普勒頻移靈敏度');xlabel('多普勒頻移fd (Hz)');ylabel('|x(0,fd)|');grid on;axis([-fdd,fdd,min(abs(out_mat(:,N))),max(abs(out_mat(:,N)))]);
subplot(2,1,2);plot(-Tp:1/fs:Tp,abs(out_mat(NN+1,:)));title('|X(td,0)| 延遲靈敏度 ');xlabel('時延td (s)');ylabel('|x(td,0)|');grid on;axis([-Tp/2,Tp/2,-200,0]);axis([-Tp,Tp,min(abs(out_mat(NN+1,:))),max(abs(out_mat(NN+1,:)))]);

%%%%%%%%%%%%%%%%%%%%%%%%%%% 模糊函數,等高線圖(2)－變換公式
% ①|χ(τ,fd)|=A2T(1-|τ/T|)exp(-jωd*τ/2)*Sa[(μτ+ωd)(1-|τ/T|*T/2)]  |τ|<T
%                ＝0		當|τ|〉T時
%   ②|χ(τ,0)|= A2T(1-|τ/T|) Sa[(μτ*T/2)(1-|τ/T|)]。
%   ③|χ(0,fd)|=A2T*Sa(ωd*T/2)。
ii=0;
for t=-Tp:1/fs:Tp;  ii=ii+1;kk=0;
    for fdx=-fdd:fdd/NN:fdd; 
        kk=kk+1;
        wd=2*pi*fd;
        out_mat2(ii,kk)=(1-abs(t)/Tp)*exp(-j*wd*t/2)*sinc(((B/Tp)*t+wd)*(1-abs(t)/Tp)*Tp/2);
    end
end
figure(11);[F,T]=meshgrid(-Tp:1/fs:Tp,-fdd:fdd/NN:fdd);mesh(F,T,abs(out_mat));axis([-Tp,Tp,-fdd,fdd,min(min(abs(out_mat))),max(max(abs(out_mat)))]);grid on;
figure(12);contour(F,T,abs(out_mat),16);title('等高線圖');xlabel('時延td (s)');ylabel('多普勒頻移fd (Hz)');grid on;
figure(13)
subplot(2,1,1);plot(-fdd:fdd/NN:fdd,abs(out_mat(:,N)));title('|X(0,fd)| 多普勒頻移靈敏度');xlabel('多普勒頻移fd (Hz)');ylabel('|x(0,fd)|');grid on;axis([-fdd,fdd,min(abs(out_mat(:,N))),max(abs(out_mat(:,N)))]);
subplot(2,1,2);plot(-Tp:1/fs:Tp,abs(out_mat(NN+1,:)));title('|X(td,0)| 延遲靈敏度 ');xlabel('時延td (s)');ylabel('|x(td,0)|');grid on;axis([-Tp/2,Tp/2,-200,0]);axis([-Tp,Tp,min(abs(out_mat(NN+1,:))),max(abs(out_mat(NN+1,:)))]);