音频如何低延时回声消除与降噪

// 定义数学常量宏，确保M_PI可用
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
 
// 手动定义M_PI作为备选，防止某些编译器不支持
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
// 配置参数
#define FRAME_SIZE     512     // 帧大小
#define FFT_SIZE       512     // FFT大小(必须为2的幂)
#define OVERLAP        256     // 重叠大小(通常为帧大小的一半)
#define NOISE_FRAMES   10      // 用于估计噪声的帧数
#define ALPHA          0.9     // 噪声更新系数
#define BETA           0.001   // 谱减过度保护系数
 
// 复数结构体
typedef struct {
    float real;
    float imag;
} Complex;
 
// 全局变量
float noise_power[FFT_SIZE/2 + 1];  // 噪声功率谱
 
// FFT相关函数声明
void fft(Complex *x, int n, int invert);
void multiply_complex(Complex *a, Complex *b, Complex *result);
void add_complex(Complex *a, Complex *b, Complex *result);
void subtract_complex(Complex *a, Complex *b, Complex *result);
float complex_abs(Complex *c);
 
// 信号处理函数声明
void hamming_window(float *frame, int size);
void estimate_noise(float *signal, int signal_len);
void spectral_subtraction(float *signal, float *output, int signal_len);
 
int main() {
    // 示例：生成带噪声的信号(1000Hz正弦波 + 高斯噪声)
    int sample_rate = 16000;    // 采样率
    int duration = 2;           // 信号时长(秒)
    int signal_len = sample_rate * duration;
    float *signal = (float*)malloc(signal_len * sizeof(float));
    float *output = (float*)malloc(signal_len * sizeof(float));
    
    // 生成带噪声信号
    for (int i = 0; i < signal_len; i++) {
        // 1000Hz正弦信号
        float clean = sin(2 * M_PI * 1000 * i / sample_rate);
        // 高斯噪声(均值0，方差0.1)
        float noise = ((float)rand() / RAND_MAX - 0.5) * 0.6;
        signal[i] = clean + noise;
    }
    
    printf("开始谱减法降噪处理...\n");
    
    // 估计噪声(假设信号开始部分是纯噪声)
    estimate_noise(signal, signal_len);
    
    // 应用谱减法降噪
    spectral_subtraction(signal, output, signal_len);
    
    printf("降噪处理完成！\n");
    
    // 保存结果到文件(可用于后续分析或播放)
    FILE *fp_in = fopen("noisy_signal.raw", "wb");
    FILE *fp_out = fopen("denoised_signal.raw", "wb");
    fwrite(signal, sizeof(float), signal_len, fp_in);
    fwrite(output, sizeof(float), signal_len, fp_out);
    fclose(fp_in);
    fclose(fp_out);
    
    printf("带噪信号已保存到noisy_signal.raw\n");
    printf("降噪后信号已保存到denoised_signal.raw\n");
    
    // 释放内存
    free(signal);
    free(output);
    
    return 0;
}
 
/**
 * 汉宁窗函数，减少频谱泄漏
 */
void hamming_window(float *frame, int size) {
    for (int i = 0; i < size; i++) {
        frame[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (size - 1));
    }
}
 
/**
 * 估计噪声功率谱
 * 假设信号开始的NOISE_FRAMES帧是纯噪声
 */
void estimate_noise(float *signal, int signal_len) {
    Complex fft_result[FFT_SIZE];
    float frame[FRAME_SIZE];
    
    // 初始化噪声功率谱为0
    memset(noise_power, 0, sizeof(noise_power));
    
    // 取前NOISE_FRAMES帧估计噪声
    int frames_used = 0;
    for (int i = 0; i < signal_len && frames_used < NOISE_FRAMES; 
         i += FRAME_SIZE - OVERLAP, frames_used++) {
        
        // 提取一帧信号
        memset(frame, 0, FRAME_SIZE * sizeof(float));
        int copy_len = (i + FRAME_SIZE <= signal_len) ? FRAME_SIZE : signal_len - i;
        memcpy(frame, &signal[i], copy_len * sizeof(float));
        
        // 加汉宁窗
        hamming_window(frame, FRAME_SIZE);
        
        // 转换为复数并执行FFT
        for (int j = 0; j < FFT_SIZE; j++) {
            fft_result[j].real = (j < FRAME_SIZE) ? frame[j] : 0;
            fft_result[j].imag = 0;
        }
        fft(fft_result, FFT_SIZE, 0);
        
        // 计算功率谱并累加
        for (int j = 0; j <= FFT_SIZE/2; j++) {
            float power = fft_result[j].real * fft_result[j].real + 
                         fft_result[j].imag * fft_result[j].imag;
            noise_power[j] += power;
        }
    }
    
    // 求平均值
    for (int j = 0; j <= FFT_SIZE/2; j++) {
        noise_power[j] /= frames_used;
    }
}
 
/**
 * 谱减法降噪主函数
 */
