自制合成孔径雷达(1) 后处理程序在Octave下运行

我最近看到一个麻省理工学院的开放课程，用一些简易电路来实现一个雷达，可以测距、测速也可以做合成孔径雷达。硬件电路用adc+单片机+usb转接实现，然后传输给电脑，电脑上c#程序做实时处理。但是这个教程资料还不是很完善，我没找到单片机和c#代码。然后我上MIT opencourseware网站上找到了老版本的资料。

链接: https://pan.baidu.com/s/1bvSZxAIw3A-79efm4WifsA  密码: 4qr7

老版本与新版本的区别是：

1.硬件上老版本全部都是模拟电路，然后通过声卡线路输入传输给电脑，也就是说adc是利用了电脑内部的，不是自己搭的。

2.软件上老版本用matlab实现，是把音频文件录下来做后处理的，但我觉得我可以改写为python，然后实时处理。

在搭硬件电路前，我想先试试软件行不行，毕竟试一下软件没啥成本，而且如果软件不行那硬件也白搭。

而且我的电脑比较慢，我不打算装庞大的matlab，我用了octave。

下面几个图就是我用octave跑的3个范例程序。多普勒测速是可以直接跑的。测距和SAR是需要改一下代码，其实也很简单只需要把&改为&&就行。有一点要注意多普勒测速程序运行起来比较快可能几秒就有结果，后面两个程序大概要等个几分钟才会有结果。

doppler:

 

range:

read_data_RTI.m

%MIT IAP Radar Course 2011
%Resource: Build a Small Radar System Capable of Sensing Range, Doppler, 
%and Synthetic Aperture Radar Imaging 
%
%Gregory L. Charvat
 
%Process Range vs. Time Intensity (RTI) plot
 
%NOTE: set up-ramp sweep from 2-3.2V to stay within ISM band
%change fstart and fstop bellow when in ISM band
 
clear all;
close all;
 
%read the raw data .wave file here
[Y,FS,NBITS] = wavread('running_outside_20ms.wav');
 
%constants
c = 3E8; %(m/s) speed of light
 
%radar parameters
Tp = 20E-3; %(s) pulse time
N = Tp*FS; %# of samples per pulse
fstart = 2260E6; %(Hz) LFM start frequency for example
fstop = 2590E6; %(Hz) LFM stop frequency for example
%fstart = 2402E6; %(Hz) LFM start frequency for ISM band
%fstop = 2495E6; %(Hz) LFM stop frequency for ISM band
BW = fstop-fstart; %(Hz) transmti bandwidth
f = linspace(fstart, fstop, N/2); %instantaneous transmit frequency
 
%range resolution
rr = c/(2*BW);
max_range = rr*N/2;
 
%the input appears to be inverted
trig = -1*Y(:,1);
s = -1*Y(:,2);
clear Y;
 
%parse the data here by triggering off rising edge of sync pulse
count = 0;
thresh = 0;
start = (trig > thresh);
for ii = 100:(size(start,1)-N)
    if start(ii) == 1 && mean(start(ii-11:ii-1)) == 0
        %start2(ii) = 1;
        count = count + 1;
        sif(count,:) = s(ii:ii+N-1);
        time(count) = ii*1/FS;
    end
end
%check to see if triggering works
% plot(trig,'.b');
% hold on;si
% plot(start2,'.r');
% hold off;
% grid on;
 
%subtract the average
ave = mean(sif,1);
for ii = 1:size(sif,1);
    sif(ii,:) = sif(ii,:) - ave;
end
 
zpad = 8*N/2;
 
%RTI plot
figure(10);
v = dbv(ifft(sif,zpad,2));
S = v(:,1:size(v,2)/2);
m = max(max(v));
imagesc(linspace(0,max_range,zpad),time,S-m,[-80, 0]);
colorbar;
ylabel('time (s)');
xlabel('range (m)');
title('RTI without clutter rejection');
 
%2 pulse cancelor RTI plot
figure(20);
sif2 = sif(2:size(sif,1),:)-sif(1:size(sif,1)-1,:);
v = ifft(sif2,zpad,2);
S=v;
R = linspace(0,max_range,zpad);
for ii = 1:size(S,1)
    %S(ii,:) = S(ii,:).*R.^(3/2); %Optional: magnitude scale to range
end
S = dbv(S(:,1:size(v,2)/2));
m = max(max(S));
imagesc(R,time,S-m,[-80, 0]);
colorbar;
ylabel('time (s)');
xlabel('range (m)');
title('RTI with 2-pulse cancelor clutter rejection');
 
% %2 pulse mag only cancelor
% figure(30);
% clear v;
% for ii = 1:size(sif,1)-1
%     v1 = abs(ifft(sif(ii,:),zpad));
%     v2 = abs(ifft(sif(ii+1,:),zpad));
%     v(ii,:) = v2-v1;
% end
% S=v;
% R = linspace(0,max_range,zpad);
% for ii = 1:size(S,1)
%     S(ii,:) = S(ii,:).*R.^(3/2); %Optional: magnitude scale to range
% end
% S = dbv(S(:,1:size(v,2)/2));
% m = max(max(S));
% imagesc(R,time,S-m,[-20, 0]);
% colorbar;
% ylabel('time (s)');
% xlabel('range (m)');
% title('RTI with 2-pulse mag only cancelor clutter rejection');

