产生lfm信号-分析其频谱-运用频域匹配滤波器.m

线性调频信号频谱及匹配滤波方面的源代码

%==产生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,:)))]);