上一篇帖子,看完了基于SDR的多普勒雷达,就可以看看硬件雷达的多普勒测速的DSP代码了。
先看一下这个图:
我们需要的多普勒频移的测量结果是从混频器(Multiply Conjugate)输出的,也就是说我们只要能得到混频器的输出就能算出我们要的东西了。
再开这个图:
混频器硬件的输出经过了Video Amp1到达了音频的右声道,而音频左声道里是锯齿波信号源的输出。
所以只是音频右声道的数据对于多普勒测速是有用的(左声道数据暂时没用,但是对于另外两个程序可能会有用)。
了解到这些信息,再看程序会方便不少。
clear all;
close all;
%read the raw data .wave file here
[Y,FS,NBITS] = wavread('Off of Newton Exit 17.wav');
%constants
c = 3E8; %(m/s) speed of light
%radar parameters
Tp = 0.250; %(s) pulse time
N = Tp*FS; %# of samples per pulse
fc = 2590E6; %(Hz) Center frequency (connected VCO Vtune to +5 for example)
%fc = 2495E6; %(Hz) Center frequency within ISM band (VCO Vtune to +3.2V)
%the input appears to be inverted
s = -1*Y(:,2);
clear Y;
%creat doppler vs. time plot data set here
for ii = 1:round(size(s,1)/N)-1
sif(ii,:) = s(1+(ii-1)*N:ii*N);
end
%subtract the average DC term here
sif = sif - mean(s);
zpad = 8*N/2;
%doppler vs. time plot:
v = dbv(ifft(sif,zpad,2));
v = v(:,1:size(v,2)/2);
mmax = max(max(v));
%calculate velocity
delta_f = linspace(0, FS/2, size(v,2)); %(Hz)
lambda=c/fc;
velocity = delta_f*lambda/2;
%calculate time
time = linspace(1,Tp*size(v,1),size(v,1)); %(sec)
%plot
imagesc(velocity,time,v-mmax,[-35, 0]);
colorbar;
xlim([0 40]); %limit velocity axis
xlabel('Velocity (m/sec)');
ylabel('time (sec)');
这个代码一开始是读取wav数据,把耳机2个通道的数据都读出来放到数组Y中。然后初始化了一下常量,比如光速,还有载波频率等。后面就是从Y数组中读取右声道的数据(可以观察到Y(:,2)表示只读取矩阵第二列的数据),因位左声道只有ramp信号源的数据,右声道才是硬件混频器输出的数据。然后还通过减去平均值(mean函数)消除了dc offset。这些都不复杂,没什么特别的,英文注释可以明白。
真正的运算在doppler vs time plot这里
%doppler vs. time plot:
v = dbv(ifft(sif,zpad,2));
v = v(:,1:size(v,2)/2);
mmax = max(max(v));
这里做的是傅里叶变换,求出了混频信号的频率输出,对应了多普勒频移的大小,也对应了速度。
但是为啥函数用了ifft而不是fft呢?
我记得以前数学上证明过,fft和ifft虽然是一堆逆运算,但是其实过程是差不多的。由于本博客不是学术论文,我就不严肃证明了,可以看看下面这个文章,从程序的角度对比了两种运算。
在MatLab中FFT和IFFT的互相转换_一个做底层的码农的博客-CSDN博客_ifft matlab
其实IFFT的计算原理之一就是将频域(注意频域是复数)数据进行取共轭复数(虚部取反),然后再进行FFT变换,这样便将频域信号转换到时域。因为FFT变换的结果是复数,所以从频域进行FFT变换过来的结果也是复数,而此时只需取复数的实部,便是原时域信号(同时此复数的虚部已经比非常小了,就省略不要了)。
而MATLAB自带函数并没有说明IFFT的所以然,所以我进行对IFFT的验证。计算了MATLAB自带的IFFT函数和进行虚部取反再计算FFT。两种计算方式表明:结果相同。所以说IFFT的计算,其实跟FFT计算原理一样。
可以发现ifft和fft只是虚部不一样, 但是我这里的fft只是针对实数的声音信号的,本就没有虚部,所以就是完全一样的了。
再看看下面这个文章。
ifft(SLM,[],dim)中的dim表示bai维度du,1表示列,2表示行。第二个参数为[]表示点数与原矩阵每列元zhi素数。如果是 ifft(SLM,N,1) 则表示列ifft的点数为N。
了解到了ifft函数的用法,第一个参数就是要用来做运算的数据,第二个参数是多少个FFT点(其实就是类似于分辨率的概念),最后是针对矩阵中的具体某一组数据做运算(因为本来矩阵里除了声音强度信息在数组里,还有另一个数组对应的是采样时间,我们要做fft的是声音强度数组)。
既然这样,我觉得我可以直接把ifft替换为fft函数来做傅里叶变换。
%doppler vs. time plot:
%v = dbv(ifft(sif,zpad,2));
v = dbv(fft(sif,zpad,2));
v = v(:,1:size(v,2)/2);
mmax = max(max(v));
这是运算结果,可以看到与第一篇文字里原始的测速程序的输出一致。
这一段程序中的dbv是dbv.m程序里实现的,只是把数值转换为分贝,比较简单,不详细展开了。
下面还有一段运算也比较重要,是把多普勒频移转换为速度。
%calculate velocity
delta_f = linspace(0, FS/2, size(v,2)); %(Hz)
lambda=c/fc;
velocity = delta_f*lambda/2;
继续参考这篇搜狐博客
技术干货:用LimeSDR Mini制作一台软件定义多普勒雷达_搜狐汽车_搜狐网
根据公式多普勒测速的速度结果等于多普勒频移大小乘以波长除以2。
这里波长是实际的无线电频率,也就是载波频率。
matlab代码里也是这样写的velocity = delta_f * lamda/2
至于delta_f是从之前的v算出来的,v只是表示在fft的数组里的位置,要把它与[0,FS/2]数组对应的值找出来,才是真正的频移。
从这里也可以看出来,delta_f不会超过FS/2,而FS是声音文件的采样率,取决于声卡的采样率。如果声卡采样率不够高,就无法分辨快速运动的物体了。
另外还有一点很重要:
这个测速实验,第一条是把信号源调到CW模式,也就是说测速的时候硬件雷达用的也是固定频率信号源,而不是ramp,这样它发出的不是chirp信号而是和SDR的测速雷达信号一样。
为啥测距和SAR要改信号源,下一篇测距的时候会详细讲。