sar:

SBAND_RMA_opendata.m

MIT IAP Radar Course 2011
%Resource: Build a Small Radar System Capable of Sensing Range, Doppler, 
%and Synthetic Aperture Radar Imaging 
%
%Gregory L. Charvat
 
%SAR algorithm from:
%Range Migration Algorithm from ch 10 of Spotlight Synthetic Aperture Radar
%Signal Processing Algorithms, Carrara, Goodman, and Majewski
 
%NOTE: set up-ramp sweep from 2-3.2V to stay within ISM band
%change fstart and fstop bellow when in ISM band
 
%-------------------------------------------%
%Process raw data here
clear all;
close all;
 
%read the raw data .wave file here
[Y,FS,NBITS] = wavread('towardswarehouse.wav');
 
%constants
c = 3E8; %(m/s) speed of light
 
%radar parameters
Tp = 20E-3; %(s) pulse time
Trp = 0.25; %(s) min range profile time duration
N = Tp*FS; %# of samples per pulse
fstart = 2260E6; %(Hz) LFM start frequency
fstop = 2590E6; %(Hz) LFM stop frequency
%fstart = 2402E6; %(Hz) LFM start frequency for ISM band
%fstop = 2495E6; %(Hz) LFM stop frequency for ISM band
BW = fstop-fstart; %(Hz) transmti bandwidth
f = linspace(fstart, fstop, N/2); %instantaneous transmit frequency
 
%the input appears to be inverted
trig = -1*Y(:,1);
s = -1*Y(:,2);
clear Y;
 
%parse data here by position (silence between recorded data)
rpstart = abs(trig)>mean(abs(trig));
count = 0;
Nrp = Trp*FS; %min # samples between range profiles
 
for ii = Nrp+1:size(rpstart,1)-Nrp
    if rpstart(ii) == 1 && sum(rpstart(ii-Nrp:ii-1)) == 0
        count = count + 1;
        RP(count,:) = s(ii:ii+Nrp-1);
        RPtrig(count,:) = trig(ii:ii+Nrp-1);
    end
end
 
%parse data by pulse
count = 0;
thresh = 0.08;
clear ii;
for jj = 1:size(RP,1)
    %clear SIF;
    SIF = zeros(N,1);
    start = (RPtrig(jj,:)> thresh);
    count = 0;
    jj
    for ii = 12:(size(start,2)-2*N)
        [Y I] =  max(RPtrig(jj,ii:ii+2*N));
        if mean(start(ii-10:ii-2)) == 0 && I == 1
            count = count + 1;
            SIF = RP(jj,ii:ii+N-1)' + SIF;
        end
    end
    %hilbert transform
    q = ifft(SIF/count);
    sif(jj,:) = fft(q(size(q,1)/2+1:size(q,1)));
end
sif(find(isnan(sif))) = 1E-30; %set all Nan values to 0
%SAR data should be ready here
clear s;
s = sif;
save routsidewarehouse2 s; %for image data
%-------------------------------------------%
%load additional varaibles and setup constants for radar here
clear all;
c = 3E8; %(m/s) speed of light
%load IQ converted data here
load routsidewarehouse2 s; %load variable sif %for image data
for ii = 1:size(s,1)
    s(ii,:) = s(ii,:) - mean(s,1);
end
%sif = s-sif_sub; %perform coherent background subtraction
%sif = sif_sub; %image just the background
sif = s; %image without background subtraction
clear s;
clear sif_sub;
%***********************************************************************
%radar parameters
fc = (2590E6 - 2260E6)/2 + 2260E6; %(Hz) center radar frequency
B = (2590E6 - 2260E6); %(hz) bandwidth
cr = B/20E-3; %(Hz/sec) chirp rate
Tp = 20E-3; %(sec) pulse width
%VERY IMPORTANT, change Rs to distance to cal target
%Rs = (12+9/12)*.3048; %(m) y coordinate to scene center (down range), make this value equal to distance to cal target
Rs = 0;
Xa = 0; %(m) beginning of new aperture length
delta_x = 2*(1/12)*0.3048; %(m) 2 inch antenna spacing
L = delta_x*(size(sif,1)); %(m) aperture length
Xa = linspace(-L/2, L/2, (L/delta_x)); %(m) cross range position of radar on aperture L
Za = 0;
Ya = Rs; %THIS IS VERY IMPORTANT, SEE GEOMETRY FIGURE 10.6
t = linspace(0, Tp, size(sif,2)); %(s) fast time, CHECK SAMPLE RATE
Kr = linspace(((4*pi/c)*(fc - B/2)), ((4*pi/c)*(fc + B/2)), (size(t,2)));
%Save background subtracted and callibrated data
save sif sif delta_x Rs Kr Xa;
%clear all;
%run IFP
SBAND_RMA_IFP;