LiuYuanchi
/
phytium-vsipl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

init3

This commit is contained in:
liuyuanchi 2024-05-22 12:25:44 +08:00
parent 795bc5f731
commit fb143565cb
16 changed files with 205 additions and 579 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -6,4 +6,5 @@ radar
*.o
*.a
*.so
data
data
img/

28
Makefile
View File

@ -1,6 +1,7 @@
TARGET = radar
CC = cc
PY = python3
# 注意库的位置
INC_DIR = -I./include -I./src -I../vsipl/include
LIB_DIR = -L./vsipl/lib
LIBS = -lvsip -lfftw3f -lm
@ -9,23 +10,28 @@ CFLAGS = $(INC_DIR) $(LIB_DIR) $(LIBS) -g
SRC = $(wildcard src/*.c)
OBJ = $(patsubst src/%.c, obj/%.o, $(SRC))
# 检查系统架构
ARCH := $(shell uname -m)
ifeq ($(ARCH), x86_64)
BIN_DIR = bin/x86
else ifeq ($(ARCH), armv7l)
BIN_DIR = bin/arm
else ifeq ($(ARCH), aarch64)
BIN_DIR = bin/arm
else
BIN_DIR = bin/unknown
endif
$(TARGET): $(OBJ)
@mkdir -p bin
@$(CC) $(OBJ) $(CFLAGS) -o bin/$(TARGET)
@mkdir -p $(BIN_DIR)
@$(CC) $(OBJ) $(CFLAGS) -o $(BIN_DIR)/$(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 obj
rm -f $(BIN_DIR)/$(TARGET) $(OBJ)
rm -rf obj

349
README.md
View File

@ -1,351 +1,4 @@
# 基于飞腾平台的远距离探测系统
## 项目内容
项目基于 VISPL 函数库实现
项目基于 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 outputRealVector(vsip_vview_f *p_vector, FILE *p_file);
/*
* 内部接口：输出复向量
* 参数：p_vector -- 输入向量
* p_file -- 输出文件
* 功能：将复向量的数据输出到文件
*/
void outputComplexVector(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 大作业。

16
include/filter.h
View File

@ -1,16 +0,0 @@
#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

7
include/hamming.h
View File

@ -3,11 +3,6 @@
#include "vsip.h"
/*
* 内部接口：生成汉明窗
* 参数：p_vector_dst -- 输出信号
* 功能：根据输出信号的长度生成汉明窗
*/
void vcreate_hamming_f(vsip_vview_f *p_vector_dst);
void hamming(vsip_vview_f *p_vector_dst);
#endif

9
include/hilbertFilter.h Executable file
View File

@ -0,0 +1,9 @@
#ifndef RADAR_FILTER_H_
#define RADAR_FILTER_H_
#include "vsip.h"
void hilbert(vsip_vview_f *p_vector_src, vsip_scalar_i n_filter_length,
vsip_cvview_f *p_vector_dst);
#endif

50
include/signal.h
View File

@ -2,74 +2,24 @@
#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);

12
include/tool.h
View File

@ -3,20 +3,8 @@
#include "vsip.h"
/*
* 输出实向量
* 参数：p_vector -- 输入向量
* p_file -- 输出文件
* 功能：将实向量的数据输出到文件
*/
void outputRealVector(vsip_vview_f *p_vector, char *p_name);
/*
* 输出复向量
* 参数：p_vector -- 输入向量
* p_file -- 输出文件
* 功能：将复向量的数据输出到文件
*/
void outputComplexVector(vsip_cvview_f *p_vector,char *p_name);

100
script/plot.py Normal file
View File

@ -0,0 +1,100 @@
import matplotlib.pyplot as plt
import numpy as np
signal = list(map(float, open('./data/signal.txt', 'r').readlines()))
import numpy as np
import matplotlib.pyplot as plt
# 读取并处理信号数据
with open('./data/signal_filtered.txt', 'r') as file:
signal_filtered = [complex(line.strip().replace(' ', '')) for line in file]
# 获取实部和虚部
signal_filtered_real = [z.real for z in signal_filtered]
signal_filtered_imag = [z.imag for z in signal_filtered]
# 绘制复平面
plt.figure(figsize=(8, 8))
# 原始信号复平面展示
plt.scatter(signal_filtered_real, signal_filtered_imag, color='b', label='Filtered Signal')
plt.axhline(0, color='grey', lw=0.5)
plt.axvline(0, color='grey', lw=0.5)
plt.grid(True, which='both', linestyle='--', lw=0.5)
plt.xlabel('Real Part')
plt.ylabel('Imaginary Part')
plt.title('Filtered Signal on Complex Plane')
plt.legend()
# 保存图表
plt.savefig('./img/signal_filtered_complex_plane.png')
plt.clf()
signal_compressed = open('./data/signal_compressed.txt', 'r').readlines()
signal_compressed = map(lambda x: x.replace(' ', '').replace('\n', ''), signal_compressed)
signal_compressed = list(map(complex, signal_compressed))
plt.plot(np.abs(signal_compressed), label='abs',mfc='w',color='b')
plt.legend()
plt.savefig('./img/signal_compressed.png')
plt.clf()
signal_reduced = open('./data/signal_reduced.txt', 'r').readlines()
signal_reduced = map(lambda x: x.replace(' ', '').replace('\n', ''), signal_reduced)
signal_reduced = list(map(complex, signal_reduced))
plt.plot(np.abs(signal_reduced), label='abs',mfc='w',color='b')
plt.legend()
plt.savefig('./img/signal_reduced.png')
plt.clf()
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))
plt.xlim(0, 10)
plt.ylim(-1.2, 1.2)
plt.plot(np.real(hilbert_coefficients), label='real', marker='x', mfc='w',color='r')
plt.plot(np.imag(hilbert_coefficients), label='imag', marker='x', mfc='w',color='b')
plt.legend()
plt.savefig('./img/hilbert_coefficients.png')
plt.clf()
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.ylim(-1.2, 1.2)
plt.plot(np.real(reference_signal), label='real' ,mfc='w',color='b')
plt.legend()
plt.savefig('./img/reference_signal.png')
plt.clf()