void spectral_subtraction(float *signal, float *output, int signal_len) {
    Complex fft_result[FFT_SIZE];
    float frame[FRAME_SIZE];
    float overlap_buf[OVERLAP] = {0};  // 重叠缓冲区，用于帧拼接
    
    // 初始化输出
    memset(output, 0, signal_len * sizeof(float));
    
    // 分帧处理
    for (int i = 0; i < signal_len; i += FRAME_SIZE - OVERLAP) {
        // 提取一帧信号
        memset(frame, 0, FRAME_SIZE * sizeof(float));
        int copy_len = (i + FRAME_SIZE <= signal_len) ? FRAME_SIZE : signal_len - i;
        memcpy(frame, &signal[i], copy_len * sizeof(float));
        
        // 加汉宁窗
        hamming_window(frame, FRAME_SIZE);
        
        // 转换为复数并执行FFT
        for (int j = 0; j < FFT_SIZE; j++) {
            fft_result[j].real = (j < FRAME_SIZE) ? frame[j] : 0;
            fft_result[j].imag = 0;
        }
        fft(fft_result, FFT_SIZE, 0);
        
        // 谱减处理
        for (int j = 0; j <= FFT_SIZE/2; j++) {
            // 计算信号功率谱
            float signal_power = fft_result[j].real * fft_result[j].real + 
                                fft_result[j].imag * fft_result[j].imag;
            
            // 谱减 (加入保护系数防止过度减法)
            float subt_result = signal_power - ALPHA * noise_power[j];
            if (subt_result < BETA * noise_power[j]) {
                subt_result = BETA * noise_power[j];
            }
            
            // 保持相位不变，幅度为谱减结果的平方根
            float amp = sqrt(subt_result);
            float phase = atan2(fft_result[j].imag, fft_result[j].real);
            
            // 更新复数频谱
            fft_result[j].real = amp * cos(phase);
            fft_result[j].imag = amp * sin(phase);
            
            // 处理共轭对称部分(除了直流和Nyquist分量)
            if (j != 0 && j != FFT_SIZE/2) {
                fft_result[FFT_SIZE - j].real = fft_result[j].real;
                fft_result[FFT_SIZE - j].imag = -fft_result[j].imag;
            }
        }
        
        // 执行逆FFT
        fft(fft_result, FFT_SIZE, 1);
        
        // 将结果转换回实数信号
        float ifft_result[FRAME_SIZE];
        for (int j = 0; j < FRAME_SIZE; j++) {
            ifft_result[j] = fft_result[j].real / FFT_SIZE;  // 除以FFT大小进行归一化
        }
        
        // 重叠相加
        int output_idx = i;
        for (int j = 0; j < OVERLAP; j++) {
            if (output_idx + j < signal_len) {
                output[output_idx + j] += overlap_buf[j] + ifft_result[j];
            }
        }
        
        for (int j = OVERLAP; j < FRAME_SIZE; j++) {
            if (output_idx + j < signal_len) {
                output[output_idx + j] += ifft_result[j];
            }
        }
        
        // 更新重叠缓冲区
        memcpy(overlap_buf, &ifft_result[FRAME_SIZE - OVERLAP], 
               OVERLAP * sizeof(float));
    }
}
 
/**
 * FFT实现 (Cooley-Tukey算法)
 * x: 输入/输出复数数组
 * n: 变换大小(必须为2的幂)
 * invert: 0=正向FFT, 1=逆向FFT
 */
void fft(Complex *x, 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) {
            Complex temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
    }
    
    // Cooley-Tukey FFT
    for (int len = 2; len <= n; len <<= 1) {
        float ang = 2 * M_PI / len * (invert ? -1 : 1);
        Complex wlen = {cos(ang), sin(ang)};  // 旋转因子
        
        for (int i = 0; i < n; i += len) {
            Complex w = {1, 0};
            for (int j = 0; j < len / 2; j++) {
                Complex u = x[i + j];
                Complex v;
                v.real = x[i + j + len/2].real * w.real - x[i + j + len/2].imag * w.imag;
                v.imag = x[i + j + len/2].real * w.imag + x[i + j + len/2].imag * w.real;
                
                x[i + j].real = u.real + v.real;
                x[i + j].imag = u.imag + v.imag;
                x[i + j + len/2].real = u.real - v.real;
                x[i + j + len/2].imag = u.imag - v.imag;
                
                // 更新旋转因子
                Complex next_w;
                next_w.real = w.real * wlen.real - w.imag * wlen.imag;
                next_w.imag = w.real * wlen.imag + w.imag * wlen.real;
                w = next_w;
            }
        }
    }
    
    // 逆FFT时进行归一化
    if (invert) {
        for (int i = 0; i < n; i++) {
            x[i].real /= n;
            x[i].imag /= n;
        }
    }
}
 
/**
 * 复数乘法
 */
void multiply_complex(Complex *a, Complex *b, Complex *result) {
    result->real = a->real * b->real - a->imag * b->imag;
    result->imag = a->real * b->imag + a->imag * b->real;
}
 
/**
 * 复数加法
 */
void add_complex(Complex *a, Complex *b, Complex *result) {
    result->real = a->real + b->real;
    result->imag = a->imag + b->imag;
}
 
/**
 * 复数减法
 */
void subtract_complex(Complex *a, Complex *b, Complex *result) {
    result->real = a->real - b->real;
    result->imag = a->imag - b->imag;
}
 
/**
 * 计算复数的模
 */
float complex_abs(Complex *c) {
    return sqrt(c->real * c->real + c->imag * c->imag);
}