时域脉冲压缩和频域脉冲压缩比较附matlab源码.doc

首先在脉冲压缩之前要搞清楚以下两点，而不是上来就自己想着设计一个匹配滤波器进行脉冲压缩，这样往往会导致错误的结果!
1. LFM波形是up-chirp的，还是down-chirp的？
2. LFM的带宽B，脉冲持续时间Te（不是Tr噢），mu = B/Te和TS(时域采样间隔)
其次，要清楚频域脉冲压缩的实质是快速计算时域脉冲压缩，这两个是一回事！
1 时域脉冲压缩，采用xcorr(等长)，或conv(不等长)
nscat = 3;
rrec = 200;
b = 100e6;
scat_range = [10, 30, 100];
scat_rcs = [1 1 2];
taup = 0.005e-3;
f0 = 5.6e9; % carrier freq ,由于exp(-i*2*pi*f0*2*range/c)与t无关，此项是个固定的相移，不起作用，建模中可以忽略
c = 3.e8; % speed of light
fs = 2*b; % sampling freq
sampling_interval = 1/fs; 
n = fix(taup/sampling_interval);
freqlimit = 0.5*fs;
freq = linspace(-freqlimit,freqlimit,n); % notice freq_sampling_interval = fs/n = 1/taup;
t = linspace(-taup/2,taup/2,n);
% initialize input, output and replica vectors
x(nscat,1:n) = 0.; % x is a nscat-by-n matrix
y(1:n) = 0.;
replica(1:n) = 0.;
% baseband lfm signal
replica = exp(i * pi * (b/taup) .* t.^2);
% 我以前不知道原来接收的一个脉冲中包含了多个目标的回波
for j = 1:1nscat
    range = scat_range(j) ;
%对于函数Y来说，平移都是加．．
%其实，平移都是针对函数Y来说的，而不是x．．．
%因为x是y的自变量,分析的时候将对函数Y的加，等效的变成了自变量的减！
%MATLAB针对自变量X的移动 (t - tau) 必须变成针对函数Y本身的移动 (t + tau)
% t +(2*range/c) means target locate at tau = 2*range/c
%  t +(2*range/c) MUST greater than t, 不是课本上自变量的移动 t - (2*range/c) ！
% f(t-t0)是自变量t的范围不动,比如观测范围永远是t=0:100ms，用于不会有t0=200ms的图像出现在f(t-t0)的图像中
% 所以f(t-t0)是通过调整自变量实现曲线的移动的，而且自变量t的范围固定，这样移动曲线y不具有物理意义。
% 要采用直接移动y的方式移动曲线，要t的范围不固定，比如将t=0:100ms的回波波形，移动到t=100:200，移动后要画图的话，t的范围要变大。
% 也就是说t凡是书上针对自变量 y = f(x - x0)移动曲线的写法在Matlab中建模雷达回波时要写成 y = f(x + x0)，才有物理意义，注意x的范围已经变大

    x(j,:) = scat_rcs(j) .*exp(-i*2*pi*f0*2*range/c).* exp(i * pi * (b/taup) .* (t +(2*range/c)).^2) ; 
    y = x(j,:)  + y; 
end
out =xcorr(replica, y); % 注意：xcorr(replica, y)中 replica和y的次序不能颠倒！
out = out ./ n; % 除以n是有偏自相关估计的要求
s = taup * c /2; %最大不模糊距离
Npoints = ceil(rrec * n /s); % 最大不模糊距离s上采了n个点，问在接收窗rrec大小的距离上采了多少点？是个比例式
% n是LFM信号taup上的采样点数，按说 Ru = c*taup/2 = s上面应该对应 n/2点才对啊，为什么脉冲压缩后这里s对应了n点呢？
% 后面做了个实验（见最后）验证了下，人家这里是对的我自己的理解是错的。先死记住这个结论：xcorr脉冲压缩后的n点对应最大不模糊距离(0~Te/2)！
% 原因是经过脉冲压缩后原来长度为n的序列变成了长度为2n-1的序列，前边的n-1点加上0点这n个点对应着(-Te/2~0)这部分结果，也是个最大不模糊距离，但是此部分由于对应负的延迟没有实际意义。而脉冲宽度为Te的LFM真正有用的是(0~Te/2)这部分结果，该部分进行xcorr或者卷积的匹配滤波后点数从n/2变为n。 阿弥陀佛
dist =linspace(0, rrec, Npoints);
figure(3)
plot(dist,abs(out(n:n+Npoints-1)),'k'); % 用xcorr匹配滤波有n-1点的延迟，看看1:n-1点被跳过
xlabel ('Target relative position in meters')
ylabel ('Compressed echo')
grid
figure; plot(abs(out(1000:end)))  % 跳过1000-1点out(1000:1999)共1000点
下图是时域脉冲压缩结果，标准答案！

2 频域脉冲压缩，用FFT以及IFFT
% 注意频域脉冲压缩要求，y 和 replica等长度，否则结果和时域脉冲压缩对不上，有点mystery！这里以前困扰过我，其实本质是数字信号处理的频域采样定理没有搞懂，不懂的请看解释。
% 这一点要特别注意，原因是若FFT的长度为L>=M+N-1，则频域FFT计算的循环卷积的结果以L延拓，若循环卷积的长度取M点，则线性卷积的结果有L-M点是混叠的，只有前面的M点是有用的数据，即正确的M点线性卷积结果。
1） 由out结果可以看出，设replica的长度为L，则时域脉冲压缩的结果有(2L-1)那么长，若不进行2L-1点的FFT而采用L点的FFT则频域脉冲压缩的结果是先将时域脉冲压缩的结果按L为周期左右延拓，由此可见主值序列延拓后两边的两个周期各有L-1点的混叠到主值序列中，然后取主值序列得到最终的频域脉冲压缩结果。频域脉冲压缩结果的第一个点对应着时域脉冲压缩结果的第L个点，去掉了时域脉冲压缩的前L-1点！这不正说明了频域脉冲压缩结果是时域脉冲压缩结果的以L为周期延拓后的主值序列。
 2）而时域脉冲压缩中的前L-1个点恰好是长度为L的匹配滤波器没有完全滑动到数据y时的暂态过程，将其去掉刚刚好。