43
verification/main.py → script/test.py
View File

@ -1,21 +1,22 @@
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
def wgn(signal, snr = 0):
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):
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')
@ -26,7 +27,6 @@ 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)
@ -36,7 +36,7 @@ def pulse_compress(signal_src, 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')
@ -60,7 +60,7 @@ signal = np.real(lfm(7e-6, 20e6, 222e6, 6e6))
noise_samples = []
for _ in range(5000):
for _ in range(100):
noise = wgn(signal, 0)
noise = hilbert(noise)
noise = pulse_compress(noise, noise)
@ -70,8 +70,6 @@ for _ in range(5000):
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)
@ -79,21 +77,24 @@ 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)
# 概率密度函数（PDF）图
plt.figure()
sns.kdeplot(noise_samples, bw_adjust=0.5)
plt.title('Probability Density Function (PDF)')
plt.savefig('./img/noise_pdf.png')
# 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}')
# 累积分布函数（CDF）图
plt.figure()
sns.ecdfplot(noise_samples)
plt.title('Cumulative Distribution Function (CDF)')
plt.savefig('./img/noise_cdf.png')
# 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}')
# 箱线图（Box Plot）
plt.figure()
sns.boxplot(x=noise_samples)
plt.title('Box Plot')
plt.savefig('./img/noise_boxplot.png')
# 显示所有图像
plt.show()

4
src/hamming.c
View File

