[实战]带限高斯白噪声的C语言实现（完整代码）

带限高斯白噪声的C语言实现（完整代码）

不依赖任何第三方库，用C语言生成带限高斯白噪声，并使用python进行数据频谱分析，确认C语言实现的正确性。
话不多说，实战开始。

---

C语言实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

#define PI 3.14159265358979323846

// 生成高斯随机数
double gauss_rand(double sigma) {
    double u1 = rand() / (RAND_MAX + 1.0);
    double u2 = rand() / (RAND_MAX + 1.0);
    return sigma * sqrt(-2.0 * log(u1)) * cos(2 * PI * u2);
}

// 计算sinc函数
double sinc(double x) {
    if (fabs(x) < 1e-6) return 1.0;
    return sin(PI * x) / (PI * x);
}

int main() {
    // 参数设置
    const double Fs = 1000.0;       // 采样率 (Hz)
    const double BW = 100.0;        // 带宽 (Hz)
    const int N = 4096;             // 样本点数
    const double sigma = 1.0;       // 高斯分布标准差
    const char* filename = "noise_signal.bin";
    
    // 滤波器参数
    const int filter_length = 101;  // 滤波器长度 (奇数)
    const double cutoff = BW / Fs;  // 归一化截止频率

    srand(time(NULL));  // 初始化随机种子
    
    // 分配内存
    double *real = (double*)malloc(N * sizeof(double));
    double *imag = (double*)malloc(N * sizeof(double));
    double *filter = (double*)malloc(filter_length * sizeof(double));
    
    // 1. 生成高斯白噪声
    for (int i = 0; i < N; i++) {
        real[i] = gauss_rand(sigma);
        imag[i] = gauss_rand(sigma);
    }
    
    // 2. 设计低通滤波器 (sinc函数)
    int center = filter_length / 2;
    for (int i = 0; i < filter_length; i++) {
        int idx = i - center;
        filter[i] = 2.0 * cutoff * sinc(2.0 * cutoff * idx);
        
        // 应用汉明窗
        filter[i] *= 0.54 - 0.46 * cos(2 * PI * i / (filter_length - 1));
    }
    
    // 3. 滤波处理 (卷积)
    double *filtered_real = (double*)malloc(N * sizeof(double));
    double *filtered_imag = (double*)malloc(N * sizeof(double));
    memset(filtered_real, 0, N * sizeof(double));
    memset(filtered_imag, 0, N * sizeof(double));
    
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < filter_length; j++) {
            int pos = i - center + j;
            if (pos >= 0 && pos < N) {
                filtered_real[i] += real[pos] * filter[j];
                filtered_imag[i] += imag[pos] * filter[j];
            }
        }
    }
    
    // 4. 写入文件
    FILE *fp = fopen(filename, "wb");
    for (int i = 0; i < N; i++) {
        float r = (float)filtered_real[i];
        float im = (float)filtered_imag[i];
        fwrite(&r, sizeof(float), 1, fp);
        fwrite(&im, sizeof(float), 1, fp);
    }
    fclose(fp);
    
    // 清理资源
    free(real);
    free(imag);
    free(filtered_real);
    free(filtered_imag);
    free(filter);

    printf("Generated %d samples with BW=%.1fHz\n", N, BW);
    printf("Used convolution method with %d-tap filter\n", filter_length);
    return 0;
}

---

Python分析代码

import numpy as np
import matplotlib.pyplot as plt

# 参数设置 (必须与C程序匹配)
Fs = 1000.0  # 采样率 (Hz)
N = 4096     # 样本点数
filename = "noise_signal.bin"

# 读取二进制文件
data = np.fromfile(filename, dtype=np.float32)
complex_data = data[::2] + 1j * data[1::2]  # 重组为复数

# 计算功率
power = np.mean(np.abs(complex_data) ** 2)
print(f"信号功率: {power:.4f} W")
print(f"实部标准差: {np.std(complex_data.real):.4f}")
print(f"虚部标准差: {np.std(complex_data.imag):.4f}")

# 计算频谱
freqs = np.fft.fftshift(np.fft.fftfreq(N, 1/Fs))
spectrum = np.fft.fftshift(np.fft.fft(complex_data))
psd = np.abs(spectrum) ** 2 / (Fs * N)  # 功率谱密度

# 绘制时域信号
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(np.real(complex_data)[:200], label='Real')
plt.plot(np.imag(complex_data)[:200], label='Imag')
plt.title('Time Domain Signal (First 200 samples)')
plt.xlabel('Sample Index')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

# 绘制频谱
plt.subplot(3, 1, 2)
plt.plot(freqs, 10 * np.log10(psd + 1e-10))  # 避免log(0)
plt.title('Power Spectral Density')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (dB/Hz)')
plt.xlim(-Fs/2, Fs/2)
plt.grid(True)

# 绘制直方图 (验证高斯分布)
plt.subplot(3, 1, 3)
plt.hist(np.real(complex_data), bins=50, alpha=0.7, label='Real')
plt.hist(np.imag(complex_data), bins=50, alpha=0.7, label='Imag')
plt.title('Amplitude Distribution')
plt.xlabel('Amplitude')
plt.ylabel('Count')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.savefig('noise_analysis.png')
plt.show()

结果分析

通过python分析，结果如下

信号功率: 0.1231 W
实部标准差: 0.2451
虚部标准差: 0.2503

时域、频域与分布情况

证实信号生成是正确可用的。

编译和运行说明

1. 编译C程序 (无需外部库):

   gcc bandlimited_noise.c -o bandlimited_noise -lm
   

2. 运行C程序生成数据:

   ./bandlimited_noise
   

3. 运行Python分析脚本:

   python analyze_noise.py
   

新方法的关键特性

1. 无需FFTW库:

   • 使用时域卷积代替频域滤波
   • 完全避免外部库依赖

2. 高斯噪声生成:

   • 使用Box-Muller方法生成高斯随机数
   • 实部和虚部独立生成

3. FIR滤波器设计:

   • 基于sinc函数设计低通滤波器
   • 应用汉明窗减少频谱泄漏
   • 滤波器长度可调 (当前为101点)

4. 卷积实现带限滤波:

   • 对实部和虚部分别进行卷积
   • 保留复数信号特性

参数调整建议

1. 滤波器长度:

   • 增加长度可获得更锐利的截止特性
   • 减少长度可提高计算速度
   • 修改filter_length变量 (应为奇数)

2. 带宽调整:

   • 修改BW参数设置所需带宽
   • 确保带宽小于采样率的一半 (Nyquist频率)

3. 信号特性:

   • sigma: 控制信号幅度
   • N: 样本点数
   • Fs: 采样率

这种方法虽然计算效率低于FFT方法，但对于中等长度的信号(如4096点)完全可行，且避免了外部库依赖问题。生成的信号仍具有带限高斯白噪声的特性，你学废了吗？

代码同步上传资源~~~~

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---