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
|
QQ好友和群
QQ空間
騰訊微博
騰訊朋友
收藏
淘帖
頂
踩
