基于谱减法的自适应消嗓算法 C语言

下面是基于谱减法的自适应消噪算法的C语言实现代码：

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

#define N 256           //每帧长度
#define PI 3.1415926535 //圆周率

//复数结构体
typedef struct
{
    double real;
    double imag;
} complex;

//分帧函数
void framing(double *signal, double **frame, int frame_num)
{
    int i, j, k;
    for (i = 0; i < frame_num; i++)
    {
        for (j = 0; j < N; j++)
        {
            k = i * N + j;
            frame[i][j] = signal[k];
        }
    }
}

//FFT变换函数
void fft(complex *x, complex *y, int n)
{
    int i, j, k;
    complex t, u;
    double p;

    for (i = 0; i < n; i++)
    {
        y[i].real = x[i].real;
        y[i].imag = x[i].imag;
    }

    for (i = 0; i < n; i++)
    {
        j = 0;
        k = i;
        while (k > 0)
        {
            j = (j << 1) | (k & 1);
            k >>= 1;
        }
        if (i < j)
        {
            t.real = y[i].real;
            t.imag = y[i].imag;
            y[i].real = y[j].real;
            y[i].imag = y[j].imag;
            y[j].real = t.real;
            y[j].imag = t.imag;
        }
    }

    for (i = 2; i <= n; i <<= 1)
    {
        p = PI / (i >> 1);
        t.real = cos(p);
        t.imag = -sin(p);
        for (j = 0; j < n; j += i)
        {
            u.real = 1.0;
            u.imag = 0.0;
            for (k = j; k < j + (i >> 1); k++)
            {
                t.real = u.real * y[k + (i >> 1)].real - u.imag * y[k + (i >> 1)].imag;
                t.imag = u.real * y[k + (i >> 1)].imag + u.imag * y[k + (i >> 1)].real;
                y[k + (i >> 1)].real = y[k].real - t.real;
                y[k + (i >> 1)].imag = y[k].imag - t.imag;
                y[k].real += t.real;
                y[k].imag += t.imag;
                t.real = u.real * cos(p) - u.imag * sin(p);
                t.imag = u.real * sin(p) + u.imag * cos(p);
                u.real = t.real;
                u.imag = t.imag;
            }
        }
    }
}

//功率谱函数
void power_spectrum(complex *y, double *power, int n)
{
    int i;
    for (i = 0; i < n; i++)
        power[i] = y[i].real * y[i].real + y[i].imag * y[i].imag;
}

//谱减函数
void spectral_subtraction(double *power, double *noise_power, double *factor, int n)
{
    int i;
    for (i = 0; i < n; i++)
    {
        if (power[i] < noise_power[i])
            factor[i] = 1.0;
        else
            factor[i] = sqrt((power[i] - noise_power[i]) / power[i]);
    }
}

//IFFT变换函数
void ifft(complex *y, complex *x, int n)
{
    int i;
    for (i = 0; i < n; i++)
        y[i].imag = -y[i].imag;
    fft(y, x, n);
    for (i = 0; i < n; i++)
    {
        x[i].real = x[i].real / n;
        x[i].imag = -x[i].imag / n;
    }
}

//重叠相加函数
void overlap_add(double **frame, double *result, int frame_num)
{
    int i, j, k;
    for (i = 0; i < frame_num; i++)
    {
        for (j = 0; j < N; j++)
        {
            k = i * N + j;
            result[k] += frame[i][j];
        }
    }
}

//消噪函数
void noise_reduction(double *signal, double *result, int signal_len, double *noise_power, double threshold)
{
    int frame_num = (int)ceil(signal_len * 1.0 / N); //总帧数
    double **frame = (double **)malloc(frame_num * sizeof(double *)); //每帧信号数组
    int i, j;

    for (i = 0; i < frame_num; i++)
        frame[i] = (double *)malloc(N * sizeof(double));

    framing(signal, frame, frame_num);

    complex *x = (complex *)malloc(N * sizeof(complex));
    complex *y = (complex *)malloc(N * sizeof(complex));
    double *power = (double *)malloc(N * sizeof(double));
    double *factor = (double *)malloc(N * sizeof(double));

    for (i = 0; i < N; i++)
    {
        x[i].real = 0.0;
        x[i].imag = 0.0;
    }

    for (i = 0; i < frame_num; i++)
    {
        for (j = 0; j < N; j++)
        {
            if (i * N + j < signal_len)
                x[j].real = frame[i][j];
            else
                x[j].real = 0.0;
            x[j].imag = 0.0;
        }

        fft(x, y, N);

        power_spectrum(y, power, N);

        spectral_subtraction(power, noise_power, factor, N);

        for (j = 0; j < N; j++)
        {
            y[j].real = y[j].real * factor[j];
            y[j].imag = y[j].imag * factor[j];
        }

        ifft(y, x, N);

        for (j = 0; j < N; j++)
        {
            if (i * N + j < signal_len)
                frame[i][j] = x[j].real;
            else
                frame[i][j] = 0.0;
        }
    }

    overlap_add(frame, result, frame_num);

    for (i = 0; i < frame_num; i++)
        free(frame[i]);
    free(frame);

    free(x);
    free(y);
    free(power);
    free(factor);
}

int main()
{
    double signal[1024]; //原始信号
    double result[1024]; //消噪后的信号
    double noise_power[N]; //噪声功率谱
    double threshold = 0.5; //谱减因子阈值
    int signal_len = 1024; //原始信号长度
    int i;

    for (i = 0; i < signal_len; i++)
        signal[i] = sin(2.0 * PI * 100.0 * i / 8000.0) + 0.5 * sin(2.0 * PI * 300.0 * i / 8000.0); //原始信号为100Hz和300Hz的正弦信号叠加

    for (i = 0; i < N; i++)
        noise_power[i] = 0.0001; //初始化噪声功率谱，设一个很小的值

    noise_reduction(signal, result, signal_len, noise_power, threshold);

    for (i = 0; i < signal_len; i++)
        printf("%f ", result[i]);

    return 0;
}

这份代码实现了基于谱减法的自适应消噪算法，使用了分帧、FFT、功率谱、谱减、IFFT和重叠相加等技术。需要注意的是，为了简化代码，本例中噪声功率谱被初始化为一个很小的值，实际应用中需要根据实际情况进行估计。另外，谱减因子阈值也需要根据实际情况进行调整，以达到较好的消噪效果。