R语言的数字信号处理

minolink

于 2022-12-22 04:52:17 发布

阅读量550

点赞数

文章标签：信号处理 r语言

本文链接：https://blog.csdn.net/minolink/article/details/128403888

版权

本文介绍了如何使用R语言处理FMCW雷达的二进制浮点数据，包括读取二进制文件、实现汉明窗函数及进行FFT变换。通过R的readBin函数读取数据，利用apply函数进行矩阵运算，并展示了R语言在数据分析和图形绘制方面的高效性。虽然Python和MATLAB也是常用工具，但作者认为R在特定场景下更具优势。

摘要由CSDN通过智能技术生成

为何用R

对于算法原型来讲，实际上MATLAB才是数字信号甚至包括CV数据处理的最好工具，不过MATLAB实在太厚重了，而且从正规使用来讲成本相当的高，公司大概是不会拨出这个预算的，个人更不可能了。考虑到一些邮件交互的时候还是要尽量避免版权问题，需要一个免费的数据处理和图形工具，所以R自然就是首选了。

实际上R的官方包一直都在本子里装着，只是平时大都拿来当计算器用（相当好用，秒算极限求导微分积分，用来辅导高中生简直神一般的存在）。不过R的命令行还是不方便做工程化的事情的，所以装了个人版RStudio，方便一边写函数脚本一边调试，还能像MATLAB一样随时看工作变量和图形。

对了，为什么不用Python，因为Python比MATLAB还厚重，哐哐的吃内存，实时绘图也慢一些。某种意义上Python还是更适合“计算”而不是“研究”。

数据说明

用的是一个3天线FMCW雷达的ADC信号，每个chirp采样128次。数据存储是二进制浮点格式：

ANT_0_float[128] | ANT_1_float[128] | ANT_2_float[128] | ...

即每一帧依次保存每个天线的128个浮点数。

数据可以从这里下载。

https://pan.baidu.com/s/1v3bEQQRT0TPK1uUKJzEOQQ?pwd=qkfk

R语言处理

语法建议去翻官方的Reference pdf链接，很全面。

读二进制文件

比较不方便的是数组没有那么灵活，不能直接写m[...] = ...，需要提前定义。不过考虑到效率问题实际上MATLAB也会疯狂提示你预先定义类型并预分配空间，只是R这里就变成强制了。

readframe <- function(filename){
  fid = file(filename, "rb")
  nrx = 3
  samples_per_chirp = 128
  frameno = 1
  m = array(dim=c(nrx,10000,samples_per_chirp))
    repeat{
    rx1 = readBin(fid, single(), n=samples_per_chirp)
    rx2 = readBin(fid, single(), n=samples_per_chirp)
    rx3 = readBin(fid, single(), n=samples_per_chirp)
    
    if(length(rx3)!=samples_per_chirp){
      break
    }
    
    m[1,frameno,] = rx1
    m[2,frameno,] = rx2
    m[3,frameno,] = rx3
    frameno = frameno+1
  }
  close(fid)
  cat("frameno ", frameno, "\n")
  m[,1:frameno,]
}

窗函数和FFT

基础包是没有窗函数的，不过自己写一个也无妨。


hann <- function(N){
  n = c(0:N-1)
  w = 0.5 - 0.5*cos(2*pi*n/(N-1))
  w
}

fft_hann <- function(x){
  samples = length(x)
  w = hann(samples-1)
  xw = x*w
  y = fft(xw)
  y
}

所有数据做Range-FFT

这里用到R语言的神器apply来调用上边的带窗fft，这里要注意的是apply后的处理结果是保存在结果数组的第1维的。也就是说c(nRx, frames, samples)在apply处理之后会变成c(fft_N,nRx, frames)。

R的apply如此好用以至于我自己写的一个c++的BLAS库也实现了一个类似的调用，相当方便。你去翻一翻MATLAB的工具箱就会发现大量重复的代码在判断输入输出的数据维度和合法性，太low了。


rangefft <- function(m){
  dims = dim(m)
  nrx = dims[1]
  frames = dims[2]
  samples = dims[3]
  n = samples
  iq = apply(m, MARGIN=c(1,2), fft_hann)
  iq
}

测试一下

我把上边的函数都写在radarfile.R，记得每次改完.R脚本，需要重新source一下，因为source(...)实际上不是“引用”而是“执行”。

> filename = 不给你们看
> source("radarfile.R")
> adc = readframe(filename)
> iq = rangefft(adc)
> iq_mean_rx = apply(iq, c(1,3), mean)
> dim(iq_mean_rx)
[1]  128 2751
> contour(abs(iq_mean_rx))

画出来是下边这样，这个是一个水流测算的雷达数据，因为ADC是实数信号，所以距离FFT做出来是对称的，实际只取一半。横坐标是距离，纵坐标是时间，还是可以看出来水流的发生和结束以及水面的变化。