matlab-(含教程)基于13位巴克码和线性调频混合调制信号MTI，MTD以及CFAR的matlab仿真

%% 雷达系统仿真 %% % 发射信号为13位巴克码和线性调频混合调制的信号，线性调频的中心频率为30MHz， % 调频带宽为4MHz，每一位码宽为10微秒，发射信号的帧周期为1毫秒 % 该雷达具有数字化正交解调、数字脉冲处理、固定目标对消、动目标检测（MTD）、 % 和恒虚警（CFAR）处理等功能 clc; clear; close all; warning off; addpath(genpath(pwd)); %%%%%%%%%%%%%%% 产生雷达发射信号 %%%%%%%%%%%%% code=[1,1,1,1,1,-1,-1,1,1,-1,1,-1,1]; % 13位巴克码 tao=10e-6; % 脉冲宽度10us fc=28e6; % 调频信号起始频率 f0=30e6; % 调频信号中心频率 fs=100e6; % 采样频率 ts=1/fs; % 采样间隔 B=4e6; % 调频信号调频带宽 t_tao=0:1/fs:tao-1/fs; % 调制信号，对于线性调频来说，调制信号就是时间序列 N=length(t_tao); k=B/fs*2*pi/max(t_tao); % 调制灵敏度，也就是线性调频的步进系数 n=length(code); pha=0; s=zeros(1,n*N); for i=1:n if code(i)==1 pha=pi; else pha = 0; end s(1,(i-1)*N+1:i*N)=cos(2*pi*fc*t_tao+k*cumsum(t_tao)+pha); end t=0:1/fs:n*tao-1/fs; figure,subplot(2,1,1),plot(t,s); xlabel('t(单位：S)'),title('混合调制信号（13为巴克码+线性调频）'); s_fft_result=abs(fft(s(1:N))); subplot(2,1,2),plot((0:fs/N:fs/2-fs/N),abs(s_fft_result(1:N/2))); xlabel('频率(单位：Hz)'),title('码内信号频谱'); %%%%%%%%%%%%%%%%%%% 产生脉冲压缩系数 %%%%%%%%%%%%%%%% %--------------------- 正交解调 --------------------% N=tao/ts; n=0:N-1; s1=s(1:N); local_oscillator_i=cos(2*pi*f0/fs*n); % I路本振信号 local_oscillator_q=sin(2*pi*f0/fs*n); % Q路本振信号 fbb_i = local_oscillator_i.*s1; % I路解调 fbb_q = local_oscillator_q.*s1; % Q路解调 window=chebwin(51,40); % 50阶cheby窗的FIR低通滤波器 [b,a]=fir1(50,2*B/fs,window); fbb_i=[fbb_i,zeros(1,25)]; % 因为该FIR滤波器有25个采样周期的延迟，为了保证 % 所有有效信息全部通过滤波器，故在信号后补25个0 fbb_q=[fbb_q,zeros(1,25)]; fbb_i=filter(b,a,fbb_i); fbb_q=filter(b,a,fbb_q); fbb_i=fbb_i(26:end); % 截取有效信息 fbb_q=fbb_q(26:end); % 截取有效信息 fbb=fbb_i+j*fbb_q; %-------- 产生理想线性调频脉冲压缩匹配系数 ---------% D = B*tao; match_filter_1=ts*fliplr(conj(fbb))*sqrt(D)*2/tao; NFFT = 131126;%2^nextpow2(length(s)); match_filter_1_fft=fft(match_filter_1,NFFT); % 第一次脉冲压缩处理匹配系数 figure; subplot(2,1,1),plot(real(match_filter_1_fft)),title('脉冲压缩系数（实部）'); subplot(2,1,2),plot(imag(match_filter_1_fft)),title('脉冲压缩系数（虚部）'); N=length(s); n=0:N-1; local_oscillator_i=cos(2*pi*f0/fs*n); % I路本振信号 local_oscillator_q=sin(2*pi*f0/fs*n); % Q路本振信号 fbb_i = local_oscillator_i.*s; % I路解调 fbb_q = local_oscillator_q.*s; % Q路解调 window=chebwin(51,40); % 50阶cheby窗的FIR低通滤波器 [b,a]=fir1(50,0.5,window); fbb_i=[fbb_i,zeros(1,25)]; % 因为该FIR滤波器有25个采样周期的延迟，为了保证 % 所有有效信息全部通过滤波器，故在信号后补25个0 fbb_q=[fbb_q,zeros(1,25)]; fbb_i=filter(b,a,fbb_i); fbb_q=filter(b,a,fbb_q); fbb_i=fbb_i(26:end); % 截取有效信息 fbb_q=fbb_q(26:end); % 截取有效信息 signal=fbb_i+j*fbb_q; clear fbb_i;clear fbb_q;clear local_oscillator_i;clear local_oscillator_q; signal_fft=fft(signal,NFFT); pc_result_fft=signal_fft.*match_filter_1_fft; pc_result=ifft(pc_result_fft,NFFT); figure,plot((0:ts:length(signal)*ts-ts),pc_result(1:length(signal))); xlabel('t(单位：S)'),title('回波脉冲压缩处理结果'); t=tao*length(code); match_filter_2=2*ts*fliplr(conj(pc_result))*2/t; match_filter_2_fft=fft(match_filter_2,NFFT); % 第二次脉冲压缩处理匹配系数 figure; subplot(2,1,1),plot(real(match_filter_2_fft)),title('脉冲压缩系数（实部）'); subplot(2,1,2),plot(imag(match_filter_2_fft)),title('脉冲压缩系数（虚部）'); %%%%%%%%%%%%%%%%%%% 产生雷达回波 %%%%%%%%%%%%%%%% f_frame=1e3; % 雷达发射信号重复频率 T_frame=1/f_frame; N_echo_frame=18; f_doppler=3.5e3; % 动目标的多普勒频率 t_frame=0:ts:T_frame-ts; t_mobj=200e-6; % 动目标位置 echo_mobj_pulse=[zeros(1,t_mobj/ts),s,zeros(1,(T_frame-t_mobj)/ts-length(s))]; echo_mobj=repmat(echo_mobj_pulse,1,N_echo_frame); t_doppler=0:ts:N_echo_frame*T_frame-ts; s_doppler=cos(2*pi*f_doppler*t_doppler); s_echo_mobj=echo_mobj.*s_doppler; t_fobj=450e-6; % 固定目标位置 echo_fobj_pulse=[zeros(1,t_fobj/ts),s,zeros(1,(T_frame-t_fobj)/ts-length(s))]; echo_fobj=repmat(echo_fobj_pulse,1,N_echo_frame); t_clutter=700e-6; % 杂波位置 t_clutter_pulse=39e-6; sigma=2; % 杂波瑞利分布参数sigma t1=0:ts:t_clutter_pulse-ts; u=rand(1,length(t1)); echo_clutter=0.08*sqrt(2*log(1./u))*sigma; % 产生瑞利分布信号 s_echo_clutter_pulse=[zeros(1,t_clutter/ts),echo_clutter,... zeros(1,round((T_frame-t_clutter)/ts-length(echo_clutter)))]; s_echo_clutter=repmat(s_echo_clutter_pulse,1,N_echo_frame); s_noise=0.1*rand(1,round(N_echo_frame*T_frame/ts)); s_echo=s_echo_mobj+echo_fobj+s_echo_clutter+s_noise; %--------------------- 正交解调 --------------------% N=N_echo_frame*T_frame/ts; n=0:N-1; local_oscillator_i=cos(2*pi*f0/fs*n); % I路本振信号 local_oscillator_q=sin(2*pi*f0/fs*n); % Q路本振信号 s_echo_i = local_oscillator_i.*s_echo; % I路解调 s_echo_q = local_oscillator_q.*s_echo; % Q路解调 window=chebwin(51,40); % 50阶cheby窗的FIR低通滤波器 [b,a]=fir1(50,2*B/fs,window); s_echo_i=[s_echo_i,zeros(1,25)]; % 因为该FIR滤波器有25个采样周期的延迟，为了保证 % 所有有效信息全部通过滤波器，故在信号后补25个0 s_echo_q=[s_echo_q,zeros(1,25)]; s_echo_i=filter(b,a,s_echo_i); s_echo_q=filter(b,a,s_echo_q); s_echo_i=s_echo_i(26:end); % 截取有效信息 s_echo_q=s_echo_q(26:end); % 截取有效信息 s_echo_mf=s_echo_i+j*s_echo_q; %%%%%%%%%%%%%%%%%%%%% 脉冲压缩处理 %%%%%%%%%%%%%%%%%% s_pc_result = zeros(N_echo_frame,NFFT); for i=1:N_echo_frame s_echo_fft_result=fft(s_echo_mf(1,(i-1)*T_frame/ts+1:i*T_frame/ts),NFFT); s_pc_fft_1=s_echo_fft_result.*match_filter_1_fft; s_pc_fft_2=s_pc_fft_1.*match_filter_2_fft; s_pc_result(i,:)=ifft(s_pc_fft_2,NFFT); end s_pc_result_1=s_pc_result'; s_pc_result_1=reshape(s_pc_result_1,1,N_echo_frame*NFFT); figure,subplot(2,1,1),plot((0:ts:N_echo_frame*NFFT*ts-ts),real(s_pc_result_1)); xlabel('t(单位：S)'),title('脉冲压缩处理后结果（实部）'); subplot(2,1,2),plot((0:ts:N_echo_frame*NFFT*ts-ts),imag(s_pc_result_1)); xlabel('t(单位：S)'),title('脉冲压缩处理后结果（虚部）'); %%%%%%%%%%%%%%%%%%%%% 固定杂波对消处理 %%%%%%%%%%%%%%%%%% S_MTI_result = zeros(N_echo_frame-2,NFFT); for i=1:N_echo_frame-2 S_MTI_result(i,:)=s_pc_result(i,:)+s_pc_result(i+2,:)-2*s_pc_result(i+1,:); end S_MTI_result_1=S_MTI_result'; S_MTI_result_1=reshape(S_MTI_result_1,1,(N_echo_frame-2)*NFFT); figure,subplot(2,1,1),plot((0:ts:(N_echo_frame-2)*NFFT*ts-ts),real(S_MTI_result_1)); xlabel('t(单位：S)'),title('固定杂波对消后结果（实部）'); subplot(2,1,2),plot((0:ts:(N_echo_frame-2)*NFFT*ts-ts),imag(S_MTI_result_1)); xlabel('t(单位：S)'),title('固定杂波对消后结果（虚部）'); clear S_MTI_result_1; %%%%%%%%%%%%%%%%%%%%% MTD处理和求模 %%%%%%%%%%%%%%%%%% S_MTD_result_1 = fft(S_MTI_result,N_echo_frame-2); S_MTD_result=abs(max(S_MTD_result_1)); figure,plot((0:ts:NFFT*ts-ts),S_MTD_result); xlabel('t(单位：S)'),title('MTD处理后求模结果（信号最大通道）'); %%%%%%%%%%%%%%%%%%%%% CFAR处理 %%%%%%%%%%%%%%%%%% cfar_result = zeros(1,NFFT); cfar_result(1,1)=S_MTD_result(1,1)/(sqrt(2)/pi*mean(S_MTD_result(1,2:17))); % 第1点恒虚警处理的噪声均值由其后面的16点的噪声决定 for i=2:16 % 第2点到第16点恒虚警处理的噪声均值由其前面和后面16点的噪声共同决定 noise_mean=sqrt(2)/pi*(mean(S_MTD_result(1,1:i-1))+mean(S_MTD_result(1,i+1:i+16)))/2; cfar_result(1,i)=S_MTD_result(1,i)/noise_mean; end for i=17:NFFT-17 % 正常的数据点的恒虚警处理的噪声均值由其前面和后面各16点的噪声中的大者决定 noise_mean=sqrt(2)/pi*max(mean(S_MTD_result(1,i-16:i