@ -1,10 +1,10 @@
#include "hamming.h"
void vcreate_hamming_f(vsip_vview_f *p_vector_dst)
void hamming(vsip_vview_f *p_vector_dst)
{
if (p_vector_dst == NULL)
{
fprintf(stderr, "vcreate_hamming_f: p_vector_dst is NULL\n");
fprintf(stderr, "hamming: p_vector_dst is NULL\n");
return;
}

2
src/filter.c → src/hilbertFilter.c
View File

@ -1,4 +1,4 @@
#include "filter.h"
#include "hilbertFilter.h"
#include "tool.h"
void hilbert(vsip_vview_f *p_vector_src, vsip_scalar_i n_filter_length, vsip_cvview_f *p_vector_dst)

98
src/main.c
View File

@ -1,117 +1,112 @@
#include "filter.h"
#include "hilbertFilter.h"
#include "hamming.h"
#include "signal.h"
#include "tool.h"
#include "vsip.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
const vsip_scalar_f pulse_width = 7e-6; // 脉冲宽度
const vsip_scalar_f lower_limit_frequency = 222e6; // 下限频率
const vsip_scalar_f band_width = 6e6; // 带宽
const vsip_scalar_f sampling_rate = 20e6; // 采样率
const vsip_scalar_f distance = 600.0f; // 目标距离
const vsip_length signal_length = (vsip_length)(0.5f + pulse_width * sampling_rate); // LFM 信号长度
const vsip_scalar_f delay_time = 2.0f * distance / 3e8; // 延迟时间
const vsip_scalar_f total_time = pulse_width + delay_time; // 总时间
const vsip_length delay_signal_length = (vsip_length)(0.5f + delay_time * sampling_rate); // 延迟信号长度（样本数量）
const vsip_length total_signal_length = (vsip_length)(0.5f + total_time * sampling_rate); // 总信号长度（样本数量）
const vsip_length hilbert_length = 11; // 希尔伯特系数
int main(int argc, char *argv[])
{
// 初始化VSIPL库
vsip_init((void *)0);
const vsip_scalar_f pulse_width = 7e-6; // 脉冲宽度
const vsip_scalar_f lower_limit_frequency = 222e6; // 下限频率
const vsip_scalar_f band_width = 6e6; // 带宽
const vsip_scalar_f sampling_rate = 20e6; // 采样率
const vsip_scalar_f distance = 600.0f; // 目标距离
const vsip_length signal_length = (vsip_length)(0.5f + pulse_width * sampling_rate); // LFM 信号长度
const vsip_scalar_f delay_time = 2.0f * distance / 3e8; // 延迟时间
const vsip_scalar_f total_time = pulse_width + delay_time; // 总时间
const vsip_length delay_signal_length = (vsip_length)(0.5f + delay_time * sampling_rate); // 延迟信号长度（样本数量）
const vsip_length total_signal_length = (vsip_length)(0.5f + total_time * sampling_rate); // 总信号长度（样本数量）
const vsip_length hilbert_length = 11; // 希尔伯特系数
// 生成雷达接收的实信号
// 创建一个用于存储雷达信号的实数向量
vsip_vview_f *p_vector_radar_signal = vsip_vcreate_f(total_signal_length, VSIP_MEM_NONE);
// 生成雷达信号
generate_radar_signal(pulse_width, sampling_rate, lower_limit_frequency, band_width, 600.0f,
p_vector_radar_signal);
// DEBUG
outputRealVector(p_vector_radar_signal, "./data/radar_signal.txt");
// 将实数雷达信号输出到文件
outputRealVector(p_vector_radar_signal, "./data/signal.txt");
// 希尔伯特滤波
// 创建一个用于存储经过希尔伯特变换的雷达信号的复数向量
vsip_cvview_f *p_vector_radar_signal_filtered =
vsip_cvcreate_f(total_signal_length, VSIP_MEM_NONE);
// 对雷达信号进行希尔伯特变换
hilbert(p_vector_radar_signal, 11, p_vector_radar_signal_filtered);
// DEBUG
outputComplexVector(p_vector_radar_signal_filtered, "./data/radar_signal_filtered.txt");
// 将经过希尔伯特变换的雷达信号输出到文件
outputComplexVector(p_vector_radar_signal_filtered, "./data/signal_filtered.txt");
// 匹配滤波参考信号
// 创建一个用于存储参考信号的复数向量
vsip_cvview_f *p_vector_signal_ref = vsip_cvcreate_f(signal_length, VSIP_MEM_NONE);
// 生成参考信号
generate_lfm_signal(pulse_width, sampling_rate, lower_limit_frequency, band_width, p_vector_signal_ref);
// 傅里叶变换长度
// 计算进行脉冲压缩的快速傅里叶变换长度
vsip_length fft_len = 2 * total_signal_length - 1;
fft_len = (vsip_length)ceil(0.5 + log2(fft_len));
fft_len = (vsip_length)floor(0.5 + pow(2, fft_len));
// DEBUG
printf("radar: fft_len = %ld\n", fft_len);
// 脉冲压缩
// 创建一个用于存储经过脉冲压缩的雷达信号的复数向量
vsip_cvview_f *p_vector_radar_signal_compressed = vsip_cvcreate_f(fft_len, VSIP_MEM_NONE);
// 对雷达信号进行脉冲压缩
pulse_compress(p_vector_radar_signal_filtered, p_vector_signal_ref,
p_vector_radar_signal_compressed);
// DEBUG
outputComplexVector(p_vector_radar_signal_compressed, "./data/radar_signal_compressed.txt");
// 将经过脉冲压缩的雷达信号输出到文件
outputComplexVector(p_vector_radar_signal_compressed, "./data/signal_compressed.txt");
// 根据 30 阈值筛选
// 创建一个用于存储经过阈值检测的雷达信号的复数向量
vsip_cvview_f *p_vector_radar_signal_reduced = vsip_cvcreate_f(fft_len, VSIP_MEM_NONE);
detect_signal(p_vector_radar_signal_compressed, 30.0f, p_vector_radar_signal_reduced);
// 对雷达信号进行阈值检测
detect_signal(p_vector_radar_signal_compressed, 20.0f, p_vector_radar_signal_reduced);
// DEBUG
outputComplexVector(p_vector_radar_signal_reduced, "./data/radar_signal_reduced.txt");
// 将经过阈值检测的雷达信号输出到文件
outputComplexVector(p_vector_radar_signal_reduced, "./data/signal_reduced.txt");
// 检测目标并且输出距离
// 采样间隔
// 计算时间差
vsip_scalar_f delta_time = 1.0f / sampling_rate;
// 上一个峰值的时间
vsip_scalar_f prev_time = 0.0f;
// 在向量中的下标
vsip_length index = 0;
// 遍历经过阈值检测的雷达信号,寻找目标
while (index < fft_len - 1)
{
// 当前时间
vsip_scalar_f curr_time = (vsip_scalar_f)index * delta_time;
// 下一个时间
vsip_scalar_f next_time = (vsip_scalar_f)(index + 1) * delta_time;
// 当前幅值
vsip_scalar_f curr_mag =
vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index));
// 下一个幅值
vsip_scalar_f next_mag =
vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index + 1));
vsip_scalar_f curr_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index));
vsip_scalar_f next_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index + 1));
// 如果当前幅值小于等于下一个幅值,则继续遍历
if (curr_mag <= next_mag)
{
// 若仍然处于递增阶段，则继续向后搜索
index++;
}
else
{
// 当前幅值大于前一个幅值，且大于后一个幅值，说明是一个峰值
// 如果当前幅值大于下一个幅值,则说明找到了一个目标
if (prev_time == 0.0f)
{
// 第一个目标，记录时间
// 如果是第一个目标,记录时间
prev_time = curr_time;
printf("detect: time = %fs\n", curr_time);
printf("First record: %fs\n", curr_time);
}
else
{
// 与上一个目标的距离
// 如果不是第一个目标,计算距离并记录时间
vsip_scalar_f distance = 3e8 * (curr_time - prev_time) / 2.0f;
// 记录时间
prev_time = curr_time;
printf("detect: time = %fs, distance = %fm\n", curr_time, distance);
printf("Second record: %fs, \nDistance: %fm\n", curr_time, distance);
}
// 继续寻找目标
while (curr_mag > next_mag)
{
// 在递减阶段，继续向后搜索
index++;
curr_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index));
next_mag = vsip_cmag_f(vsip_cvget_f(p_vector_radar_signal_reduced, index + 1));
@ -126,6 +121,7 @@ int main(int argc, char *argv[])
vsip_valldestroy_f(p_vector_radar_signal);
vsip_cvalldestroy_f(p_vector_radar_signal_filtered);
// 结束VSIPL库
vsip_finalize((void *)0);
return 0;

