首页/数据结构与算法
FMCW信号面阵rma成像算法
算法
matlab
最近在复现mimo面阵的RMA成像代码,但是一直不明白将如何表述Kx、Ky波数域该怎么用matlab代码实现,然后我也写了一篇博客记录了复现这篇代码的过程,但是中途有几步结果不太对,我博客的链接在最下方,如下为matlab代码部分:
%复现IEEE上的fmcw成像
clc;clear;close all;
c=3e8; %光速
Fc=77e9; %起始频率
Fs=5e6; %采样频率(Hz)
Ts=1/Fs; %采样间隔
K=80e12; %调频斜率
gamma=2*pi*K; %调频斜率的角频率形式
num_sample=256; %快时间采样点数
B=K*(Ts*(num_sample-1));%扫频带宽
F_end = Fc+B; %截止频率
F0 = (F_end+Fc)/2; %中心频率
Tp=(num_sample-1)*Ts;%扫频时间
% t=linspace(0,Tp,num_sample); %两种时间表述
t=linspace(-Tp/2,Tp/2,num_sample);
Rres=c/(2*B);%距离分辨率
range=(0:num_sample/2-1)*Rres; %距离
%% 单通道成像
%% 生成回波数据
target=[0.1 0 1;
-0.1 -0.1 1;
0.1 0.1 1;
-0.1 0.1 1
-0.15 0.1 1];%目标的位置
size_T=size(target);
target_num=size_T(1);
Xm=0.4;Nx=100; %x方向发射天线数量,也是成像网格的数量
Ym=0.4;Ny=100; %y方向发射天线数量,也是成像网格的数量
T_x = [-0.4 0.4]; %发射天线x坐标
T_y = [-0.4 0.4]; %发射天线y坐标
T_z = [ 0 0]; %发射天线z坐标
imag_area_x=linspace(-0.2,0.2,Nx); %x方向成像区域 也是接收天线X坐标
imag_area_y=linspace(-0.2,0.2,Ny); %y方向成像区域 也是接收天下Y坐标
Echo=zeros(length(T_x),length(T_y),length(imag_area_x),length(imag_area_y),num_sample);%五维回波矩阵
%% 生成回波矩阵 回波生成完全按照论文中的回波信号公式来表述
disp('生成回波数据中...');
tic;
for l = 1 : length(T_y)
for m = 1 : length(T_x)
for i=1:length(imag_area_y)
for j=1:length(imag_area_x)
for k=1:target_num
% R_st = Distance(target(k,:),imag_area_x(i),imag_area_y(j));
% R_sr = Distance(target(k,:),T_x(m),T_y(l));
S_r = exp(-1j*2*pi*(F0+K*t)/c*sqrt((T_x(m)-target(k,1))^2 +(T_y(l)-target(k,2))^2+(0 - target(k,3))^2)) .*...
exp(-1j*2*pi*(F0+K*t)/c*sqrt((imag_area_x(i)-target(k,1))^2 +(imag_area_y(j)-target(k,2))^2+(0 - target(k,3))^2));
% tx_signal1=mmwave(F0,Fs,S,num_sample,0);
% rx_signal1=mmwave(F0,Fs,S,num_sample,R);
% Echo(l,m,i,j,:)=squeeze(Echo(l,m,i,j,:))'+real(tx_signal1).*real(rx_signal1)+1j*(imag(tx_signal1).*imag(rx_signal1));
% Echo(l,m,i,j,:) = squeeze(Echo(l,m,i,j,:))' + tx_signal1 .* conj(rx_signal1);
Echo(l,m,i,j,:) = Echo(l,m,i,j,:) + reshape(S_r,[1,1,1,1,num_sample]);
end
end
end
end
end
time1=toc;
fprintf('回波数据生成完毕,共用时: %.2f 秒\n', time1);
%% RMA
%%方位向波数分辨率
Echo_1 = squeeze(Echo(1,1,:,:,:)); %对应论文中提取其中一个发射阵元对应所有接收阵元的数据
disp('RMA算法成像中...');
tic;
% t=linspace(-Tp/2,Tp/2,num_sample);
% t=linspace(0,Tp,num_sample);
Rxsm=0.2; %方位x轴最大采集距离-Rxsm~Rxsm
Rysm=0.2; %方位y轴最大采集距离-Rysm~Rysm
wSx=2*pi/(2*Rxsm); %波束分辨率 = 2*pi/(成像区域孔径长度)
wSy=2*pi/(2*Rysm);
% kX=wSx*(-(Nx-1)/2:(Nx-1)/2); %波数分量x方向
% kY=wSy*(-(Ny-1)/2:(Ny-1)/2); %波数分量y方向
kX=wSx*linspace(-Nx/2,Nx/2,Nx); %波数分量x方向
kY=wSy*linspace(-Ny/2,Ny/2,Ny); %波数分量y方向
kX3D=repmat(kX,[Nx,1,num_sample]); %kX3D为沿着列方向为x方向,所以需要沿着行和Z进行拓展
kY3D=repmat(kY.',[1,Ny,num_sample]); %kY3D为沿着列方向所以首先需要进行转置为沿着行方向,然后沿着列方向和Z方向进行拓展
% kX3D=repmat(kX.',[1,Ny,num_sample]);
% kY3D=repmat(kY,[Nx,1,num_sample]);
% K=gamma*t/c+W0/c; %K为空间中的波数 第一种写法
K = 2*pi*(F0+K*t)/c; %第二种写法
k1D(1,1,:)=K;
k3D=repmat(k1D,[Nx,Ny,1]); %空间波数矩阵
%天线位置 和 切片位置
Za = 0;
Zs = 1;
%矫正相位
disp('计算相位因子...');
Phase0=zeros(Nx,Ny,num_sample);
phase0 = exp(-1j*sqrt(k3D.^2 - kX3D.^2 - kY3D.^2)*(Za - Zs)); %论文中的第一个相位补偿因子
%方位fft
disp('二维方位向fft...');
% data_fft2D= fftshift(fft(ifftshift((fftshift(fft(ifftshift(Echo_1,1),[],1),1)),2),[],2),2);%2Dfft
data_fft2D= fftshift(fft2(Echo_1));
figure(1);
imagesc(kX,kY,real(phase0(:,:,1)));
title('第一个矫正因子');
%% 相位矫正
disp('相位补偿一...');
data_phase_correct = data_fft2D .* phase0;
disp('二维方位向ifft...');
% data_ifft2D=ifftshift(ifft(fftshift((ifftshift(ifft(fftshift(data_phase_correct,1),[],1),1)),2),[],2),2);
data_ifft2D = ifft2(ifftshift(data_phase_correct));
figure(2);
imagesc(kX,kY,abs(data_ifft2D(:,:,25)));
title('矫正后二维方位iFFT结果');
disp('距离维fft...');
data_ifft_3D = fft(data_ifft2D,[],3);
figure(3);
imagesc(kX,kY,abs(data_ifft_3D(:,:,25)));
title('距离维fft');
disp('相位补偿二...');
Rti = zeros(length(imag_area_x),length(imag_area_y));
for i = 1 : length(imag_area_y) %计算Rti到当前切片的所有成像网格的距离
for j=1:length(imag_area_x)
Rt1 = sqrt((T_x(1)-imag_area_x(j))^2 + (T_y(1)-imag_area_y(i))^2 + (Za-Zs)^2);
Rti(i,j) = Rt1;
end
end
data_O = data_ifft_3D .* exp(1j*2*pi*F0* Rti/c);
figure(4);
imagesc(imag_area_x,imag_area_y,abs(data_O(:,:,1)));
title('最终结果');
% for i = 1 : 1 : 256
% figure(i);
% imagesc(imag_area_x,imag_area_y,abs(data_O(:,:,i)));
% end
time2=toc;
fprintf('成像完成,共用时: %.2f 秒\n', time2);
https://blog.youkuaiyun.com/qq_34997585/article/details/153648340?spm=1011.2124.3001.6209
收起
最新发布