使用 fopen_s读取一个wav文件并对其进行fft和STFT并输出结果c语言代码不使用fftw库函数并解释代码

以下是使用 C 语言实现读取 WAV 文件并对其进行 FFT 和 STFT 的示例代码：

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

#define BUFFER_SIZE 1024
#define PI 3.14159265358979323846

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

// 读取 WAV 文件
int read_wav_file(const char *filename, short **samples, int *sample_rate, int *num_samples) {
    FILE *fp;
    int header_size, bytes_per_sample, num_channels, sample_rate_, num_samples_;
    short *samples_;

    fp = fopen(filename, "rb");
    if (!fp) {
        return -1;
    }

    // 读取 WAV 头信息
    fseek(fp, 22, SEEK_SET);
    fread(&num_channels, sizeof(short), 1, fp);
    fseek(fp, 34, SEEK_SET);
    fread(&sample_rate_, sizeof(int), 1, fp);
    fseek(fp, 40, SEEK_SET);
    fread(&bytes_per_sample, sizeof(short), 1, fp);
    fseek(fp, 44, SEEK_SET);
    fread(&num_samples_, sizeof(int), 1, fp);
    header_size = 44 + bytes_per_sample * num_samples_ - 8;

    // 分配空间并读取数据
    samples_ = (short *)malloc(num_samples_ * sizeof(short));
    fseek(fp, header_size, SEEK_SET);
    fread(samples_, sizeof(short), num_samples_, fp);

    fclose(fp);

    *samples = samples_;
    *sample_rate = sample_rate_;
    *num_samples = num_samples_;

    return 0;
}

// 计算 FFT
void fft(complex *x, int n) {
    int i, j, k;
    int m, r;
    complex t, u;

    for (i = 0, j = 1; j < n - 1; j++) {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) {
            t = x[j];
            x[j] = x[i];
            x[i] = t;
        }
    }

    for (m = 2; m <= n; m <<= 1) {
        r = m >> 1;
        for (i = 0; i < n; i += m) {
            for (j = i, k = 0; j < i + r; j++, k += n / r) {
                t.real = x[j + r].real * cos(-2 * PI * k / n) + x[j + r].imag * sin(-2 * PI * k / n);
                t.imag = x[j + r].imag * cos(-2 * PI * k / n) - x[j + r].real * sin(-2 * PI * k / n);
                u.real = x[j].real;
                u.imag = x[j].imag;
                x[j].real = u.real + t.real;
                x[j].imag = u.imag + t.imag;
                x[j + r].real = u.real - t.real;
                x[j + r].imag = u.imag - t.imag;
            }
        }
    }
}

// 计算 STFT
void stft(short *samples, int sample_rate, int num_samples, int window_size, int step, int fft_size, double ***stft_result, int *num_frames) {
    int i, j, k;
    int num_frames_, num_windows, num_bins;
    double *window, *windowed_samples;
    complex *x;
    double **stft_result_;

    // 计算窗函数
    window = (double *)malloc(window_size * sizeof(double));
    for (i = 0; i < window_size; i++) {
        window[i] = 0.5 - 0.5 * cos(2 * PI * i / (window_size - 1));
    }

    // 计算 STFT
    num_frames_ = (num_samples - window_size) / step + 1;
    num_windows = window_size / 2 + 1;
    num_bins = fft_size / 2 + 1;
    stft_result_ = (double **)malloc(num_windows * sizeof(double *));
    x = (complex *)malloc(fft_size * sizeof(complex));
    windowed_samples = (double *)malloc(window_size * sizeof(double));
    for (i = 0; i < num_windows; i++) {
        stft_result_[i] = (double *)malloc(num_bins * num_frames_ * sizeof(double));
    }
    for (i = 0; i < num_samples - window_size; i += step) {
        // 加窗
        for (j = 0; j < window_size; j++) {
            windowed_samples[j] = (double)samples[i + j] * window[j];
        }
        // 补零
        for (j = 0; j < fft_size - window_size; j++) {
            windowed_samples[j + window_size] = 0;
        }
        // 计算 FFT
        for (j = 0; j < fft_size; j++) {
            x[j].real = windowed_samples[j];
            x[j].imag = 0;
        }
        fft(x, fft_size);
        // 计算 STFT
        for (j = 0; j < num_windows; j++) {
            for (k = 0; k < num_bins; k++) {
                stft_result_[j][i / step * num_bins + k] = sqrt(x[j * (num_bins - 1) / (num_windows - 1)].real * x[j * (num_bins - 1) / (num_windows - 1)].real + x[j * (num_bins - 1) / (num_windows - 1)].imag * x[j * (num_bins - 1) / (num_windows - 1)].imag);
            }
        }
    }

    free(window);
    free(windowed_samples);
    free(x);

    *stft_result = stft_result_;
    *num_frames = num_frames_;
}

int main() {
    short *samples;
    int sample_rate, num_samples;
    int window_size, step, fft_size;
    double **stft_result;
    int num_frames;
    int i, j;

    // 读取 WAV 文件
    if (read_wav_file("test.wav", &samples, &sample_rate, &num_samples)) {
        printf("Failed to read WAV file.\n");
        return -1;
    }

    // 设置参数
    window_size = 1024;
    step = 512;
    fft_size = 1024;

    // 计算 STFT
    stft(samples, sample_rate, num_samples, window_size, step, fft_size, &stft_result, &num_frames);

    // 输出结果
    for (i = 0; i < num_frames; i++) {
        for (j = 0; j < fft_size / 2 + 1; j++) {
            printf("%f ", stft_result[j][i]);
        }
        printf("\n");
    }

    // 释放内存
    for (i = 0; i < fft_size / 2 + 1; i++) {
        free(stft_result[i]);
    }
    free(stft_result);
    free(samples);

    return 0;
}

上述代码中，`read_wav_file` 函数用于读取 WAV 文件，返回值为 0 表示读取成功，-1 表示读取失败。`fft` 函数用于计算 FFT，`stft` 函数用于计算 STFT。在 `stft` 函数中，首先计算窗函数，然后对每个窗口进行加窗、补零、FFT 和 STFT 计算。计算 STFT 时使用的是平方和谐平均值（Root Mean Square）。

需要注意的是，在本代码中，FFT 的实现采用了 Cooley-Tukey 算法。由于 FFT 算法的复杂度为 O(N log N)，其中 N 为采样点数，因此 FFT 的计算速度非常快。而且，由于 FFT 的特殊性质，可以将 STFT 计算中的 FFT 部分高度优化，进一步提高计算速度。