c语言对wav语音谱减法降噪例子

下面是一个简单的C语言例子，演示如何使用谱减法对WAV语音进行降噪。

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

#define PI 3.14159265358979323846

// WAV文件头结构体
typedef struct wav_header {
    char riff[4];            // RIFF标志
    unsigned int size;       // 文件大小
    char wave[4];            // WAVE标志
    char fmt[4];             // fmt标志
    unsigned int fmt_size;   // fmt块大小
    unsigned short format;   // 格式类型
    unsigned short channels; // 声道数
    unsigned int sample_rate;// 采样率
    unsigned int byte_rate;  // 每秒传输字节数
    unsigned short block_align; // 数据块对齐
    unsigned short bits_per_sample; // 每个采样的位数
    char data[4];            // 数据标志
    unsigned int data_size;  // 数据大小
} wav_header;

// 加载WAV文件
int load_wav(char* filename, short** data, int* size, int* channels, int* sample_rate) {
    FILE* fp = fopen(filename, "rb");
    if (!fp) {
        printf("Error: cannot open file '%s'!\n", filename);
        return -1;
    }

    wav_header header;
    fread(&header, sizeof(wav_header), 1, fp);

    if (header.fmt_size != 16) {
        printf("Error: fmt size is not 16!\n");
        fclose(fp);
        return -1;
    }

    if (header.format != 1) {
        printf("Error: format is not PCM!\n");
        fclose(fp);
        return -1;
    }

    if (header.bits_per_sample != 16) {
        printf("Error: bits per sample is not 16!\n");
        fclose(fp);
        return -1;
    }

    if (header.data_size % 2 != 0) {
        printf("Error: data size is not even!\n");
        fclose(fp);
        return -1;
    }

    int num_samples = header.data_size / 2;
    short* samples = (short*)malloc(num_samples * sizeof(short));
    fread(samples, sizeof(short), num_samples, fp);

    *data = samples;
    *size = num_samples;
    *channels = header.channels;
    *sample_rate = header.sample_rate;

    fclose(fp);
    return 0;
}

// 保存WAV文件
int save_wav(char* filename, short* data, int size, int channels, int sample_rate) {
    FILE* fp = fopen(filename, "wb");
    if (!fp) {
        printf("Error: cannot create file '%s'!\n", filename);
        return -1;
    }

    wav_header header;
    header.riff[0] = 'R';
    header.riff[1] = 'I';
    header.riff[2] = 'F';
    header.riff[3] = 'F';
    header.size = size * channels * 2 + sizeof(wav_header) - 8;
    header.wave[0] = 'W';
    header.wave[1] = 'A';
    header.wave[2] = 'V';
    header.wave[3] = 'E';
    header.fmt[0] = 'f';
    header.fmt[1] = 'm';
    header.fmt[2] = 't';
    header.fmt[3] = ' ';
    header.fmt_size = 16;
    header.format = 1;
    header.channels = channels;
    header.sample_rate = sample_rate;
    header.byte_rate = sample_rate * channels * 2;
    header.block_align = channels * 2;
    header.bits_per_sample = 16;
    header.data[0] = 'd';
    header.data[1] = 'a';
    header.data[2] = 't';
    header.data[3] = 'a';
    header.data_size = size * channels * 2;

    fwrite(&header, sizeof(wav_header), 1, fp);
    fwrite(data, sizeof(short), size * channels, fp);

    fclose(fp);
    return 0;
}

// 计算FFT
void fft(short* data, int size) {
    int i, j, k, n;
    double xr, xi, tr, ti, ur, ui, sr, si;

    n = size / 2;
    j = n / 2;

    for (i = 1; i < n - 1; i++) {
        if (i < j) {
            tr = data[j * 2];
            data[j * 2] = data[i * 2];
            data[i * 2] = tr;
            ti = data[j * 2 + 1];
            data[j * 2 + 1] = data[i * 2 + 1];
            data[i * 2 + 1] = ti;
        }

        k = n;

        while (k <= j) {
            j -= k;
            k /= 2;
        }

        j += k;
    }

    for (i = 0; i < log2(n); i++) {
        k = 1 << i;
        ur = 1;
        ui = 0;

        sr = cos(PI / k);
        si = -sin(PI / k);

        for (j = 0; j < k; j++) {
            for (n = j; n < size; n += 2 * k) {
                xr = data[n];
                xi = data[n + 1];
                tr = xr * ur - xi * ui;
                ti = xr * ui + xi * ur;
                data[n] = xr - tr;
                data[n + 1] = xi - ti;
                data[n + k * 2] = tr;
                data[n + k * 2 + 1] = ti;
            }

            xr = ur;
            ur = xr * sr - ui * si;
            ui = xr * si + ui * sr;
        }
    }
}

