From b9143a8db34eba7d9f861d4a5712d6218083f14e Mon Sep 17 00:00:00 2001 From: liuyuanchi Date: Tue, 21 May 2024 15:02:33 +0800 Subject: [PATCH] init --- .gitignore | 8 + Makefile | 31 ++++ README.md | 351 ++++++++++++++++++++++++++++++++++++++++++ include/filter.h | 16 ++ include/hamming.h | 13 ++ include/signal.h | 76 ++++++++++ include/utils.h | 54 +++++++ plot.py | 57 +++++++ src/filter.c | 99 ++++++++++++ src/hamming.c | 28 ++++ src/main.c | 143 ++++++++++++++++++ src/signal.c | 353 +++++++++++++++++++++++++++++++++++++++++++ src/utils.c | 147 ++++++++++++++++++ verification/main.py | 99 ++++++++++++ 14 files changed, 1475 insertions(+) create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 README.md create mode 100755 include/filter.h create mode 100755 include/hamming.h create mode 100755 include/signal.h create mode 100755 include/utils.h create mode 100644 plot.py create mode 100755 src/filter.c create mode 100755 src/hamming.c create mode 100755 src/main.c create mode 100755 src/signal.c create mode 100755 src/utils.c create mode 100644 verification/main.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7ffc170 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +.vscode/ +vsipl/ + +radar + +*.o +*.a +*.so \ No newline at end of file diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..98f2a1b --- /dev/null +++ b/Makefile @@ -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 \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..278c10a --- /dev/null +++ b/README.md @@ -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 大作业。 diff --git a/include/filter.h b/include/filter.h new file mode 100755 index 0000000..655cd5d --- /dev/null +++ b/include/filter.h @@ -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 \ No newline at end of file diff --git a/include/hamming.h b/include/hamming.h new file mode 100755 index 0000000..c1e8631 --- /dev/null +++ b/include/hamming.h @@ -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 \ No newline at end of file diff --git a/include/signal.h b/include/signal.h new file mode 100755 index 0000000..984425f --- /dev/null +++ b/include/signal.h @@ -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 \ No newline at end of file diff --git a/include/utils.h b/include/utils.h new file mode 100755 index 0000000..85e2d91 --- /dev/null +++ b/include/utils.h @@ -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 \ No newline at end of file diff --git a/plot.py b/plot.py new file mode 100644 index 0000000..3ef6648 --- /dev/null +++ b/plot.py @@ -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() diff --git a/src/filter.c b/src/filter.c new file mode 100755 index 0000000..ffef97c --- /dev/null +++ b/src/filter.c @@ -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); +} \ No newline at end of file diff --git a/src/hamming.c b/src/hamming.c new file mode 100755 index 0000000..f8cdc43 --- /dev/null +++ b/src/hamming.c @@ -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); + } +} \ No newline at end of file diff --git a/src/main.c b/src/main.c new file mode 100755 index 0000000..d63095b --- /dev/null +++ b/src/main.c @@ -0,0 +1,143 @@ +#include "filter.h" +#include "hamming.h" +#include "signal.h" +#include "utils.h" +#include "vsip.h" +#include +#include +#include + +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; +} \ No newline at end of file diff --git a/src/signal.c b/src/signal.c new file mode 100755 index 0000000..8092427 --- /dev/null +++ b/src/signal.c @@ -0,0 +1,353 @@ +#include "signal.h" +#include "hamming.h" +#include "utils.h" +#include + +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)); + } + } +} \ No newline at end of file diff --git a/src/utils.c b/src/utils.c new file mode 100755 index 0000000..4e5f2bc --- /dev/null +++ b/src/utils.c @@ -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)); + } +} \ No newline at end of file diff --git a/verification/main.py b/verification/main.py new file mode 100644 index 0000000..6b86ea6 --- /dev/null +++ b/verification/main.py @@ -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() + +