init
This commit is contained in:
commit
b9143a8db3
|
|
@ -0,0 +1,8 @@
|
|||
.vscode/
|
||||
vsipl/
|
||||
|
||||
radar
|
||||
|
||||
*.o
|
||||
*.a
|
||||
*.so
|
||||
|
|
@ -0,0 +1,31 @@
|
|||
TARGET = radar
|
||||
CC = cc
|
||||
PY = python3
|
||||
INC_DIR = -I./include -I./src -I../vsipl/include
|
||||
LIB_DIR = -L./vsipl/lib
|
||||
LIBS = -lvsip -lfftw3f -lm
|
||||
CFLAGS = $(INC_DIR) $(LIB_DIR) $(LIBS) -g
|
||||
|
||||
SRC = $(wildcard src/*.c)
|
||||
OBJ = $(patsubst src/%.c, obj/%.o, $(SRC))
|
||||
|
||||
$(TARGET): $(OBJ)
|
||||
@mkdir -p bin
|
||||
@$(CC) $(OBJ) $(CFLAGS) -o bin/$(TARGET)
|
||||
|
||||
obj/%.o: src/%.c include/*.h
|
||||
@mkdir -p obj
|
||||
@$(CC) -c $< $(CFLAGS) -o $@
|
||||
|
||||
.PHONY: clean all run
|
||||
|
||||
all: $(TARGET)
|
||||
|
||||
run: $(TARGET)s
|
||||
@mkdir -p data
|
||||
@./$(TARGET)
|
||||
@$(PY) ./plot.py
|
||||
|
||||
clean:
|
||||
rm -f $(TARGET) $(OBJ)
|
||||
rm -rf data obj
|
||||
|
|
@ -0,0 +1,351 @@
|
|||
# 基于飞腾平台的远距离探测系统
|
||||
|
||||
## 项目内容
|
||||
|
||||
项目基于 VISPL 函数库实现以下内容
|
||||
|
||||
### 封装汉明窗（Hamming Window）
|
||||
|
||||
$$
|
||||
w(n)=a_0-(1-a_0)\cos\left(\frac{2\pi n}{N-1}\right), \quad 0\leq n\leq N-1
|
||||
$$
|
||||
|
||||
### 利用余弦函数构造一个雷达回波脉冲（实信号）
|
||||
|
||||
线性调频信号的相位可以表示为：
|
||||
|
||||
$$
|
||||
\varphi(t)=2\pi f_{\text{low}}t+\frac{\pi\text{BW}}{\tau}t^2
|
||||
$$
|
||||
|
||||
线性调频的复信号则可以用欧拉公式表示为：
|
||||
|
||||
$$
|
||||
s(t) = e^{i\varphi(t)}
|
||||
$$
|
||||
|
||||
由于雷达回波得到的是实信号，所以我们只需要取复信号的实部即可，即应用余弦函数
|
||||
|
||||
$$
|
||||
s_{\text{real}}(t) = \cos\left(2\pi f_{\text{low}}t+\frac{\pi\text{BW}}{\tau}t^2\right)
|
||||
$$
|
||||
|
||||
设两个物体相聚为 $d$，则两个物体的回波相差的时间为
|
||||
|
||||
$$
|
||||
\Delta t = \frac{2d}{c}
|
||||
$$
|
||||
|
||||
其中 $c$ 为光速， $d$ 为两个物体的距离。假设接收到第一个物体反射信号的时间为 $t_0$，则雷达接收到的两个物体的信号分别为
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
s_1(t) &= \cos\left(2\pi f_{\text{low}}(t-t_0)+\frac{\pi\text{BW}}{\tau}(t-t_0)^2\right)\\
|
||||
s_2(t) &= \cos\left(2\pi f_{\text{low}}(t-t_0-\Delta t)+\frac{\pi\text{BW}}{\tau}(t-t_0-\Delta t)^2\right)
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
叠加的信号为
|
||||
|
||||
$$
|
||||
s(t) = s_1(t)+s_2(t)
|
||||
$$
|
||||
|
||||
### 在回波上叠加使其信噪比为 0dB 的高斯白噪声
|
||||
|
||||
信噪比的定义如下（单位为分贝）：
|
||||
|
||||
$$
|
||||
\text{SNR} = 10\log_{10}\frac{P_{\text{signal}}}{P_{\text{noise}}}
|
||||
$$
|
||||
|
||||
对于离散的采样，信噪比可以表示为：
|
||||
|
||||
$$
|
||||
\text{SNR} = 10\log_{10}\frac{\sum_{i=0}^{N-1}x_i^2}{\sum_{i=0}^{N-1}n_i^2}
|
||||
$$
|
||||
|
||||
其中 $x_i$ 为信号，$n_i$ 为噪声，且振幅满足分布：
|
||||
|
||||
$$
|
||||
n \sim \mathcal{N}\left(0, \frac{\sum x_i^2}{10^{\left(\frac{\text{SNR}}{10}\right)}N} \right)
|
||||
$$
|
||||
|
||||
### 设计希尔伯特滤波器
|
||||
|
||||
定义符号函数：
|
||||
|
||||
$$
|
||||
\text{sgn}(x)=\begin{cases}
|
||||
1, & x>0\\
|
||||
0, & x=0\\
|
||||
-1, & x<0
|
||||
\end{cases}
|
||||
$$
|
||||
|
||||
对于时域信号 $x(t)$，设其频域表示为 $X(\Omega)$，则在频域下的希尔伯特变换可以表示为：
|
||||
|
||||
$$
|
||||
\widehat{X}(\Omega) = [-j\text{sgn}(\Omega)]X(\Omega)
|
||||
$$
|
||||
|
||||
假设希尔伯特变换之后信号的时域表示为 $\widehat{x}(t)$，则希尔伯特滤波结果的时域表示为：
|
||||
|
||||
$$
|
||||
s(t) = x(t)+j\widehat{x}(t)
|
||||
$$
|
||||
|
||||
其频域表示为：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
S(\Omega) &= X(\Omega)+j\widehat{X}(\Omega)\\
|
||||
&= X(\Omega)+j[-j\text{sgn}(\Omega)]X(\Omega)\\
|
||||
&= [1+\text{sgn}(\Omega)]X(\Omega)
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
则希尔伯特滤波器的频域表示为：
|
||||
|
||||
$$
|
||||
H(\Omega) = 1+\text{sgn}(\Omega)
|
||||
$$
|
||||
|
||||
表示为时域下的卷积运算：
|
||||
|
||||
$$
|
||||
s(t) = h(t) * x(t)
|
||||
$$
|
||||
|
||||
其中 $h(t)$ 为希尔伯特滤波器 $H(\Omega)$ 的时域表示，通过逆傅里叶变换可以得到：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
h(t) &= \mathcal{F}^{-1}\{H(\Omega)\}\\
|
||||
&= \mathcal{F}^{-1}\{1+\text{sgn}(\Omega)\}\\
|
||||
&= \delta(t)+\frac{1}{\pi t}
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
但是由于所采用的是离散时间，所以需要对 $H(\Omega)$ 做离散时间傅里叶逆变换，得到 $h(t)$ 的离散表示：
|
||||
|
||||
$$
|
||||
\begin{aligned}
|
||||
h(n) &= \frac{1}{2\pi}\int_{-\pi}^{\pi}[1+\text{sgn}(\omega)]e^{j\omega n}\mathrm{d}\omega \\
|
||||
& = \frac{1}{2\pi}\int_{0}^{\pi}2e^{j\omega n}\mathrm{d}\omega
|
||||
= \frac{1}{\pi}\int_{0}^{\pi}e^{j\omega n}\mathrm{d}\omega \\
|
||||
& = \frac{1}{\pi}\frac{e^{j\pi n} - 1}{jn} = \frac{\cos\pi n - 1}{j\pi n} = \frac{j(1-\cos\pi n)}{\pi n}\\
|
||||
& = \begin{cases}
|
||||
1, & n=0\\
|
||||
\frac{2j}{\pi n}, & n \text{ 为奇数} \\
|
||||
0, & n \text{ 为偶数}
|
||||
\end{cases}
|
||||
\end{aligned}
|
||||
$$
|
||||
|
||||
### 应用脉冲压缩
|
||||
|
||||
### 检测目标数量和间距
|
||||
|
||||
根据一定的阈值对脉冲压缩之后所得到的信号进行筛选，之后线性遍历信号中离散的点，依据对应的下标和采样率得到时间，从而得到目标的距离。首先根据采样率 $f_{\mathrm{sample}}$ 得到采样的间隔
|
||||
$$
|
||||
\Delta t_{\mathrm{sample}} = \frac{1}{f_{\mathrm{sample}}}
|
||||
$$
|
||||
假设检测到了两个相邻的峰值，且二者的下标相差 $n$，则二者的时间间隔为
|
||||
$$
|
||||
\Delta t = n\Delta t_{\mathrm{sample}}
|
||||
$$
|
||||
最后根据光速计算得到距离
|
||||
$$
|
||||
d = \frac{c\Delta t}{2}
|
||||
$$
|
||||
|
||||
## 流程图与关键接口
|
||||
|
||||
### 流程图
|
||||
|
||||
```mermaid
|
||||
flowchart
|
||||
A --> O
|
||||
O[参考信号] --> F
|
||||
subgraph 回波信号生成
|
||||
A[生成信号] --> B[叠加高斯白噪声]
|
||||
B --> C[得到实信号]
|
||||
end
|
||||
subgraph 希尔伯特滤波
|
||||
C --> D[希尔伯特滤波]
|
||||
D --> E[得到复信号]
|
||||
E --> L
|
||||
end
|
||||
subgraph 脉冲压缩
|
||||
F[参考信号的共轭与翻转] --> G[应用汉明窗]
|
||||
G --> K
|
||||
H[计算傅里叶变换长度] -- "补零" --> K[参考信号]
|
||||
H -- "补零" --> L[回波信号]
|
||||
K -- "傅里叶变换" --> J[频域相乘]
|
||||
L -- "傅里叶变换" --> J
|
||||
J --> M[脉冲压缩结果]
|
||||
end
|
||||
M --> N[检测目标及距离]
|
||||
P[对噪声信号单独采样变换仿真得到阈值] --> N
|
||||
```
|
||||
|
||||
### 关键接口
|
||||
|
||||
希尔伯特滤波
|
||||
|
||||
```c
|
||||
/*
|
||||
* 内部接口：希尔伯特滤波
|
||||
* 参数：p_vector_src -- 输入信号
|
||||
* n_filter_length -- 滤波器长度
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：对输入信号进行希尔伯特滤波
|
||||
*/
|
||||
void hilbert(vsip_vview_f *p_vector_src, vsip_scalar_i n_filter_length,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
```
|
||||
|
||||
汉明窗
|
||||
|
||||
```c
|
||||
/*
|
||||
* 内部接口：生成汉明窗
|
||||
* 参数：p_vector_dst -- 输出信号
|
||||
* 功能：根据输出信号的长度生成汉明窗
|
||||
*/
|
||||
void vcreate_hamming_f(vsip_vview_f *p_vector_dst);
|
||||
```
|
||||
|
||||
信号生成以及处理
|
||||
|
||||
```c
|
||||
/*
|
||||
* 内部接口：生成线性调频信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成线性调频信号（复信号）
|
||||
*/
|
||||
void generate_lfm_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成线性调频信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成线性调频信号（实信号）
|
||||
*/
|
||||
void generate_lfm_signal_real(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成雷达回波信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* f_disatance -- 两个物体之间的距离
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成两个有一定距离的物体反射叠加得到的雷达回波信号
|
||||
*/
|
||||
void generate_radar_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_scalar_f f_distance, vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成雷达回波信号
|
||||
* 参数：p_vector_signal -- 输入信号
|
||||
* f_snr -- 目标信号信噪比
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成可以叠加到原信号上的给定信噪比的高斯白噪声
|
||||
*/
|
||||
void generate_wgn_signal(vsip_vview_f *p_vector_signal, vsip_scalar_f f_snr,
|
||||
vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：脉冲压缩
|
||||
* 参数：p_vector_signal_src -- 输入信号
|
||||
* p_vector_signal_ref -- 参考信号
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：使用给定的参考信号对输入信号进行脉冲压缩
|
||||
*/
|
||||
void pulse_compress(vsip_cvview_f *p_vector_signal_src, vsip_cvview_f *p_vector_signal_ref,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
|
||||
/*
|
||||
* 内部接口：检测信号
|
||||
* 参数：p_vector_signal -- 脉冲压缩之后得到的信号
|
||||
* f_threshold -- 阈值
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：对脉冲压缩之后的信号进行进一步检测，依据阈值进行筛选
|
||||
*/
|
||||
void detect_signal(vsip_cvview_f *p_vector_signal, vsip_scalar_f f_threshold,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
```
|
||||
|
||||
用于输出和调试的函数
|
||||
|
||||
```c
|
||||
/*
|
||||
* 内部接口：输出实向量
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_file -- 输出文件
|
||||
* 功能：将实向量的数据输出到文件
|
||||
*/
|
||||
void vdump_f(vsip_vview_f *p_vector, FILE *p_file);
|
||||
|
||||
/*
|
||||
* 内部接口：输出复向量
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_file -- 输出文件
|
||||
* 功能：将复向量的数据输出到文件
|
||||
*/
|
||||
void cvdump_f(vsip_cvview_f *p_vector, FILE *p_file);
|
||||
|
||||
/*
|
||||
* 内部接口：实向量调试
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_name -- 输出文件名
|
||||
* 功能：将实向量的数据输出到指定文件名的文件
|
||||
*/
|
||||
void vdebug_f(vsip_vview_f *p_vector, char *p_name);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量调试
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_name -- 输出文件名
|
||||
* 功能：将复向量的数据输出到指定文件名的文件
|
||||
*/
|
||||
void cvdebug_f(vsip_cvview_f *p_vector, char *p_name);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量翻转
|
||||
* 参数：p_vector_src -- 输入向量
|
||||
* p_vector_dst -- 输出向量
|
||||
* 功能：将输入向量的数据翻转后输出到输出向量
|
||||
*/
|
||||
void cvflip_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量填充
|
||||
* 参数：p_vector_src -- 输入向量
|
||||
* p_vector_dst -- 输出向量
|
||||
* 功能：根据输出向量的长度对输入向量进行零填充得到输出
|
||||
*/
|
||||
void cvpad_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst);
|
||||
|
||||
```
|
||||
|
||||
## 关于本项目
|
||||
|
||||
本项目为 NKU 2023 暑期实习实训飞腾课程 VSIPL 大作业。
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
#ifndef RADAR_FILTER_H_
|
||||
#define RADAR_FILTER_H_
|
||||
|
||||
#include "vsip.h"
|
||||
|
||||
/*
|
||||
* 内部接口：希尔伯特滤波
|
||||
* 参数：p_vector_src -- 输入信号
|
||||
* n_filter_length -- 滤波器长度
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：对输入信号进行希尔伯特滤波
|
||||
*/
|
||||
void hilbert(vsip_vview_f *p_vector_src, vsip_scalar_i n_filter_length,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,13 @@
|
|||
#ifndef RADAR_HAMMING_H_
|
||||
#define RADAR_HAMMING_H_
|
||||
|
||||
#include "vsip.h"
|
||||
|
||||
/*
|
||||
* 内部接口：生成汉明窗
|
||||
* 参数：p_vector_dst -- 输出信号
|
||||
* 功能：根据输出信号的长度生成汉明窗
|
||||
*/
|
||||
void vcreate_hamming_f(vsip_vview_f *p_vector_dst);
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,76 @@
|
|||
#ifndef RADAR_SIGNAL_H_
|
||||
#define RADAR_SIGNAL_H_
|
||||
|
||||
#include "vsip.h"
|
||||
|
||||
/*
|
||||
* 内部接口：生成线性调频信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成线性调频信号（复信号）
|
||||
*/
|
||||
void generate_lfm_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成线性调频信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成线性调频信号（实信号）
|
||||
*/
|
||||
void generate_lfm_signal_real(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成雷达回波信号
|
||||
* 参数：f_tau -- 脉冲宽度
|
||||
* f_freq_sampling -- 采样频率
|
||||
* f_freq_low -- 起始频率
|
||||
* f_band_width -- 带宽
|
||||
* f_distance -- 两个物体之间的距离
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成两个有一定距离的物体反射叠加得到的雷达回波信号
|
||||
*/
|
||||
void generate_radar_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_scalar_f f_distance, vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：生成雷达回波信号
|
||||
* 参数：p_vector_signal -- 输入信号
|
||||
* f_snr -- 目标信号信噪比
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：生成可以叠加到原信号上的给定信噪比的高斯白噪声
|
||||
*/
|
||||
void generate_wgn_signal(vsip_vview_f *p_vector_signal, vsip_scalar_f f_snr,
|
||||
vsip_vview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：脉冲压缩
|
||||
* 参数：p_vector_signal_src -- 输入信号
|
||||
* p_vector_signal_ref -- 参考信号
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：使用给定的参考信号对输入信号进行脉冲压缩
|
||||
*/
|
||||
void pulse_compress(vsip_cvview_f *p_vector_signal_src, vsip_cvview_f *p_vector_signal_ref,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：检测信号
|
||||
* 参数：p_vector_signal -- 脉冲压缩之后得到的信号
|
||||
* f_threshold -- 阈值
|
||||
* p_vector_dst -- 输出信号
|
||||
* 功能：对脉冲压缩之后的信号进行进一步检测，依据阈值进行筛选
|
||||
*/
|
||||
void detect_signal(vsip_cvview_f *p_vector_signal, vsip_scalar_f f_threshold,
|
||||
vsip_cvview_f *p_vector_dst);
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
#ifndef RADAR_UTILS_H_
|
||||
#define RADAR_UTILS_H_
|
||||
|
||||
#include "vsip.h"
|
||||
|
||||
/*
|
||||
* 内部接口：输出实向量
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_file -- 输出文件
|
||||
* 功能：将实向量的数据输出到文件
|
||||
*/
|
||||
void vdump_f(vsip_vview_f *p_vector, FILE *p_file);
|
||||
|
||||
/*
|
||||
* 内部接口：输出复向量
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_file -- 输出文件
|
||||
* 功能：将复向量的数据输出到文件
|
||||
*/
|
||||
void cvdump_f(vsip_cvview_f *p_vector, FILE *p_file);
|
||||
|
||||
/*
|
||||
* 内部接口：实向量调试
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_name -- 输出文件名
|
||||
* 功能：将实向量的数据输出到指定文件名的文件
|
||||
*/
|
||||
void vdebug_f(vsip_vview_f *p_vector, char *p_name);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量调试
|
||||
* 参数：p_vector -- 输入向量
|
||||
* p_name -- 输出文件名
|
||||
* 功能：将复向量的数据输出到指定文件名的文件
|
||||
*/
|
||||
void cvdebug_f(vsip_cvview_f *p_vector, char *p_name);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量翻转
|
||||
* 参数：p_vector_src -- 输入向量
|
||||
* p_vector_dst -- 输出向量
|
||||
* 功能：将输入向量的数据翻转后输出到输出向量
|
||||
*/
|
||||
void cvflip_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst);
|
||||
|
||||
/*
|
||||
* 内部接口：复向量填充
|
||||
* 参数：p_vector_src -- 输入向量
|
||||
* p_vector_dst -- 输出向量
|
||||
* 功能：根据输出向量的长度对输入向量进行零填充得到输出
|
||||
*/
|
||||
void cvpad_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst);
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,57 @@
|
|||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
|
||||
radar_signal = list(map(float, open('./data/radar_signal.txt', 'r').readlines()))
|
||||
|
||||
radar_signal_filtered = open('./data/radar_signal_filtered.txt', 'r').readlines()
|
||||
radar_signal_filtered = map(lambda x: x.replace(' ', '').replace('\n', ''), radar_signal_filtered)
|
||||
radar_signal_filtered = list(map(complex, radar_signal_filtered))
|
||||
|
||||
radar_signal_compressed = open('./data/radar_signal_compressed.txt', 'r').readlines()
|
||||
radar_signal_compressed = map(lambda x: x.replace(' ', '').replace('\n', ''), radar_signal_compressed)
|
||||
radar_signal_compressed = list(map(complex, radar_signal_compressed))
|
||||
|
||||
radar_signal_reduced = open('./data/radar_signal_reduced.txt', 'r').readlines()
|
||||
radar_signal_reduced = map(lambda x: x.replace(' ', '').replace('\n', ''), radar_signal_reduced)
|
||||
radar_signal_reduced = list(map(complex, radar_signal_reduced))
|
||||
|
||||
hilbert_coefficients = open('./data/hilbert_coefficients.txt', 'r').readlines()
|
||||
hilbert_coefficients = map(lambda x: x.replace(' ', '').replace('\n', ''), hilbert_coefficients)
|
||||
hilbert_coefficients = list(map(complex, hilbert_coefficients))
|
||||
|
||||
reference_signal = open('./data/reference_signal.txt', 'r').readlines()
|
||||
reference_signal = map(lambda x: x.replace(' ', '').replace('\n', ''), reference_signal)
|
||||
reference_signal = list(map(complex, reference_signal))
|
||||
|
||||
plt.plot(radar_signal, label='radar-signal')
|
||||
plt.plot(np.abs(radar_signal_filtered), label='radar-signal-filtered')
|
||||
plt.legend()
|
||||
plt.savefig('./data/radar_signal_filtered.png')
|
||||
|
||||
plt.clf()
|
||||
|
||||
plt.plot(np.abs(radar_signal_compressed), label='radar-signal-compressed-abs')
|
||||
plt.legend()
|
||||
plt.savefig('./data/radar_signal_compressed.png')
|
||||
|
||||
plt.clf()
|
||||
|
||||
plt.plot(np.abs(radar_signal_reduced), label='radar-signal-reduced-abs')
|
||||
plt.legend()
|
||||
plt.savefig('./data/radar_signal_reduced.png')
|
||||
|
||||
plt.clf()
|
||||
|
||||
plt.plot(np.real(hilbert_coefficients), label='hilbert-coefficients-real', marker='o', mfc='w')
|
||||
plt.plot(np.imag(hilbert_coefficients), label='hilbert-coefficients-imag', marker='o', mfc='w')
|
||||
plt.legend()
|
||||
plt.savefig('./data/hilbert_coefficients.png')
|
||||
|
||||
plt.clf()
|
||||
|
||||
plt.plot(np.real(reference_signal), label='reference-signal-real')
|
||||
plt.legend()
|
||||
plt.savefig('./data/reference_signal.png')
|
||||
|
||||
plt.clf()
|
||||
|
|
@ -0,0 +1,99 @@
|
|||
#include "filter.h"
|
||||
#include "utils.h"
|
||||
|
||||
void hilbert(vsip_vview_f *p_vector_src, vsip_scalar_i n_filter_length, vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
vsip_length n_signal_length = vsip_vgetlength_f(p_vector_src);
|
||||
|
||||
if (p_vector_src == NULL)
|
||||
{
|
||||
fprintf(stderr, "hilbert: p_vector_src is NULL\n");
|
||||
return;
|
||||
}
|
||||
else if (vsip_cvgetlength_f(p_vector_dst) != n_signal_length)
|
||||
{
|
||||
fprintf(stderr, "hilbert: p_vector_dst length is not equal to p_vector_src\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_cvview_f *p_vector_hilbert_transformer = vsip_cvcreate_f(n_filter_length, VSIP_MEM_NONE);
|
||||
|
||||
// 正负半轴元素个数
|
||||
vsip_scalar_i n_half_filter_length = n_filter_length / 2;
|
||||
|
||||
for (vsip_length n_index = 0; n_index < n_filter_length; n_index++)
|
||||
{
|
||||
// 计算索引偏移
|
||||
vsip_scalar_i n_index_offset = (vsip_scalar_i)n_index - n_half_filter_length;
|
||||
// 默认 0
|
||||
vsip_cscalar_f c_value = vsip_cmplx_f(0.0f, 0.0f);
|
||||
|
||||
if (n_index_offset == 0)
|
||||
{
|
||||
// 0 处的冲激响应
|
||||
c_value = vsip_cmplx_f(1.0f, 0.0f);
|
||||
}
|
||||
else if (n_index_offset % 2 != 0)
|
||||
{
|
||||
// 2j / (pi * n)
|
||||
c_value = vsip_cmplx_f(0.0f, 2.0f / (VSIP_PI * n_index_offset));
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
vsip_cvput_f(p_vector_hilbert_transformer, n_index, c_value);
|
||||
}
|
||||
|
||||
cvdebug_f(p_vector_hilbert_transformer, "./data/hilbert_coefficients.txt");
|
||||
|
||||
// 得到滤波器系数的实部和虚部
|
||||
vsip_vview_f *p_vector_hilbert_transformer_real =
|
||||
vsip_vcreate_f(n_filter_length, VSIP_MEM_NONE);
|
||||
vsip_vview_f *p_vector_hilbert_transformer_imag =
|
||||
vsip_vcreate_f(n_filter_length, VSIP_MEM_NONE);
|
||||
vsip_vreal_f(p_vector_hilbert_transformer, p_vector_hilbert_transformer_real);
|
||||
vsip_vimag_f(p_vector_hilbert_transformer, p_vector_hilbert_transformer_imag);
|
||||
|
||||
// 用于卷积结果的实部和虚部
|
||||
vsip_vview_f *p_vector_cmplx_signal_real = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
vsip_vview_f *p_vector_cmplx_signal_imag = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
|
||||
// 当前 VSIPL 版本中没有复数卷积的操作，由于希尔伯特滤波器系数为实数或虚数而没有复合，
|
||||
// 所以可以分别对实部和虚部进行卷积操作。
|
||||
|
||||
// 实部卷积核
|
||||
vsip_conv1d_f *p_filter_real =
|
||||
vsip_conv1d_create_f(p_vector_hilbert_transformer_real, VSIP_NONSYM, n_signal_length, 1,
|
||||
VSIP_SUPPORT_SAME, 0, VSIP_ALG_TIME);
|
||||
// 虚部卷积核
|
||||
vsip_conv1d_f *p_filter_imag =
|
||||
vsip_conv1d_create_f(p_vector_hilbert_transformer_imag, VSIP_NONSYM, n_signal_length, 1,
|
||||
VSIP_SUPPORT_SAME, 0, VSIP_ALG_TIME);
|
||||
|
||||
// 实部卷积
|
||||
vsip_convolve1d_f(p_filter_real, p_vector_src, p_vector_cmplx_signal_real);
|
||||
// 虚部卷积
|
||||
vsip_convolve1d_f(p_filter_imag, p_vector_src, p_vector_cmplx_signal_imag);
|
||||
|
||||
// 合成复信号
|
||||
for (vsip_length n_index = 0; n_index < n_signal_length; n_index++)
|
||||
{
|
||||
vsip_cscalar_f c_value = vsip_cmplx_f(vsip_vget_f(p_vector_cmplx_signal_real, n_index),
|
||||
vsip_vget_f(p_vector_cmplx_signal_imag, n_index));
|
||||
vsip_cvput_f(p_vector_dst, n_index, c_value);
|
||||
}
|
||||
|
||||
// 释放
|
||||
vsip_conv1d_destroy_f(p_filter_real);
|
||||
vsip_conv1d_destroy_f(p_filter_imag);
|
||||
vsip_vdestroy_f(p_vector_cmplx_signal_imag);
|
||||
vsip_vdestroy_f(p_vector_cmplx_signal_real);
|
||||
vsip_vdestroy_f(p_vector_hilbert_transformer_imag);
|
||||
vsip_vdestroy_f(p_vector_hilbert_transformer_real);
|
||||
vsip_cvalldestroy_f(p_vector_hilbert_transformer);
|
||||
}
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
#include "hamming.h"
|
||||
|
||||
void vcreate_hamming_f(vsip_vview_f *p_vector_dst)
|
||||
{
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "vcreate_hamming_f: p_vector_dst is NULL\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_length = vsip_vgetlength_f(p_vector_dst);
|
||||
|
||||
// 汉明窗有关参数
|
||||
const vsip_scalar_f f_a0 = (float)0.53836;
|
||||
const vsip_scalar_f f_a1 = (float)(1.0f - f_a0);
|
||||
|
||||
// 生成汉明窗
|
||||
for (vsip_length n_index = 0; n_index < n_length; n_index++)
|
||||
{
|
||||
vsip_scalar_f f_element =
|
||||
f_a0 - f_a1 * (float)vsip_cos_f(2.0 * VSIP_PI * n_index / (n_length - 1));
|
||||
vsip_vput_f(p_vector_dst, n_index, f_element);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,143 @@
|
|||
#include "filter.h"
|
||||
#include "hamming.h"
|
||||
#include "signal.h"
|
||||
#include "utils.h"
|
||||
#include "vsip.h"
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
vsip_init((void *)0);
|
||||
|
||||
// 脉冲宽度
|
||||
const vsip_scalar_f f_tau = 7e-6;
|
||||
// 采样率
|
||||
const vsip_scalar_f f_freq_sampling = 20e6;
|
||||
// 下限频率
|
||||
const vsip_scalar_f f_freq_low = 222e6;
|
||||
// 带宽
|
||||
const vsip_scalar_f f_band_width = 6e6;
|
||||
// 目标距离
|
||||
const vsip_scalar_f f_distance = 600.0f;
|
||||
// LFM 信号长度
|
||||
const vsip_length n_signal_length = (vsip_length)(0.5f + f_tau * f_freq_sampling);
|
||||
// 延迟时间
|
||||
const vsip_scalar_f f_delay_time = 2.0f * f_distance / 3e8;
|
||||
// 总时间
|
||||
const vsip_scalar_f f_total_time = f_tau + f_delay_time;
|
||||
// 延迟信号长度（样本数量）
|
||||
const vsip_length n_delay_signal_length = (vsip_length)(0.5f + f_delay_time * f_freq_sampling);
|
||||
// 总信号长度（样本数量）
|
||||
const vsip_length n_total_signal_length = (vsip_length)(0.5f + f_total_time * f_freq_sampling);
|
||||
// 希尔伯特系数
|
||||
const vsip_length n_hilbert_length = 11;
|
||||
|
||||
// 生成雷达接收的实信号
|
||||
vsip_vview_f *p_vector_radar_signal = vsip_vcreate_f(n_total_signal_length, VSIP_MEM_NONE);
|
||||
generate_radar_signal(f_tau, f_freq_sampling, f_freq_low, f_band_width, 600.0f,
|
||||
p_vector_radar_signal);
|
||||
|
||||
// DEBUG
|
||||
vdebug_f(p_vector_radar_signal, "./data/radar_signal.txt");
|
||||
|
||||
// 希尔伯特滤波
|
||||
vsip_cvview_f *p_vector_radar_signal_filtered =
|
||||
vsip_cvcreate_f(n_total_signal_length, VSIP_MEM_NONE);
|
||||
hilbert(p_vector_radar_signal, 11, p_vector_radar_signal_filtered);
|
||||
|
||||
// DEBUG
|
||||
cvdebug_f(p_vector_radar_signal_filtered, "./data/radar_signal_filtered.txt");
|
||||
|
||||
// 匹配滤波参考信号
|
||||
vsip_cvview_f *p_vector_signal_ref = vsip_cvcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
generate_lfm_signal(f_tau, f_freq_sampling, f_freq_low, f_band_width, p_vector_signal_ref);
|
||||
|
||||
// 傅里叶变换长度
|
||||
vsip_length n_fft_len = 2 * n_total_signal_length - 1;
|
||||
n_fft_len = (vsip_length)ceil(0.5 + log2(n_fft_len));
|
||||
n_fft_len = (vsip_length)floor(0.5 + pow(2, n_fft_len));
|
||||
|
||||
// DEBUG
|
||||
printf("radar: n_fft_len = %ld\n", n_fft_len);
|
||||
|
||||
// 脉冲压缩
|
||||
vsip_cvview_f *p_vector_radar_signal_compressed = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
pulse_compress(p_vector_radar_signal_filtered, p_vector_signal_ref,
|
||||
p_vector_radar_signal_compressed);
|
||||
|
||||
// DEBUG
|
||||
cvdebug_f(p_vector_radar_signal_compressed, "./data/radar_signal_compressed.txt");
|
||||
|
||||
// 根据 30 阈值筛选
|
||||
vsip_cvview_f *p_vector_radar_signal_reduced = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
detect_signal(p_vector_radar_signal_compressed, 30.0f, p_vector_radar_signal_reduced);
|
||||
|
||||
// DEBUG
|
||||
cvdebug_f(p_vector_radar_signal_reduced, "./data/radar_signal_reduced.txt");
|
||||
|
||||
// 检测目标并且输出距离
|
||||
// 采样间隔
|
||||
vsip_scalar_f f_delta_time = 1.0f / f_freq_sampling;
|
||||
// 上一个峰值的时间
|
||||
vsip_scalar_f f_prev_time = 0.0f;
|
||||
// 在向量中的下标
|
||||
vsip_length n_index = 0;
|
||||
|
||||
while (n_index < n_fft_len - 1)
|
||||
{
|
||||
// 当前时间
|
||||
vsip_scalar_f f_curr_time = (vsip_scalar_f)n_index * f_delta_time;
|
||||
// 下一个时间
|
||||
vsip_scalar_f f_next_time = (vsip_scalar_f)(n_index + 1) * f_delta_time;
|
||||
// 当前幅值
|
||||
vsip_scalar_f f_curr_mag =
|
||||
vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, n_index));
|
||||
// 下一个幅值
|
||||
vsip_scalar_f f_next_mag =
|
||||
vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, n_index + 1));
|
||||
|
||||
if (f_curr_mag <= f_next_mag)
|
||||
{
|
||||
// 若仍然处于递增阶段，则继续向后搜索
|
||||
n_index++;
|
||||
}
|
||||
else
|
||||
{
|
||||
// 当前幅值大于前一个幅值，且大于后一个幅值，说明是一个峰值
|
||||
if (f_prev_time == 0.0f)
|
||||
{
|
||||
// 第一个目标，记录时间
|
||||
f_prev_time = f_curr_time;
|
||||
printf("detect: time = %fs\n", f_curr_time);
|
||||
}
|
||||
else
|
||||
{
|
||||
// 与上一个目标的距离
|
||||
vsip_scalar_f f_distance = 3e8 * (f_curr_time - f_prev_time) / 2.0f;
|
||||
// 记录时间
|
||||
f_prev_time = f_curr_time;
|
||||
printf("detect: time = %fs, distance = %fm\n", f_curr_time, f_distance);
|
||||
}
|
||||
while (f_curr_mag > f_next_mag)
|
||||
{
|
||||
// 在递减阶段，继续向后搜索
|
||||
n_index++;
|
||||
f_curr_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, n_index));
|
||||
f_next_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, n_index + 1));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 释放内存
|
||||
vsip_cvalldestroy_f(p_vector_radar_signal_reduced);
|
||||
vsip_cvalldestroy_f(p_vector_radar_signal_compressed);
|
||||
vsip_cvalldestroy_f(p_vector_signal_ref);
|
||||
vsip_valldestroy_f(p_vector_radar_signal);
|
||||
vsip_cvalldestroy_f(p_vector_radar_signal_filtered);
|
||||
|
||||
vsip_finalize((void *)0);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,353 @@
|
|||
#include "signal.h"
|
||||
#include "hamming.h"
|
||||
#include "utils.h"
|
||||
#include <time.h>
|
||||
|
||||
void generate_lfm_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
// 信号长度（样本数量）
|
||||
vsip_length n_signal_length = (vsip_length)(0.5f + f_tau * f_freq_sampling);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "generate_lfm_signal: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (vsip_cvgetlength_f(p_vector_dst) != n_signal_length)
|
||||
{
|
||||
fprintf(stderr,
|
||||
"generate_lfm_signal: p_vector_dst length is not equal to n_signal_length.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
// 初始化时间向量
|
||||
vsip_vview_f *p_vector_time = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
// 初始化中间变量
|
||||
vsip_vview_f *p_vector_temp_0 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
// 初始化中间变量
|
||||
vsip_vview_f *p_vector_temp_1 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
|
||||
// 生成时间向量
|
||||
vsip_vramp_f(0.0f, 1.0f / f_freq_sampling, p_vector_time);
|
||||
|
||||
vsip_svmul_f(2.0 * VSIP_PI * f_freq_low, p_vector_time, p_vector_temp_0);
|
||||
|
||||
vsip_vsq_f(p_vector_time, p_vector_temp_1);
|
||||
vsip_svmul_f(VSIP_PI * f_band_width / f_tau, p_vector_temp_1, p_vector_temp_1);
|
||||
|
||||
// 相加得到相位
|
||||
vsip_vadd_f(p_vector_temp_0, p_vector_temp_1, p_vector_temp_0);
|
||||
|
||||
// 应用欧拉函数
|
||||
vsip_veuler_f(p_vector_temp_0, p_vector_dst);
|
||||
|
||||
// 释放内存
|
||||
vsip_valldestroy_f(p_vector_time);
|
||||
vsip_valldestroy_f(p_vector_temp_0);
|
||||
vsip_valldestroy_f(p_vector_temp_1);
|
||||
}
|
||||
|
||||
void generate_lfm_signal_real(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_vview_f *p_vector_dst)
|
||||
{
|
||||
// 信号长度（样本数量）
|
||||
vsip_length n_signal_length = (vsip_length)(0.5f + f_tau * f_freq_sampling);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "generate_lfm_signal: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (vsip_vgetlength_f(p_vector_dst) != n_signal_length)
|
||||
{
|
||||
fprintf(stderr,
|
||||
"generate_lfm_signal: p_vector_dst length is not equal to n_signal_length.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
// 初始化时间向量
|
||||
vsip_vview_f *p_vector_time = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
// 初始化中间变量
|
||||
vsip_vview_f *p_vector_temp_0 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
// 初始化中间变量
|
||||
vsip_vview_f *p_vector_temp_1 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
|
||||
// 生成时间向量
|
||||
vsip_vramp_f(0.0f, 1.0f / f_freq_sampling, p_vector_time);
|
||||
|
||||
vsip_svmul_f(2.0 * VSIP_PI * f_freq_low, p_vector_time, p_vector_temp_0);
|
||||
|
||||
vsip_vsq_f(p_vector_time, p_vector_temp_1);
|
||||
vsip_svmul_f(VSIP_PI * f_band_width / f_tau, p_vector_temp_1, p_vector_temp_1);
|
||||
|
||||
// 相加得到相位
|
||||
vsip_vadd_f(p_vector_temp_0, p_vector_temp_1, p_vector_temp_0);
|
||||
|
||||
// 应用余弦函数
|
||||
vsip_vcos_f(p_vector_temp_0, p_vector_dst);
|
||||
|
||||
// 释放内存
|
||||
vsip_valldestroy_f(p_vector_time);
|
||||
vsip_valldestroy_f(p_vector_temp_0);
|
||||
vsip_valldestroy_f(p_vector_temp_1);
|
||||
}
|
||||
|
||||
void generate_radar_signal(vsip_scalar_f f_tau, vsip_scalar_f f_freq_sampling,
|
||||
vsip_scalar_f f_freq_low, vsip_scalar_f f_band_width,
|
||||
vsip_scalar_f f_distance, vsip_vview_f *p_vector_dst)
|
||||
{
|
||||
|
||||
// 延迟时间
|
||||
vsip_scalar_f f_delay_time = 2.0f * f_distance / 3e8;
|
||||
// 总时间
|
||||
vsip_scalar_f f_total_time = f_tau + f_delay_time;
|
||||
|
||||
// 延迟信号长度（样本数量）
|
||||
vsip_length n_delay_signal_length = (vsip_length)(0.5f + f_delay_time * f_freq_sampling);
|
||||
// 总信号长度（样本数量）
|
||||
vsip_length n_total_signal_length = (vsip_length)(0.5f + f_total_time * f_freq_sampling);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "generate_radar_signal: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (vsip_vgetlength_f(p_vector_dst) != n_total_signal_length)
|
||||
{
|
||||
fprintf(
|
||||
stderr,
|
||||
"generate_radar_signal: p_vector_dst length is not equal to n_total_signal_length.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
// LFM 信号长度
|
||||
vsip_length n_signal_length = (vsip_length)(0.5f + f_tau * f_freq_sampling);
|
||||
// 第一个目标返回的信号
|
||||
vsip_vview_f *p_vector_signal_0 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
// 第二个目标返回的信号
|
||||
vsip_vview_f *p_vector_signal_1 = vsip_vcreate_f(n_signal_length, VSIP_MEM_NONE);
|
||||
|
||||
generate_lfm_signal_real(f_tau, f_freq_sampling, f_freq_low, f_band_width, p_vector_signal_0);
|
||||
generate_lfm_signal_real(f_tau, f_freq_sampling, f_freq_low, f_band_width, p_vector_signal_1);
|
||||
|
||||
// 延长原始信号向量
|
||||
vsip_vview_f *p_vector_signal_0_extended = vsip_vcreate_f(n_total_signal_length, VSIP_MEM_NONE);
|
||||
// 延长原始信号向量
|
||||
vsip_vview_f *p_vector_signal_1_extended = vsip_vcreate_f(n_total_signal_length, VSIP_MEM_NONE);
|
||||
|
||||
// 零填充
|
||||
vsip_vfill_f(0.0f, p_vector_signal_0_extended);
|
||||
vsip_vfill_f(0.0f, p_vector_signal_1_extended);
|
||||
|
||||
// 拷贝数据
|
||||
for (vsip_length n_index = 0; n_index < n_signal_length; n_index++)
|
||||
{
|
||||
// 原始信号
|
||||
vsip_vput_f(p_vector_signal_0_extended, n_index, vsip_vget_f(p_vector_signal_0, n_index));
|
||||
// 延迟信号
|
||||
vsip_vput_f(p_vector_signal_1_extended, n_index + n_delay_signal_length,
|
||||
vsip_vget_f(p_vector_signal_1, n_index));
|
||||
}
|
||||
|
||||
// 叠加信号
|
||||
vsip_vadd_f(p_vector_signal_0_extended, p_vector_signal_1_extended, p_vector_dst);
|
||||
|
||||
// 高斯白噪声，SNR=0dB
|
||||
vsip_vview_f *p_vector_wgn_signal = vsip_vcreate_f(n_total_signal_length, VSIP_MEM_NONE);
|
||||
generate_wgn_signal(p_vector_dst, 0.0f, p_vector_wgn_signal);
|
||||
|
||||
// 叠加噪声
|
||||
vsip_vadd_f(p_vector_dst, p_vector_wgn_signal, p_vector_dst);
|
||||
|
||||
// 释放内存
|
||||
vsip_valldestroy_f(p_vector_signal_0);
|
||||
vsip_valldestroy_f(p_vector_signal_1);
|
||||
vsip_valldestroy_f(p_vector_signal_0_extended);
|
||||
vsip_valldestroy_f(p_vector_signal_1_extended);
|
||||
vsip_valldestroy_f(p_vector_wgn_signal);
|
||||
}
|
||||
|
||||
void generate_wgn_signal(vsip_vview_f *p_vector_signal, vsip_scalar_f f_snr,
|
||||
vsip_vview_f *p_vector_dst)
|
||||
{
|
||||
vsip_length n_length = vsip_vgetlength_f(p_vector_signal);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "generate_wgn_signal: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (vsip_vgetlength_f(p_vector_dst) != n_length)
|
||||
{
|
||||
fprintf(stderr, "generate_wgn_signal: p_vector_dst length is not equal to n_length.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
// 功率信噪比
|
||||
vsip_scalar_f f_snr_power = vsip_pow_f(10.0f, f_snr / 10.0f);
|
||||
// 噪声功率
|
||||
vsip_scalar_f f_noise_power = vsip_vsumsqval_f(p_vector_signal) / (n_length * f_snr_power);
|
||||
// 随机生成
|
||||
vsip_randstate *p_rand_state = vsip_randcreate(time(NULL), 1, 1, VSIP_PRNG);
|
||||
vsip_vrandn_f(p_rand_state, p_vector_dst);
|
||||
vsip_svmul_f(vsip_sqrt_f(f_noise_power), p_vector_dst, p_vector_dst);
|
||||
|
||||
// 释放内存
|
||||
vsip_randdestroy(p_rand_state);
|
||||
}
|
||||
|
||||
/// 脉冲压缩
|
||||
void pulse_compress(vsip_cvview_f *p_vector_signal_src, vsip_cvview_f *p_vector_signal_ref,
|
||||
vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
|
||||
if (p_vector_signal_src == NULL)
|
||||
{
|
||||
fprintf(stderr, "pulse_compress: p_vector_signal_src is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (p_vector_signal_ref == NULL)
|
||||
{
|
||||
fprintf(stderr, "pulse_compress: p_vector_signal_ref is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "pulse_compress: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_signal_src_length = vsip_cvgetlength_f(p_vector_signal_src);
|
||||
vsip_length n_signal_ref_length = vsip_cvgetlength_f(p_vector_signal_ref);
|
||||
|
||||
// 求参考信号共轭
|
||||
vsip_cvconj_f(p_vector_signal_ref, p_vector_signal_ref);
|
||||
// 翻转参考信号
|
||||
vsip_cvview_f *p_vector_signal_ref_flipped =
|
||||
vsip_cvcreate_f(n_signal_ref_length, VSIP_MEM_NONE);
|
||||
cvflip_f(p_vector_signal_ref, p_vector_signal_ref_flipped);
|
||||
// 加汉明窗
|
||||
vsip_vview_f *p_vector_window = vsip_vcreate_f(n_signal_ref_length, VSIP_MEM_NONE);
|
||||
vcreate_hamming_f(p_vector_window);
|
||||
|
||||
vsip_rcvmul_f(p_vector_window, p_vector_signal_ref_flipped, p_vector_signal_ref);
|
||||
|
||||
cvdebug_f(p_vector_signal_ref, "./data/reference_signal.txt");
|
||||
|
||||
// 傅里叶变换长度
|
||||
vsip_length n_fft_len = 2 * n_signal_src_length - 1;
|
||||
n_fft_len = (vsip_length)ceil(0.5 + log2(n_fft_len));
|
||||
n_fft_len = (vsip_length)floor(0.5 + pow(2, n_fft_len));
|
||||
|
||||
if (vsip_cvgetlength_f(p_vector_dst) != n_fft_len)
|
||||
{
|
||||
fprintf(stderr, "pulse_compress: p_vector_dst length is not equal to n_fft_len.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
// 补零向量
|
||||
vsip_cvview_f *p_vector_signal_src_padded = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
vsip_cvview_f *p_vector_signal_ref_padded = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
|
||||
// 补零
|
||||
cvpad_f(p_vector_signal_src, p_vector_signal_src_padded);
|
||||
cvpad_f(p_vector_signal_ref, p_vector_signal_ref_padded);
|
||||
|
||||
// 傅里叶变换
|
||||
vsip_fft_f *p_fft = vsip_ccfftop_create_f(n_fft_len, 1.0, VSIP_FFT_FWD, 1, 0);
|
||||
vsip_fft_f *p_ifft = vsip_ccfftop_create_f(n_fft_len, 1.0 / n_fft_len, VSIP_FFT_INV, 1, 0);
|
||||
vsip_cvview_f *p_vector_signal_src_fft = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
vsip_cvview_f *p_vector_signal_ref_fft = vsip_cvcreate_f(n_fft_len, VSIP_MEM_NONE);
|
||||
|
||||
vsip_ccfftop_f(p_fft, p_vector_signal_src_padded, p_vector_signal_src_fft);
|
||||
vsip_ccfftop_f(p_fft, p_vector_signal_ref_padded, p_vector_signal_ref_fft);
|
||||
|
||||
// 相乘
|
||||
vsip_cvmul_f(p_vector_signal_src_fft, p_vector_signal_ref_fft, p_vector_signal_src_fft);
|
||||
|
||||
// 傅里叶逆变换
|
||||
vsip_ccfftop_f(p_ifft, p_vector_signal_src_fft, p_vector_dst);
|
||||
|
||||
// 释放内存
|
||||
vsip_cvalldestroy_f(p_vector_signal_src_fft);
|
||||
vsip_cvalldestroy_f(p_vector_signal_ref_fft);
|
||||
vsip_fft_destroy_f(p_fft);
|
||||
vsip_fft_destroy_f(p_ifft);
|
||||
vsip_cvalldestroy_f(p_vector_signal_src_padded);
|
||||
vsip_cvalldestroy_f(p_vector_signal_ref_padded);
|
||||
vsip_cvalldestroy_f(p_vector_signal_ref_flipped);
|
||||
vsip_valldestroy_f(p_vector_window);
|
||||
}
|
||||
|
||||
void detect_signal(vsip_cvview_f *p_vector_signal, vsip_scalar_f f_threshold,
|
||||
vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
if (p_vector_signal == NULL)
|
||||
{
|
||||
fprintf(stderr, "detect_signal: p_vector_signal is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "detect_signal: p_vector_dst is NULL.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_signal_length = vsip_cvgetlength_f(p_vector_signal);
|
||||
|
||||
if (vsip_cvgetlength_f(p_vector_dst) != n_signal_length)
|
||||
{
|
||||
fprintf(stderr, "detect_signal: p_vector_dst length is not equal to n_signal_length.\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
for (vsip_length n_index = 0; n_index < n_signal_length; n_index++)
|
||||
{
|
||||
vsip_scalar_f f_abs = vsip_cmag_f(vsip_cvget_f(p_vector_signal, n_index));
|
||||
if (fabs(f_abs) < f_threshold && fabs(f_abs) < f_threshold)
|
||||
{
|
||||
// 低于阈值的信号置零
|
||||
vsip_cvput_f(p_vector_dst, n_index, vsip_cmplx_f(0.0f, 0.0f));
|
||||
}
|
||||
else
|
||||
{
|
||||
// 大于阈值的信号保留
|
||||
vsip_cvput_f(p_vector_dst, n_index, vsip_cvget_f(p_vector_signal, n_index));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,147 @@
|
|||
#include "utils.h"
|
||||
|
||||
void vdump_f(vsip_vview_f *p_vector, FILE *p_file)
|
||||
{
|
||||
if (p_vector == NULL)
|
||||
{
|
||||
fprintf(stderr, "vdump_f: p_vector is NULL\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
vsip_length n_length = vsip_vgetlength_f(p_vector);
|
||||
for (vsip_index n_index = 0; n_index < n_length; n_index++)
|
||||
{
|
||||
fprintf(p_file, "%f\n", vsip_vget_f(p_vector, n_index));
|
||||
}
|
||||
}
|
||||
|
||||
void cvdump_f(vsip_cvview_f *p_vector, FILE *p_file)
|
||||
{
|
||||
if (p_vector == NULL)
|
||||
{
|
||||
fprintf(stderr, "vdump_f: p_vector is NULL\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
vsip_length n_length = vsip_cvgetlength_f(p_vector);
|
||||
for (vsip_index n_index = 0; n_index < n_length; n_index++)
|
||||
{
|
||||
vsip_cscalar_f c_value = vsip_cvget_f(p_vector, n_index);
|
||||
fprintf(p_file, "%+9.4f%+9.4fj\n", vsip_real_f(c_value), vsip_imag_f(c_value));
|
||||
}
|
||||
}
|
||||
|
||||
void vdebug_f(vsip_vview_f *p_vector, char *p_name)
|
||||
{
|
||||
FILE *p_file = fopen(p_name, "w");
|
||||
if (p_file == NULL)
|
||||
{
|
||||
fprintf(stderr, "vdebug_f: unable to open file `%s`\n", p_name);
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
vdump_f(p_vector, p_file);
|
||||
fclose(p_file);
|
||||
}
|
||||
|
||||
void cvdebug_f(vsip_cvview_f *p_vector, char *p_name)
|
||||
{
|
||||
FILE *p_file = fopen(p_name, "w");
|
||||
if (p_file == NULL)
|
||||
{
|
||||
fprintf(stderr, "cvdebug_f: unable to open file `%s`\n", p_name);
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
cvdump_f(p_vector, p_file);
|
||||
fclose(p_file);
|
||||
}
|
||||
|
||||
void cvflip_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
if (p_vector_src == NULL)
|
||||
{
|
||||
fprintf(stderr, "cvflip_f: p_vector is NULL\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_length = vsip_cvgetlength_f(p_vector_src);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "cvflip_f: p_vector_dst is NULL\n");
|
||||
return;
|
||||
}
|
||||
else if (n_length != vsip_cvgetlength_f(p_vector_dst))
|
||||
{
|
||||
fprintf(stderr, "cvflip_f: p_vector_src and p_vector_dst have different lengths\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
for (vsip_index n_index = 0; n_index < n_length; n_index++)
|
||||
{
|
||||
vsip_cvput_f(p_vector_dst, n_index, vsip_cvget_f(p_vector_src, n_length - n_index - 1));
|
||||
}
|
||||
}
|
||||
|
||||
void cvpad_f(vsip_cvview_f *p_vector_src, vsip_cvview_f *p_vector_dst)
|
||||
{
|
||||
if (p_vector_src == NULL)
|
||||
{
|
||||
fprintf(stderr, "cvpad_f: p_vector is NULL\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_length_src = vsip_cvgetlength_f(p_vector_src);
|
||||
|
||||
if (p_vector_dst == NULL)
|
||||
{
|
||||
fprintf(stderr, "cvpad_f: p_vector_dst is NULL\n");
|
||||
return;
|
||||
}
|
||||
else if (n_length_src > vsip_cvgetlength_f(p_vector_dst))
|
||||
{
|
||||
fprintf(stderr, "cvpad_f: p_vector_src is longer than p_vector_dst\n");
|
||||
return;
|
||||
}
|
||||
else
|
||||
{
|
||||
;
|
||||
}
|
||||
|
||||
vsip_length n_length_dst = vsip_cvgetlength_f(p_vector_dst);
|
||||
|
||||
for (vsip_index n_index = 0; n_index < n_length_src; n_index++)
|
||||
{
|
||||
vsip_cvput_f(p_vector_dst, n_index, vsip_cvget_f(p_vector_src, n_index));
|
||||
}
|
||||
|
||||
for (vsip_index n_index = n_length_src; n_index < n_length_dst; n_index++)
|
||||
{
|
||||
vsip_cvput_f(p_vector_dst, n_index, vsip_cmplx_f(0.0f, 0.0f));
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,99 @@
|
|||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def wgn(signal, snr = 0):
|
||||
snr = 10 ** (snr / 10)
|
||||
xpower = np.sum(signal ** 2) / len(signal)
|
||||
npower = xpower / snr
|
||||
noise = np.random.randn(len(signal)) * np.sqrt(npower)
|
||||
return noise
|
||||
|
||||
def hilbert(signal, num_hilbert = 11):
|
||||
hilbert_transformer = np.zeros(num_hilbert, dtype=np.complex128)
|
||||
for i in range(-num_hilbert // 2, num_hilbert // 2 + 1, 1):
|
||||
if i % 2 == 0:
|
||||
hilbert_transformer[i + num_hilbert // 2] = 0
|
||||
else:
|
||||
hilbert_transformer[i + num_hilbert // 2] = 2j / (np.pi * i)
|
||||
|
||||
hilbert_transformer[0 + num_hilbert // 2] = 1
|
||||
|
||||
signal_hilbert = np.convolve(signal, hilbert_transformer, mode='same')
|
||||
|
||||
return signal_hilbert
|
||||
|
||||
def pulse_compress(signal_src, signal_ref):
|
||||
signal_src_length = len(signal_src)
|
||||
signal_ref_length = len(signal_ref)
|
||||
|
||||
|
||||
fft_len = 2 * signal_src_length - 1
|
||||
fft_len = np.ceil(0.5 + np.log2(fft_len))
|
||||
fft_len = np.floor(0.5 + 2 ** fft_len)
|
||||
|
||||
signal_ref = np.conjugate(signal_ref)
|
||||
signal_ref = np.flip(signal_ref)
|
||||
|
||||
window = np.hamming(signal_ref_length)
|
||||
# signal_ref = signal_ref * window
|
||||
|
||||
signal_src = np.pad(signal_src, (0, int(fft_len - signal_src_length)), 'constant')
|
||||
signal_ref = np.pad(signal_ref, (0, int(fft_len - signal_ref_length)), 'constant')
|
||||
|
||||
signal_src_fft = np.fft.fft(signal_src, fft_len)
|
||||
signal_ref_fft = np.fft.fft(signal_ref, fft_len)
|
||||
|
||||
signal_dst_fft = signal_src_fft * signal_ref_fft
|
||||
|
||||
signal_dst = np.fft.ifft(signal_dst_fft, fft_len)
|
||||
|
||||
return signal_dst
|
||||
|
||||
def lfm(tau, fs, fl, bw):
|
||||
times = np.linspace(0, tau - 1 / fs, int(tau * fs))
|
||||
phase = 2 * np.pi * fl * times + np.pi * bw / tau * times ** 2
|
||||
signal = np.exp(1j * phase)
|
||||
|
||||
return signal
|
||||
|
||||
signal = np.real(lfm(7e-6, 20e6, 222e6, 6e6))
|
||||
|
||||
noise_samples = []
|
||||
|
||||
for _ in range(5000):
|
||||
noise = wgn(signal, 0)
|
||||
noise = hilbert(noise)
|
||||
noise = pulse_compress(noise, noise)
|
||||
|
||||
noise_samples.extend(noise)
|
||||
|
||||
xpower = np.sum(signal ** 2) / len(signal)
|
||||
npower = np.sum(np.abs(noise_samples) ** 2) / len(noise_samples)
|
||||
|
||||
# noise_samples = list(filter(lambda x: np.abs(x) > 0.1, noise_samples))
|
||||
|
||||
noise_samples = np.abs(noise_samples)
|
||||
|
||||
mean = np.mean(noise_samples)
|
||||
var = np.var(noise_samples)
|
||||
print(mean, var)
|
||||
print(xpower, npower)
|
||||
print(np.max(noise_samples))
|
||||
plt.hist(noise_samples, bins=100)
|
||||
|
||||
percentile_95 = np.percentile(noise_samples, 95)
|
||||
percentile_99 = np.percentile(noise_samples, 99)
|
||||
|
||||
# draw a line at 95% of the samples and mark the corresponding value
|
||||
plt.axvline(percentile_95, color='k', linestyle='dashed', linewidth=1)
|
||||
min_ylim, max_ylim = plt.ylim()
|
||||
plt.text(percentile_95 * 1.1, max_ylim * 0.9, f'95% of samples are below {percentile_95:.2f}')
|
||||
|
||||
# 99%
|
||||
plt.axvline(percentile_99, color='k', linestyle='dashed', linewidth=1)
|
||||
min_ylim, max_ylim = plt.ylim()
|
||||
plt.text(percentile_99 * 1.1, max_ylim * 0.8, f'99% of samples are below {percentile_99:.2f}')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
Loading…
Reference in New Issue