所以这里有点Magic！以下验证之
out1 out';
out2 = out1 + circshift(out1, 1000)+ circshift(out1, -1000);
plot(abs(out2))

nscat = 3;
rrec = 200;
b = 100e6;
scat_range = [10, 30, 100];
scat_rcs = [1 1 2];
taup = 0.005e-3;
f0 = 5.6e9; % carrier freq ,由于exp(-i*2*pi*f0*2*range/c)与t无关，此项是个固定的相移，不起作用，建模中可以忽略
c = 3.e8; % speed of light
fs = 2*b; % sampling freq
sampling_interval = 1/fs; 
n = fix(taup/sampling_interval);
nfft =n % 不能等于这个 2.^ceil(log2(n))，会造出误差的
freqlimit = 0.5*fs;
freq = linspace(-freqlimit,freqlimit,n); % notice freq_sampling_interval = fs/n = 1/taup;
t = linspace(-taup/2,taup/2,n);
% initialize input, output and replica vectors
x(nscat,1:n) = 0.; % x is a nscat-by-n matrix
y(1:n) = 0.;
replica(1:n) = 0.;
% baseband lfm signal
replica = exp(i * pi * (b/taup) .* t.^2);
% 我以前不知道原来接收的一个脉冲中包含了多个目标的回波
for j = 1:1:nscat
    range = scat_range(j) ;
%对于函数Y来说，平移都是加．．
%其实，平移都是针对函数Y来说的，而不是x．．．
%因为x是y的自变量,分析的时候将对函数Y的加，等效的变成了自变量的减！
%MATLAB针对自变量X的移动 (t - tau) 必须变成针对函数Y本身的移动 (t + tau)
% t +(2*range/c) means target locate at tau = 2*range/c
%  t +(2*range/c) MUST greater than t, 不是课本上自变量的移动 t - (2*range/c) ！
% f(t-t0)是自变量t的范围不动,比如观测范围永远是t=0:100ms，用于不会有t0=200ms的图像出现在f(t-t0)的图像中
% 所以f(t-t0)是通过调整自变量实现曲线的移动的，而且自变量t的范围固定，这样移动曲线y不具有物理意义。
% 要采用直接移动y的方式移动曲线，要t的范围不固定，比如将t=0:100ms的回波波形，移动到t=100:200，移动后要画图的话，t的范围要变大。
% 也就是说t凡是书上针对自变量 y = f(x - x0)移动曲线的写法在Matlab中建模雷达回波时要写成 y = f(x + x0)，才有物理意义，注意x的范围已经变大
    x(j,:) = scat_rcs(j) .*exp(-i*2*pi*f0*2*range/c).* exp(i * pi * (b/taup) .* (t +(2*range/c)).^2) ; 
    y = x(j,:)  + y; 
end
deltat = taup / nfft;
rfft = fft(replica,nfft);
yfft = fft(y,nfft);
out= abs(ifft((rfft .* conj(yfft)))) ./ (nfft);
s = taup * c /2; Npoints = ceil(rrec * nfft /s);
dist =linspace(0, rrec, Npoints);
plot(dist, out(1:Npoints))
xlabel ('Target relative position in meters')
ylabel ('Compressed echo')
grid
 plot(abs(out)) 可以看出结果与时域脉冲压缩结果是一致的！

3 频域脉冲压缩，用FFT以及IFFT时h(n)和x(n)的长度不相等时怎么办？
1）要补零计算循环卷积
% circular convolution length L >= M + N - 1
nLength = 2.^(ceil(log2(M + length(sig_lfm) -1)));
2）截取与时域卷积相等的部分
mxFreqPulsCmprs = mxFreqPulsCmprs(length(sig_lfm):length(sig_lfm) + M - 1,:);

----------------------------------------------------------------------------------------------------------

rrec = 2.4000e+004; % taup * c /2
b = 0.8e6;
taup = 160E-6;
R0 = 20e3;
f0 = 1e10; 
c = 3.e8; % speed of light
sampling_interval = 1e-6; 
n = fix(taup/sampling_interval);
t = linspace(-taup/2,taup/2,n);
% baseband lfm signal
replica = exp(i * pi * (b/taup) .* t.^2);

y = exp(-i*2*pi*f0*2*R0/c).* exp(i * pi * (b/taup) .* (t +(2*R0/c)).^2) ;
out =xcorr(replica, y);% 注意：xcorr(replica, y)中 replica和y的次序不能颠倒！
out = out ./ n; % 除以n是有偏自相关估计的要求
s = taup * c /2; %最大不模糊距离
Npoints = ceil(rrec * n /s); % 最大不模糊距离s上采了n个点，问在接收窗rrec大小的距离上采了多少点？是个比例式
dist =linspace(0, rrec, Npoints);
figure(3)
plot(dist,abs(out(n:n+Npoints-1)),'k'); % 用xcorr匹配滤波有n-1点的延迟，看看1:n-1点被跳过
xlabel ('Target relative position in meters')
ylabel ('Compressed echo')
grid