2
src/signal.c
View File

@ -329,7 +329,7 @@ void pulse_compress(vsip_cvview_f *p_vector_signal_src, vsip_cvview_f *p_vector_
cvflip_f(p_vector_signal_ref, p_vector_signal_ref_flipped);
// 加汉明窗
vsip_vview_f *p_vector_window = vsip_vcreate_f(signal_ref_length, VSIP_MEM_NONE);
vcreate_hamming_f(p_vector_window);
hamming(p_vector_window);
vsip_rcvmul_f(p_vector_window, p_vector_signal_ref_flipped, p_vector_signal_ref);

4
src/tool.c
View File

@ -7,7 +7,7 @@ void outputRealVector(vsip_vview_f *p_vector, char *p_name)
fprintf(stderr, "outputRealVector: p_vector is NULL\n");
return;
}
FILE *p_file = fopen(p_name, "w");
FILE *p_file = fopen(p_name, "w+");
if (p_file == NULL)
{
fprintf(stderr, "vdebug_f: unable to open file `%s`\n", p_name);
@ -29,7 +29,7 @@ void outputComplexVector(vsip_cvview_f *p_vector,char *p_name)
fprintf(stderr, "outputRealVector: p_vector is NULL\n");
return;
}
FILE *p_file = fopen(p_name, "w");
FILE *p_file = fopen(p_name, "w+");
if (p_file == NULL)
{
fprintf(stderr, "cvdebug_f: unable to open file `%s`\n", p_name);

57
verification/plot.py
View File

@ -1,57 +0,0 @@
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()