使用 fopen_s读取一个wav文件并对其求WVD并输出WVD结果不使用库函数C语言并解释代码

以下是使用 C 语言读取一个 WAV 文件并对其求 WVD（Wigner Ville Distribution，即维格纳分布）的代码，并输出 WVD 结果：

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

#define PI 3.14159265358979323846

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

// 计算 WVD
void wvd(complex *x, int N, double *W) {
    int i, j;
    double t;
    complex *X = (complex *)malloc(N * sizeof(complex));
    complex *Y = (complex *)malloc(N * sizeof(complex));
    complex *Z = (complex *)malloc(N * sizeof(complex));

    // FFT
    for (i = 0; i < N; i++) {
        X[i].real = x[i].real;
        X[i].imag = x[i].imag;
        Y[i].real = 0.0;
        Y[i].imag = 0.0;
        Z[i].real = 0.0;
        Z[i].imag = 0.0;
    }
    fft(X, N, 1); // 做正向 FFT
    for (i = 0; i < N; i++) {
        t = 2.0 * PI * i / N;
        for (j = 0; j < N; j++) {
            Y[j].real = x[j].real * cos(t) + x[j].imag * sin(t);
            Y[j].imag = -x[j].real * sin(t) + x[j].imag * cos(t);
        }
        fft(Y, N, 1); // 做正向 FFT
        for (j = 0; j < N; j++) {
            Z[j].real += X[j].real * Y[j].real - X[j].imag * Y[j].imag;
            Z[j].imag += X[j].real * Y[j].imag + X[j].imag * Y[j].real;
        }
    }
    fft(Z, N, -1); // 做反向 FFT
    for (i = 0; i < N; i++) {
        W[i] = Z[i].real / (double)N;
    }

    free(X);
    free(Y);
    free(Z);
}

// 做 DFT/FFT
void fft(complex *x, int N, int inv) {
    int i, j, k, s, p;
    int *bit = (int *)malloc(N * sizeof(int));
    complex t, u, v;

    for (i = 0; i < N; i++) {
        bit[i] = 0;
        j = i;
        for (k = 0; k < log2(N); k++) {
            bit[i] <<= 1;
            bit[i] |= (j & 1);
            j >>= 1;
        }
    }
    for (i = 0; i < N; i++) {
        if (i < bit[i]) {
            t.real = x[i].real;
            t.imag = x[i].imag;
            x[i].real = x[bit[i]].real;
            x[i].imag = x[bit[i]].imag;
            x[bit[i]].real = t.real;
            x[bit[i]].imag = t.imag;
        }
    }

    for (s = 1; s <= log2(N); s++) {
        p = pow(2, s);
        for (k = 0; k < N; k += p) {
            for (j = 0; j < p / 2; j++) {
                t.real = cos(-2 * PI * j / p);
                t.imag = sin(-2 * PI * j / p);
                u.real = x[k + j].real;
                u.imag = x[k + j].imag;
                v.real = t.real * x[k + j + p / 2].real - t.imag * x[k + j + p / 2].imag;
                v.imag = t.real * x[k + j + p / 2].imag + t.imag * x[k + j + p / 2].real;
                x[k + j].real = u.real + v.real;
                x[k + j].imag = u.imag + v.imag;
                x[k + j + p / 2].real = u.real - v.real;
                x[k + j + p / 2].imag = u.imag - v.imag;
            }
        }
    }

    if (inv == -1) { // 反向 FFT
        for (i = 0; i < N; i++) {
            x[i].real /= (double)N;
            x[i].imag /= (double)N;
        }
    }

    free(bit);
}

int main() {
    FILE *fp;
    char filename[] = "test.wav";
    char riff[4];
    unsigned int chunk_size, subchunk1_size, subchunk2_size;
    unsigned short audio_format, num_channels, bits_per_sample;
    unsigned int sample_rate, byte_rate;
    unsigned int i, j, k, N;
    double duration, *data, *W;

    // 打开 WAV 文件
    if (fopen_s(&fp, filename, "rb") != 0) {
        printf("Cannot open file %s!\n", filename);
        return 1;
    }

    // 读取 WAV 文件头
    fread(riff, sizeof(char), 4, fp);
    fread(&chunk_size, sizeof(unsigned int), 1, fp);
    fread(riff, sizeof(char), 4, fp);
    fread(riff, sizeof(char), 4, fp);
    fread(riff, sizeof(char), 4, fp);
    fread(&subchunk1_size, sizeof(unsigned int), 1, fp);
    fread(&audio_format, sizeof(unsigned short), 1, fp);
    fread(&num_channels, sizeof(unsigned short), 1, fp);
    fread(&sample_rate, sizeof(unsigned int), 1, fp);
    fread(&byte_rate, sizeof(unsigned int), 1, fp);
    fread(riff, sizeof(char), 2, fp);
    fread(&bits_per_sample, sizeof(unsigned short), 1, fp);
    fread(riff, sizeof(char), 4, fp);
    fread(&subchunk2_size, sizeof(unsigned int), 1, fp);

    // 计算采样点数和采样时间
    N = subchunk2_size * 8 / bits_per_sample / num_channels;
    duration = (double)N / sample_rate;

    // 读取数据
    data = (double *)malloc(N * sizeof(double));
    for (i = 0; i < N; i++) {
        if (bits_per_sample == 8) {
            data[i] = (double)(fgetc(fp) - 128);
        } else if (bits_per_sample == 16) {
            short s;
            fread(&s, sizeof(short), 1, fp);
            data[i] = (double)s;
        } else {
            printf("Unsupported bit depth %d!\n", bits_per_sample);
            return 1;
        }
        if (num_channels == 2) {
            if (bits_per_sample == 8) {
                fgetc(fp);
            } else if (bits_per_sample == 16) {
                short s;
                fread(&s, sizeof(short), 1, fp);
            }
        }
    }

    // 关闭 WAV 文件
    fclose(fp);

    // 计算 WVD
    W = (double *)malloc(N * sizeof(double));
    complex *x = (complex *)malloc(N * sizeof(complex));
    for (i = 0; i < N; i++) {
        x[i].real = data[i];
        x[i].imag = 0.0;
    }
    wvd(x, N, W);

    // 输出 WVD 结果
    fp = fopen("wvd.txt", "wt");
    for (i = 0; i < N; i++) {
        fprintf(fp, "%.8lf\n", W[i]);
    }
    fclose(fp);

    free(data);
    free(W);
    free(x);

    return 0;
}

该代码首先打开指定的 WAV 文件，读取 WAV 文件头并计算采样点数和采样时间。然后，它读取 WAV 文件中的数据并存储在一个 double 类型的数组中。接下来，它将数据转换为复数形式并调用 wvd 函数计算 WVD。最后，它将 WVD 结果输出到一个文本文件中。

其中，wvd 函数是计算 WVD 的核心函数。它首先对输入的复数数组进行 FFT 变换，然后将每个频率分量进行位移和幅度调整，最后再对变换后的结果进行反向 FFT 变换。这个过程可以用一个三重循环实现，其中第一重循环遍历所有频率分量，第二重循环遍历所有时间点，第三重循环用于进行 DFT/FFT 变换。由于 FFT 变换的复杂度为 O(N log N)，因此 wvd 函数的时间复杂度为 O(N^2 log N)。