2. 自适应噪声估计与谱减法降噪
任务说明
本案例实现一个面向嵌入式环境的单通道语音降噪算法：对带噪语音帧进行短时傅里叶变换（STFT），在频域通过递归平滑进行噪声功率谱估计，并采用谱减法计算增益函数，最终通过逆变换恢复增强语音。算法需在给定C代码中完整实现，不依赖外部DSP库，可直接在PC或嵌入式平台编译运行。
算法设计
设带噪语音为 \( y(n) = s(n) + d(n) \)，分帧加窗后做N点FFT得到频域表示 \( Y(k,l) \)，其中k为频点索引，l为帧索引。噪声估计采用最小值统计的递归平滑方法：
\[ P(k,l) = \alpha \cdot P(k,l-1) + (1-\alpha) \cdot |Y(k,l)|^2 \]
\[ P_{min}(k,l) = \min\{P(k,l), P(k,l-1), ..., P(k,l-L+1)\} \]
其中 \(\alpha=0.8\) 为平滑因子，\(L=50\) 为搜索窗长度。增益函数采用谱减规则：
\[ G(k,l) = \max\left( \sqrt{1 - \frac{\beta \cdot P_{min}(k,l)}{|Y(k,l)|^2}}, \delta \right) \]
\(\beta=1.2\) 为过减因子，\(\delta=0.01\) 为谱底限制。增强语音为 \( \hat{S}(k,l) = G(k,l) \cdot Y(k,l) \)，经IFFT和重叠相加得到时域输出。
代码实现
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define FRAME_SIZE 256 // 帧长
#define FFT_SIZE 256 // FFT点数
#define HOP_SIZE 128 // 帧移
#define SMOOTH_ALPHA 0.8f // 平滑因子
#define SEARCH_LEN 50 // 最小值搜索窗
#define OVER_SUB 1.2f // 过减因子
#define FLOOR_VAL 0.01f // 谱底限制
// 简易FFT实现（基2，需FFT_SIZE为2的幂）
void fft(float* re, float* im, int n, int invert) {
for (int i = 1, j = 0; i < n; i++) {
int bit = n >> 1;
for (; j & bit; bit >>= 1) j ^= bit;
j ^= bit;
if (i < j) {
float tr = re[i], ti = im[i];
re[i] = re[j]; im[i] = im[j];
re[j] = tr; im[j] = ti;
}
}
for (int len = 2; len <= n; len <<= 1) {
float ang = 2 * M_PI / len * (invert ? 1 : -1);
float w_re = cosf(ang), w_im = sinf(ang);
for (int i = 0; i < n; i += len) {
float cur_re = 1, cur_im = 0;
for (int j = 0; j < len / 2; j++) {
float u_re = re[i + j], u_im = im[i + j];
float v_re = re[i + j + len / 2] * cur_re - im[i + j + len / 2] * cur_im;
float v_im = re[i + j + len / 2] * cur_im + im[i + j + len / 2] * cur_re;
re[i + j] = u_re + v_re;
im[i + j] = u_im + v_im;
re[i + j + len / 2] = u_re - v_re;
im[i + j + len / 2] = u_im - v_im;
float nxt_re = cur_re * w_re - cur_im * w_im;
cur_im = cur_re * w_im + cur_im * w_re;
cur_re = nxt_re;
}
}
}
if (invert) {
for (int i = 0; i < n; i++) { re[i] /= n; im[i] /= n; }
}
}
// 噪声估计与谱减降噪处理
void noise_reduction(const float* input, float* output, int num_samples) {
int num_frames = (num_samples - FRAME_SIZE) / HOP_SIZE + 1;
float* window = (float*)malloc(FRAME_SIZE * sizeof(float));
float* re = (float*)malloc(FFT_SIZE * sizeof(float));
float* im = (float*)malloc(FFT_SIZE * sizeof(float));
float* pow_spec = (float*)malloc(FFT_SIZE * sizeof(float));
float* smooth_pow = (float*)malloc(FFT_SIZE * sizeof(float));
float* min_pow = (float*)malloc(FFT_SIZE * sizeof(float));
float* buffer = (float*)calloc(num_samples + FRAME_SIZE, sizeof(float));
float* history = (float*)calloc(FFT_SIZE * SEARCH_LEN, sizeof(float));
int* hist_idx = (int*)calloc(FFT_SIZE, sizeof(int));
float* hist_min = (float*)malloc(FFT_SIZE * sizeof(float));
// 初始化
for (int i = 0; i < FFT_SIZE; i++) {
smooth_pow[i] = 0.0f;
min_pow[i] = 1e10f;
hist_min[i] = 1e10f;
}
// 汉宁窗
for (int i = 0; i < FRAME_SIZE; i++)
window[i] = 0.5f * (1 - cosf(2 * M_PI * i / (FRAME_SIZE - 1)));
for (int f = 0; f < num_frames; f++) {
int start = f * HOP_SIZE;
// 分帧加窗
memset(re, 0, FFT_SIZE * sizeof(float));
memset(im, 0, FFT_SIZE * sizeof(float));
for (int i = 0; i < FRAME_SIZE; i++)
re[i] = input[start + i] * window[i];
// FFT
fft(re, im, FFT_SIZE, 0);
// 功率谱与递归平滑
for (int k = 0; k < FFT_SIZE; k++) {
pow_spec[k] = re[k] * re[k] + im[k] * im[k];
smooth_pow[k] = SMOOTH_ALPHA * smooth_pow[k] + (1 - SMOOTH_ALPHA) * pow_spec[k];
}
// 最小值统计噪声估计
for (int k = 0; k < FFT_SIZE; k++) {
int idx = hist_idx[k];
history[k * SEARCH_LEN + idx] = smooth_pow[k];
hist_idx[k] = (idx + 1) % SEARCH_LEN;
// 更新最小值
hist_min[k] = history[k * SEARCH_LEN];
for (int j = 1; j < SEARCH_LEN; j++) {
if (history[k * SEARCH_LEN + j] < hist_min[k])
hist_min[k] = history[k * SEARCH_LEN + j];
}
min_pow[k] = hist_min[k];
}
// 谱减增益
for (int k = 0; k < FFT_SIZE; k++) {
float gain;
if (pow_spec[k] > 1e-10f) {
float ratio = OVER_SUB * min_pow[k] / pow_spec[k];
gain = sqrtf(fmaxf(1.0f - ratio, 0.0f));
gain = fmaxf(gain, FLOOR_VAL);
} else {
gain = FLOOR_VAL;
}
re[k] *= gain;
im[k] *= gain;
}
// IFFT
fft(re, im, FFT_SIZE, 1);
// 重叠相加
for (int i = 0; i < FRAME_SIZE; i++)
buffer[start + i] += re[i] * window[i];
}
// 输出（去除延迟）
memcpy(output, buffer, num_samples * sizeof(float));
free(window); free(re); free(im); free(pow_spec);
free(smooth_pow); free(min_pow); free(buffer);
free(history); free(hist_idx); free(hist_min);
}
// 测试主函数
int main() {
// 生成测试信号：1kHz正弦 + 白噪声
int num_samples = 16000; // 1秒@16kHz
float* noisy = (float*)malloc(num_samples * sizeof(float));
float* clean = (float*)malloc(num_samples * sizeof(float));
srand(42);
for (int i = 0; i < num_samples; i++) {
float tone = 0.5f * sinf(2 * M_PI * 1000 * i / 16000.0f);
float noise = 0.2f * ((float)rand() / RAND_MAX * 2 - 1);
noisy[i] = tone + noise;
clean[i] = tone;
}
float* enhanced = (float*)malloc(num_samples * sizeof(float));
noise_reduction(noisy, enhanced, num_samples);
// 计算降噪前后信噪比（分段估计）
float noise_pow_before = 0, noise_pow_after = 0, sig_pow = 0;
for (int i = 0; i < num_samples; i++) {
float err_before = noisy[i] - clean[i];
float err_after = enhanced[i] - clean[i];
noise_pow_before += err_before * err_before;
noise_pow_after += err_after * err_after;
sig_pow += clean[i] * clean[i];
}
float snr_before = 10 * log10f(sig_pow / noise_pow_before);
float snr_after = 10 * log10f(sig_pow / noise_pow_after);
printf("输入SNR: %.2f dB\n", snr_before);
printf("输出SNR: %.2f dB\n", snr_after);
printf("SNR提升: %.2f dB\n", snr_after - snr_before);
free(noisy); free(clean); free(enhanced);
return 0;
}
编译与运行
gcc -o noise_reduction noise_reduction.c -lm
./noise_reduction
验证方法
程序生成1kHz正弦信号叠加均匀白噪声的测试数据，调用降噪函数处理后，计算处理前后信噪比。判定标准：输出SNR应高于输入SNR，且SNR提升量在5～15dB范围内（取决于噪声随机种子）。若输出SNR低于输入SNR，说明噪声估计或谱减参数存在问题，需检查平滑因子、搜索窗长度或过减因子设置。