STM32生成高斯白噪声
白噪声
功率谱密度均匀分布在整个频率范围内的随机过程成为白噪声。
白噪声是一种特殊的随机信号，具有以下关键特性：
频谱平坦性：在所有频率上具有相等的功率密度
不相关性：任意两个不同时刻的噪声值之间没有相关性
均匀分布：瞬时值通常服从均匀分布或高斯分布
类比：白噪声类似于电视无信号时的"雪花"画面，或收音机调频时的"嘶嘶"声。
白噪声的频率特性
频率范围：理论上包含所有频率成分（从0到无穷大）
功率分布：各频率分量功率相等（理想状态）
实际限制：受采样率和硬件限制，实际白噪声有上限频率
白噪声 vs 粉噪声 vs 布朗噪声
类型 功率谱特性 听觉感受 典型应用
白噪声 所有频率功率相等 "嘶嘶"声 音频测试、耳鸣治疗
粉噪声 功率与频率成反比 "雨声"感 声学测试、睡眠辅助
布朗噪声 功率与频率平方成反比 "雷声"感 低频噪声模拟
高斯白噪声
服从高斯分布的白噪声即为高斯白噪声。
1.生成高斯白噪声
Box-Muller 变换是一种从均匀分布随机数高效生成高斯分布随机数的算法，相比 “中心极限定理近似法”，它生成的随机数更接近理想高斯分布，统计特性更稳定。
因此使用 Box-Muller 变换算法实现白噪声，用于生成均值为 0、标准差为 1 的高斯分布（正态分布）随机数。
// Box-Muller高斯噪声生成所需静态变量
static uint8_t has_spare = 0; // 标记是否有"备用"的高斯随机数
static float spare; // 存储"备用"的高斯随机数
// 生成高斯分布随机数 - Box-Muller变换法
// 该算法每次生成两个独立的高斯随机数，第一个直接返回，第二个保存至变量以供下次直接调用
float Random_Gaussian(void) {
float u, v, s; // 临时变量：u、v是均匀分布随机数，s是中间计算值
// 检查是否有"备用"的高斯随机数
if (has_spare) {
has_spare = 0; // 清除备用标志
return spare; // 返回上次生成的备用随机数
}
// 生成符合条件的均匀分布随机数对
do {
// 生成[-1, 1)范围的均匀分布随机数
// Random_Generate()返回0~0x7FFFFFFF的均匀数，除以0x3FFFFFFF后范围是[0, 2)，减1后得到[-1, 1)
u = (float)Random_Generate() / 0x3FFFFFFF - 1.0f;
v = (float)Random_Generate() / 0x3FFFFFFF - 1.0f;
// 计算u和v的平方和（用于判断是否在单位圆内）
s = u * u + v * v;
} while (s >= 1.0f || s == 0.0f); // 确保(u,v)在单位圆内（s < 1）且不为原点（s ≠ 0）
// Box-Muller变换核心计算：将均匀分布转换为高斯分布
s = sqrtf(-2.0f * logf(s) / s); // 计算缩放因子
// 生成两个独立的高斯随机数
spare = v * s; // 保存一个作为"备用"，供下一次调用
has_spare = 1; // 标记有备用随机数
return u * s; // 返回当前的高斯随机数
}
// 生成随机数-为高斯分布噪声提供原始的均匀分布随机数
uint32_t Random_Generate(void) {
static uint32_t seed = 1; // 静态种子变量，保存当前状态
seed = (1103515245 * seed + 12345) % 0x7FFFFFFF; // 生成一个范围在0 ~ 0x7FFFFFFF（即 0 到 2147483647）之间的伪随机数，且每次调用会更新内部状态，确保后续生成的数值具有随机性。
//线性同余生成器（LCG） 算法的实现，公式为：新seed = (a × 旧seed + c) % m
//其中，a = 1103515245（乘数）：经过数学验证的最优值，确保生成的随机数分布均匀。c = 12345（增量）：非零值，避免生成序列退化。m = 0x7FFFFFFF（模数）：即 2³¹-1，一个大质数，决定随机数的范围（0 到 m-1），同时保证序列周期足够长（理论上可达 m-1）。
return seed; // 返回生成的均匀随机数
}
在本地通过python检验生成的高斯白噪声是否符合规范
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <time.h>
#include <string.h>
#define PI 3.14159265358979323846
...
int main()
{
const int N = 4096; // 样本点数
const char* filename = "noise_signal.bin";
// 分配内存
double *gauss_rand = (double*)malloc(N * sizeof(double));
// 1. 生成高斯白噪声
for (int i = 0; i < N; i++) {
gauss_rand[i] = Random_Gaussian();
}
// 4. 写入文件
FILE *fp = fopen(filename, "wb");
for (int i = 0; i < N; i++) {
float gauss = (float)gauss_rand[i];
fwrite(&gauss, sizeof(float), 1, fp);
}
fclose(fp);
return 0;
}
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.stats import norm
import os
# 设置中文字体，确保中文正常显示
plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False # 正确显示负号
def read_gaussian_noise(filename):
"""读取bin文件中的高斯噪声数据（float32格式）"""
if not os.path.exists(filename):
raise FileNotFoundError(f"文件 {filename} 不存在")
# 读取bin文件，数据格式为float32
noise_data = np.fromfile(filename, dtype=np.float32)
print(f"成功读取数据：共 {len(noise_data)} 个样本")
return noise_data
def plot_noise_analysis(noise_data, sample_rate=1000):
"""绘制时域图、频域图和分布直方图，调整标题对齐方式"""
n = len(noise_data)
duration = n / sample_rate # 总时长（秒）
t = np.linspace(0, duration, n, endpoint=False) # 时间轴
# 创建3行1列的图表
fig, axs = plt.subplots(3, 1, figsize=(10, 15))
# 调整主标题位置，避免与子图标题重叠（y值越大越靠上）
#fig.suptitle('高斯白噪声特性分析', fontsize=16, y=0.98)
# 1. 时域图（标题左对齐）
axs[0].plot(t, noise_data, linewidth=0.8)
# loc='left' 左对齐；loc='right' 右对齐；pad控制标题与子图顶部的距离
axs[0].set_title('时域波形', loc='left', pad=0, fontsize=12)
axs[0].set_xlabel('时间 (秒)')
axs[0].set_ylabel('幅度')
axs[0].grid(True, linestyle='--', alpha=0.7)
# 2. 频域图（FFT，标题左对齐）
yf = fft(noise_data)
xf = fftfreq(n, 1 / sample_rate)[:n//2] # 正频率轴
amplitude = 2.0 / n * np.abs(yf[:n//2]) # 幅度归一化
axs[1].plot(xf, amplitude, linewidth=0.8)
axs[1].set_title('频域特性（FFT）', loc='left', pad=0, fontsize=12)
axs[1].set_xlabel('频率 (Hz)')
axs[1].set_ylabel('幅度')
axs[1].set_xlim(0, sample_rate/2) # 仅显示0到奈奎斯特频率
axs[1].grid(True, linestyle='--', alpha=0.7)
# 3. 分布直方图（高斯验证，标题左对齐）
axs[2].hist(noise_data, bins=50, density=True, alpha=0.7, color='skyblue', label='数据分布')
# 拟合理论高斯分布曲线
mu, sigma = np.mean(noise_data), np.std(noise_data)
# 处理标准差为0的极端情况（避免绘图错误）
if sigma < 1e-9:
sigma = 1e-9
x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
axs[2].plot(x, norm.pdf(x, mu, sigma), 'r-', linewidth=2,
label=f'高斯拟合 (μ={mu:.2f}, σ={sigma:.2f})')
axs[2].set_title('幅度分布直方图', loc='left', pad=0, fontsize=12)
axs[2].set_xlabel('幅度值')
axs[2].set_ylabel('概率密度')
axs[2].grid(True, linestyle='--', alpha=0.7)
axs[2].legend()
plt.tight_layout()
# 调整子图与顶部的间距（配合suptitle的y值）
plt.subplots_adjust(top=0.98)
plt.show()
if __name__ == "__main__":
# 配置参数（根据实际情况修改）
filename = "noise_signal.bin" # bin文件路径
sample_rate = 1000 # 采样率（Hz），可根据实际生成时的配置修改
try:
# 读取数据
noise_data = read_gaussian_noise(filename)
# 绘制分析图表
plot_noise_analysis(noise_data, sample_rate)
except Exception as e:
print(f"出错：{e}")
可以看到符合要求
2.低通/带通滤波器
截取指定频段的高斯白噪声，就需要用到滤波器。常见的有IIR型和FIR型滤波器
IIR更适合在嵌入式平台，尤其像stm32这样计算资源有限的设备上。（IIR无法用在要求相位线性的情况下，但白噪声无需相位线性）
| IIR型 | 特点 |
|---|---|
| 巴特沃斯滤波器 | 通带、阻带平坦，几乎无波纹，过渡带最宽 |
| 切比雪夫 I 型滤波器 | 通带有波纹，阻带平坦，过渡带较窄 |
| 切比雪夫 II 型滤波器 | 通带平坦，阻带有波纹，过渡带较窄 |
| 椭圆滤波器 | 通带、阻带有波纹，过渡带最窄 |
巴特沃斯在通带平坦、几乎无波纹，很符合高斯白噪声 目标频段内功率谱均匀 的要求。此外，使用高阶巴特沃斯可减小过渡带。因此此处选用巴特沃斯滤波器。
3.STM32实现
#include "stm32f10x.h"
#include "nGVS.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
// 系统参数配置
#define SAMPLE_RATE 1000 // 采样率1000Hz
#define DEFAULT_LOW_CUT 0 // 默认低频截止0Hz
#define DEFAULT_HIGH_CUT 30 // 默认高频截止30Hz
#define DAC_OFFSET 2047 // 对应0mA的中间值
#define MAX_DAC_VALUE 4095 // 对应+4mA
#define MIN_DAC_VALUE 0 // 对应-4mA
#define TARGET_AMPLITUDE 0.4f // 目标峰峰值振幅1mA
#define MAX_CURRENT 4.0f // 最大电流范围±4mA
#define CALIBRATION_SAMPLES 1000 // 校准样本数量
// 滤波器结构体定义
typedef struct {
float low_cut; // 低频截止频率
float high_cut; // 高频截止频率
float sample_rate; // 采样率
float a[3]; // 滤波器系数a
float b[3]; // 滤波器系数b
float x[3]; // 输入历史
float y[3]; // 输出历史
} FilterTypeDef;
// 全局变量
uint16_t low_cut = DEFAULT_LOW_CUT;
uint16_t high_cut = DEFAULT_HIGH_CUT;
FilterTypeDef iirFilter; // 滤波器实例
float noise_gain = 1.0f; // 噪声幅度增益（动态校准）
float target_peak = TARGET_AMPLITUDE / 2.0f; // 目标峰值（峰峰值/2）
bool is_calibrated = false; // 校准状态标志
// Box-Muller高斯噪声生成所需静态变量
static uint8_t has_spare = 0;
static float spare;
// 函数声明
void RCC_Configuration(void);
void GPIO_Configuration(void);
void TIM_Configuration(void);
void NVIC_Configuration(void);
uint32_t Random_Generate(void);
float Random_Gaussian(void);
void IIR_Filter_Init(FilterTypeDef *filter, float sample_rate, float low_cut, float high_cut);
float IIR_Filter_Process(FilterTypeDef *filter, float input);
void Set_Noise_Amplitude(float amplitude_mA);
void Calibrate_Filter_Gain(void);
void Send_Text(float);
// 主函数
int nGVS(void) {
// 初始化系统
RCC_Configuration();
GPIO_Configuration();
TIM_Configuration();
NVIC_Configuration();
// 初始化随机数种子
srand((uint32_t)TIM4->CNT);
// 初始化滤波器
IIR_Filter_Init(&iirFilter, SAMPLE_RATE, low_cut, high_cut);
// 执行滤波器增益校准
Calibrate_Filter_Gain();
// 设置目标振幅
Set_Noise_Amplitude(TARGET_AMPLITUDE);
// 启动定时器
TIM_Cmd(TIM4, ENABLE);
// 主循环
while (1) {
__WFI(); // 等待中断
}
}
// 定时器4中断服务程序
void TIM4_IRQHandler(void) {
static uint32_t sample_count = 0;
static float max_sample = 0.0f;
if (TIM_GetITStatus(TIM4, TIM_IT_Update) != RESET) {
TIM_ClearITPendingBit(TIM4, TIM_IT_Update);
// 1. 生成高斯分布随机数
float noise = Random_Gaussian();
// 2. 应用滤波器
float filtered_noise = IIR_Filter_Process(&iirFilter, noise);
// 3. 应用校准增益
float scaled_noise = filtered_noise * noise_gain;
// 4. 限幅保护
if (scaled_noise > target_peak) scaled_noise = target_peak;
if (scaled_noise < -target_peak) scaled_noise = -target_peak;
// 5. 转换为0-4095范围值
uint16_t dac_equivalent = (uint16_t)(
scaled_noise / MAX_CURRENT * DAC_OFFSET + DAC_OFFSET
);
// 确保在有效范围内
if (dac_equivalent > MAX_DAC_VALUE) dac_equivalent = MAX_DAC_VALUE;
if (dac_equivalent < MIN_DAC_VALUE) dac_equivalent = MIN_DAC_VALUE;
// 6. 输出数据
Send_Text(scaled_noise);
}
}
// 滤波器增益校准：计算不同带宽下的幅度衰减并补偿
void Calibrate_Filter_Gain(void) {
float max_amplitude = 0.0f;
FilterTypeDef temp_filter;
uint32_t i;
// 初始化临时滤波器（与主滤波器配置相同）
IIR_Filter_Init(&temp_filter, SAMPLE_RATE, low_cut, high_cut);
// 生成校准样本并测量滤波后的最大幅度
for (i = 0; i < CALIBRATION_SAMPLES; i++) {
float noise = Random_Gaussian();
float filtered = IIR_Filter_Process(&temp_filter, noise);
// 跟踪最大值（绝对值）
if (fabs(filtered) > max_amplitude) {
max_amplitude = fabs(filtered);
}
// 简单延时，确保校准过程稳定
for (volatile uint32_t d = 0; d < 100; d++);
}
// 计算校准增益：目标峰值 / 测量到的最大幅度
// 增加安全系数1.2防止过度放大
if (max_amplitude > 0.001f) {
noise_gain = (target_peak / max_amplitude) / 1.2f;
is_calibrated = true;
} else {
// 校准失败时使用默认增益
noise_gain = target_peak / 3.0f; // 基于高斯分布3σ特性
is_calibrated = false;
}
}
// 设置噪声幅度(mA峰峰值)
void Set_Noise_Amplitude(float amplitude_mA) {
if (amplitude_mA <= 0 || amplitude_mA > 2 * MAX_CURRENT) {
amplitude_mA = TARGET_AMPLITUDE; // 限制在有效范围
}
target_peak = amplitude_mA / 2.0f; // 峰峰值转峰值
// 如果已校准，保持增益不变；否则重新校准
if (!is_calibrated) {
Calibrate_Filter_Gain();
}
}
// 更新滤波器参数并重新校准
void Update_Filter(uint16_t new_low, uint16_t new_high) {
low_cut = new_low;
high_cut = new_high;
IIR_Filter_Init(&iirFilter, SAMPLE_RATE, low_cut, high_cut);
Calibrate_Filter_Gain(); // 带宽改变后必须重新校准
}
// 系统时钟配置
void RCC_Configuration(void) {
RCC_HSEConfig(RCC_HSE_ON);
while (RCC_GetFlagStatus(RCC_FLAG_HSERDY) == RESET);
RCC_PLLConfig(RCC_PLLSource_HSE_Div1, RCC_PLLMul_9);
RCC_PLLCmd(ENABLE);
while (RCC_GetFlagStatus(RCC_FLAG_PLLRDY) == RESET);
RCC_SYSCLKConfig(RCC_SYSCLKSource_PLLCLK);
RCC_HCLKConfig(RCC_SYSCLK_Div1);
RCC_PCLK2Config(RCC_HCLK_Div1);
RCC_PCLK1Config(RCC_HCLK_Div2);
RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA | RCC_APB2Periph_AFIO, ENABLE);
RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM4, ENABLE);
}
// GPIO配置
void GPIO_Configuration(void) {
// 配置串口相关GPIO（根据实际硬件修改）
}
// 定时器配置
void TIM_Configuration(void) {
TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;
TIM_TimeBaseStructure.TIM_Period = (1000000 / SAMPLE_RATE) - 1;
TIM_TimeBaseStructure.TIM_Prescaler = 72 - 1;
TIM_TimeBaseStructure.TIM_ClockDivision = 0;
TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
TIM_TimeBaseInit(TIM4, &TIM_TimeBaseStructure);
TIM_ITConfig(TIM4, TIM_IT_Update, ENABLE);
}
// 中断优先级配置
void NVIC_Configuration(void) {
NVIC_InitTypeDef NVIC_InitStructure;
NVIC_PriorityGroupConfig(NVIC_PriorityGroup_2);
NVIC_InitStructure.NVIC_IRQChannel = TIM4_IRQn;
NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 1;
NVIC_InitStructure.NVIC_IRQChannelSubPriority = 0;
NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;
NVIC_Init(&NVIC_InitStructure);
}
// 生成均匀分布随机数
uint32_t Random_Generate(void) {
static uint32_t seed = 1;
seed = (1103515245 * seed + 12345) % 0x7FFFFFFF;
return seed;
}
// 生成高斯分布随机数
float Random_Gaussian(void) {
float u, v, s;
if (has_spare) {
has_spare = 0;
return spare;
}
do {
u = (float)Random_Generate() / 0x3FFFFFFF - 1.0f;
v = (float)Random_Generate() / 0x3FFFFFFF - 1.0f;
s = u * u + v * v;
} while (s >= 1.0f || s == 0.0f);
s = sqrtf(-2.0f * logf(s) / s);
spare = v * s;
has_spare = 1;
return u * s;
}
// 初始化IIR滤波器
void IIR_Filter_Init(FilterTypeDef *filter, float sample_rate, float low_cut, float high_cut) {
filter->sample_rate = sample_rate;
filter->low_cut = low_cut;
filter->high_cut = high_cut;
filter->x[0] = filter->x[1] = filter->x[2] = 0.0f;
filter->y[0] = filter->y[1] = filter->y[2] = 0.0f;
if (low_cut <= 0.001f) {
// 低通滤波器配置
float wc = 2 * PI * high_cut / sample_rate;
float alpha = sinf(wc) / (2 * 0.707f);
filter->b[0] = (1 - cosf(wc)) / 2;
filter->b[1] = 1 - cosf(wc);
filter->b[2] = (1 - cosf(wc)) / 2;
filter->a[0] = 1 + alpha;
filter->a[1] = -2 * cosf(wc);
filter->a[2] = 1 - alpha;
} else {
// 带通滤波器配置
float w0 = 2 * PI * (low_cut + high_cut) / (2 * sample_rate);
float bw = 2 * PI * (high_cut - low_cut) / sample_rate;
float alpha = sinf(w0) * sinhf(logf(2)/2 * bw / sinf(w0));
filter->b[0] = alpha;
filter->b[1] = 0;
filter->b[2] = -alpha;
filter->a[0] = 1 + alpha;
filter->a[1] = -2 * cosf(w0);
filter->a[2] = 1 - alpha;
}
// 归一化系数
filter->b[0] /= filter->a[0];
filter->b[1] /= filter->a[0];
filter->b[2] /= filter->a[0];
filter->a[1] /= filter->a[0];
filter->a[2] /= filter->a[0];
}
// IIR滤波器处理
float IIR_Filter_Process(FilterTypeDef *filter, float input) {
filter->x[2] = filter->x[1];
filter->x[1] = filter->x[0];
filter->x[0] = input;
float output = filter->b[0] * filter->x[0] +
filter->b[1] * filter->x[1] +
filter->b[2] * filter->x[2] -
filter->a[1] * filter->y[1] -
filter->a[2] * filter->y[2];
filter->y[2] = filter->y[1];
filter->y[1] = output;
return output;
}
// 发送数据到串口
void Send_Text(float dac){
char sendtext[20];
sprintf(sendtext, "{\"INA219\":{\"ma\":%.2f, \"v\":%.2f}}", dac, INA219_Voltage_V());
Serial_SendString(sendtext);
}
浙公网安备 33010602011771号