// 谱减法降噪
void spectral_subtraction(short* data, int size, int sample_rate) {
    int i, j, k, n;
    double* magnitude = (double*)malloc(size * sizeof(double));
    double* noise = (double*)malloc(size * sizeof(double));
    double* power = (double*)malloc(size * sizeof(double));
    double* threshold = (double*)malloc(size * sizeof(double));

    for (i = 0; i < size; i++) {
        magnitude[i] = sqrt(data[i * 2] * data[i * 2] + data[i * 2 + 1] * data[i * 2 + 1]);
        power[i] = magnitude[i] * magnitude[i];
    }

    fft(data, size);

    for (i = 0; i < size; i++) {
        magnitude[i] = sqrt(data[i * 2] * data[i * 2] + data[i * 2 + 1] * data[i * 2 + 1]);
        power[i] = magnitude[i] * magnitude[i];
    }

    for (i = 0; i < size; i++) {
        noise[i] = power[i];
    }

    for (j = 1; j < 10; j++) {
        fft(noise, size);

        for (i = 0; i < size; i++) {
            magnitude[i] = sqrt(noise[i * 2] * noise[i * 2] + noise[i * 2 + 1] * noise[i * 2 + 1]);
            power[i] = magnitude[i] * magnitude[i];
        }

        for (i = 0; i < size; i++) {
            noise[i] = power[i];
        }
    }

    for (i = 0; i < size; i++) {
        threshold[i] = noise[i] * 2;
    }

    for (i = 0; i < size; i++) {
        if (power[i] < threshold[i]) {
            power[i] = 0;
        }
    }

    for (i = 0; i < size; i++) {
        magnitude[i] = sqrt(power[i]);
    }

    for (i = 0; i < size; i++) {
        data[i * 2] = magnitude[i] * cos(atan2(data[i * 2 + 1], data[i * 2]));
        data[i * 2 + 1] = magnitude[i] * sin(atan2(data[i * 2 + 1], data[i * 2]));
    }

    fft(data, size);

    for (i = 0; i < size; i++) {
        data[i * 2] /= size;
        data[i * 2 + 1] /= size;
    }

    free(magnitude);
    free(noise);
    free(power);
    free(threshold);
}

int main(int argc, char* argv[]) {
    if (argc < 3) {
        printf("Usage: %s <input_file> <output_file>\n", argv[0]);
        return -1;
    }

    char* input_file = argv[1];
    char* output_file = argv[2];

    short* data;
    int size, channels, sample_rate;

    if (load_wav(input_file, &data, &size, &channels, &sample_rate) != 0) {
        return -1;
    }

    spectral_subtraction(data, size, sample_rate);

    if (save_wav(output_file, data, size, channels, sample_rate) != 0) {
        return -1;
    }

    free(data);

    return 0;
}

这个例子中，我们首先定义了一个WAV文件头结构体，并定义了一个`load_wav`函数和一个`save_wav`函数，用来加载和保存WAV文件。然后我们定义了一个`fft`函数，用来计算FFT。最后，我们定义了一个`spectral_subtraction`函数，用来进行谱减法降噪。

在`spectral_subtraction`函数中，我们首先计算音频信号的功率谱，并进行FFT变换。然后我们计算噪声信号的功率谱，通过多次FFT变换得到一个比较准确的噪声功率谱。接着，我们计算一个阈值，如果信号功率谱低于阈值，则把信号功率谱设为0。最后，我们根据降噪后的信号功率谱计算出降噪后的信号，并进行IFFT变换，得到降噪后的音频信号。

使用该例子时，需要指定输入WAV文件和输出WAV文件的文件名作为程序的两个参数。