ISAR成像算法仿真

ID:168257 · 发表于 2017-3-5 10:26

ISAR成像算法仿真
本程序是包括从ISAR 雷达系统获取数据  并且进行距离向脉冲压缩  得到脉压之后的距离时间图像
程序源码如下

tic
close all
clear all
%% file
file_rawdata0 = 'G:\ID118_8H76M';
%% parameter
Tr       = 500e-6;
Br       = 1200e6;
f0       = 34.6e9; %载频
PRF    = 1/Tr;
Fsr    = 50e6;                      %A/D采样率
head_len  = 64;                      %包头信息长度128个字节
pulse_len =  round(Tr*Fsr);             %采样点数，每个点2个Byte =16 bit：;  pulse_len = round(Tr*Fsr);
pkg_len = pulse_len+head_len;
rg_offset = 0;                      %距离向偏移
az_offset = PRF*74;    %125e3*1                %方位向偏移

c = 299792458;
Kr = Br/Tr;
Rbin = c/2*Fsr/pulse_len/abs(Kr);    % 去调频接收的距离门计算方法
rg_len = round(Tr*Fsr); %500us 时候值为25000

%% reading and processing data

pulse_num = PRF*3;%125e3*5;
fid_in0 = fopen(file_rawdata0,'rb','b'); %二进制方式打开文件 b？
fseek(fid_in0,az_offset*pkg_len*2,'bof');    %将文件指针指向pkg_len*2位置，bof从文件开头偏移  可以自己设置方位向偏移az_offset，跳过前面的脉冲数；若取某段时间，计算读取该时间的脉冲数
data0 = fread(fid_in0,pkg_len*pulse_num,'int16');  %读二进制文件16位整型数
fclose(fid_in0);
% figure;plot((data0(1:2000)))
data0 = reshape(data0,pkg_len,pulse_num);  %读二进制文件16位整型数
size(data0)
data0(1:64,:) = [];
size(data0)
%%
dt=1/PRF;
df=PRF/pulse_num;
f=df*(-pulse_num/2:pulse_num/2-1);
t=dt*(0:pulse_num-1);
data0 = fft(data0);
data0 = data0(1:rg_len/2,:);
%  data0(1:30,:) = 0;
%%
figure
%下面这句是专门用于 ID49 对于ID49 螺旋桨位置位于180:200  后面微多普勒分析也是对应这个位置
%imagesc(t,Rbin*(179:pulse_len/2-1) ,20*log10(abs(data0(180:200,:))));
%ID41 对于ID41 右侧螺旋桨位置位于120:140  左侧对应155:165 后面微多普勒分析也是对应这个位置
% ID45 对于ID45 螺旋桨位置位于280:330  后面微多普勒分析也是对应这个位置
% ID35 对于ID35 铁锹位置位于80:110  后面微多普勒分析也是对应这个位置
% ID61 对于ID45 右侧螺旋桨位置位于145:155  后面微多普勒分析也是对应这个位置
% ID76 对于ID76 右侧螺旋桨位置位于640:960  后面微多普勒分析也是对应这个位置
max_data0 = max(max(abs(data0)));
%  imagesc(t,Rbin*(79:109) ,20*log10(abs(data0(80:110,:))/(max_data0)));  %t,Rbin*(0:pulse_len/2-1) ,
% figure; imagesc(t,Rbin*(639:959) ,20*log10(abs(data0(640:960,:))/(max_data0)));set(gca,'CLim', [-50 0]);
imagesc(t,Rbin*(0:rg_len/2-1),abs(data0)); % t,Rbin*(0:rg_len/2-1), abs(data0)
axis xy
clim = get(gca,'CLim');
%  set(gca,'CLim',clim(2) + [-30 0]);
colorbar
title('时间-距离像');
xlabel('时间 (s)')
ylabel('距离 (m)')

%%
%时频分析

  data1  = sum(data0(180:195,:));    figure;plot(abs(data1)); figure;plot(abs(fftshift(fft(data1))));
% data1  = data0(289,:);figure;plot(abs(data1)); figure;plot(abs(fftshift(fft(data1))));
% % clear data0
% % % [mm,nn] = size(data11)
% % % Nn    = 10 ;  % 降采样倍数
% % % data1 = data11(1,1:10:nn); size(data1)
% % % figure;plot(abs(data1));
% % h1=window('hanning',65);
% % [TF,t,f]=tfrstft(data1',t,2^15,h1);
% % % TF    = tfrstft(data1');
% % % TF    = fftshift(TF,1); %将频率零频移到中心fftshift
% % figure
% %  imagesc(abs(TF));
% % % imagesc(20*log10(abs(TF)));
% % axis xy;
% % % title('时频分析')
%%
%时频分析
x  = data1;
np = length(x);
wd = 512;
wdd2 = wd/2;
wdd8 = wd/8;
ns = floor(np/wd);

% calculate time-frequency micro-Doppler signature
disp('Calculating segments of TF distribution ...')
for k = 1:ns
disp(strcat('  segment progress: ',num2str(k),'/',num2str(round(ns))))
sig(1:wd,1) = x(1,(k-1)*wd+1:(k-1)*wd+wd);
   TMP = stft(sig,16);
%    TMP = tfrspwv(sig);
TF2(:,(k-1)*wdd8+1:(k-1)*wdd8+wdd8) = TMP(:,1:8:wd);
end
TF = TF2;
disp('Calculating shifted segments of TF distribution ...')
TF1 = zeros(size(TF));
for k = 1:ns-1
disp(strcat('  shift progress: ',num2str(k),'/',num2str(round(ns-1))))
sig(1:wd,1) = x(1,(k-1)*wd+1+wdd2:(k-1)*wd+wd+wdd2);
   TMP = stft(sig,16);
%    TMP = tfrspwv(sig);
TF1(:,(k-1)*wdd8+1:(k-1)*wdd8+wdd8) = TMP(:,1:8:wd);
end
disp('Removing edge effects ...')
for k = 1:ns-1
TF(:,k*wdd8-8:k*wdd8+8) = ...
      TF1(:,(k-1)*wdd8+wdd8/2-8:(k-1)*wdd8+wdd8/2+8);
end

% display final time-frequency signature
figure;
colormap(jet(256))
imagesc(t,f,20*log10(fftshift(abs(TF),1)+eps))
xlabel('时间 (s)')
ylabel('多普勒频率 (Hz)')
title('旋转螺旋桨三个叶片微多普勒特征')
axis xy
clim = get(gca,'CLim');
set(gca,'CLim',clim(2) + [-50 0]);
colorbar

帐号		自动登录	找回密码
密码			立即注册