使用 fopen_s读取一个wav文件并对其求WVD和STFTS并输出WVD和stft的结不使用库函数

读取 WAV 文件并解析其头部信息的代码示例如下：

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

typedef struct {
    char ChunkID[4];
    uint32_t ChunkSize;
    char Format[4];
} RIFF_Chunk;

typedef struct {
    char Subchunk1ID[4];
    uint32_t Subchunk1Size;
    uint16_t AudioFormat;
    uint16_t NumChannels;
    uint32_t SampleRate;
    uint32_t ByteRate;
    uint16_t BlockAlign;
    uint16_t BitsPerSample;
} FMT_Subchunk;

typedef struct {
    char Subchunk2ID[4];
    uint32_t Subchunk2Size;
} DATA_Subchunk;

typedef struct {
    RIFF_Chunk riff;
    FMT_Subchunk fmt;
    DATA_Subchunk data;
} WAV_Header;

WAV_Header read_wav_header(FILE* fp) {
    WAV_Header header;
    fread(&header.riff, sizeof(RIFF_Chunk), 1, fp);
    fread(&header.fmt, sizeof(FMT_Subchunk), 1, fp);
    fread(&header.data, sizeof(DATA_Subchunk), 1, fp);
    return header;
}

int main() {
    FILE* fp;
    fopen_s(&fp, "test.wav", "rb");
    if (fp == NULL) {
        printf("Failed to open file.\n");
        return -1;
    }
    WAV_Header header = read_wav_header(fp);
    printf("Sample rate: %u\n", header.fmt.SampleRate);
    printf("Bits per sample: %u\n", header.fmt.BitsPerSample);
    uint32_t num_samples = header.data.Subchunk2Size / (header.fmt.BitsPerSample / 8);
    printf("Number of samples: %u\n", num_samples);
    fclose(fp);
    return 0;
}

接下来，我们可以使用 FFT 算法来计算 STFT，用 Wigner-Ville 分布来计算 WVD。FFT 的代码可以使用现成的库函数（比如 FFTW），但是 WVD 的计算需要自己编写代码。以下是一个简单的 WVD 计算程序，供参考：

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

#define PI 3.14159265358979323846

typedef struct {
    char ChunkID[4];
    uint32_t ChunkSize;
    char Format[4];
} RIFF_Chunk;

typedef struct {
    char Subchunk1ID[4];
    uint32_t Subchunk1Size;
    uint16_t AudioFormat;
    uint16_t NumChannels;
    uint32_t SampleRate;
    uint32_t ByteRate;
    uint16_t BlockAlign;
    uint16_t BitsPerSample;
} FMT_Subchunk;

typedef struct {
    char Subchunk2ID[4];
    uint32_t Subchunk2Size;
} DATA_Subchunk;

typedef struct {
    RIFF_Chunk riff;
    FMT_Subchunk fmt;
    DATA_Subchunk data;
} WAV_Header;

WAV_Header read_wav_header(FILE* fp) {
    WAV_Header header;
    fread(&header.riff, sizeof(RIFF_Chunk), 1, fp);
    fread(&header.fmt, sizeof(FMT_Subchunk), 1, fp);
    fread(&header.data, sizeof(DATA_Subchunk), 1, fp);
    return header;
}

double* read_wav_data(FILE* fp, uint32_t num_samples) {
    double* data = (double*)malloc(num_samples * sizeof(double));
    if (data == NULL) {
        return NULL;
    }
    WAV_Header header = read_wav_header(fp);
    fseek(fp, sizeof(header), SEEK_SET);
    uint16_t bytes_per_sample = header.fmt.BitsPerSample / 8;
    double max_amplitude = pow(2, header.fmt.BitsPerSample - 1);
    for (uint32_t i = 0; i < num_samples; i++) {
        uint8_t buffer[4];
        fread(buffer, bytes_per_sample, 1, fp);
        int16_t sample = 0;
        for (uint16_t j = 0; j < bytes_per_sample; j++) {
            sample |= buffer[j] << (j * 8);
        }
        data[i] = sample / max_amplitude;
    }
    return data;
}

double* wvd(double* data, uint32_t num_samples, uint32_t window_size, uint32_t hop_size) {
    double* result = (double*)malloc(num_samples * window_size * sizeof(double));
    if (result == NULL) {
        return NULL;
    }
    double* window = (double*)malloc(window_size * sizeof(double));
    if (window == NULL) {
        free(result);
        return NULL;
    }
    for (uint32_t i = 0; i < window_size; i++) {
        window[i] = sin(PI * (i - (window_size - 1) / 2.0) / window_size);
    }
    for (uint32_t t = 0; t < num_samples - window_size; t += hop_size) {
        for (uint32_t i = 0; i < window_size; i++) {
            for (uint32_t j = 0; j < window_size; j++) {
                uint32_t index = t * window_size + i * window_size + j;
                double real = 0.0;
                double imag = 0.0;
                for (uint32_t k = 0; k < window_size; k++) {
                    double x = data[t + k] * window[k];
                    double y = data[t + i + k] * window[k];
                    double phase = 2 * PI * j * k / window_size;
                    double cos_phase = cos(phase);
                    double sin_phase = sin(phase);
                    real += x * y * cos_phase;
                    imag += x * y * sin_phase;
                }
                result[index] = real * real + imag * imag;
            }
        }
    }
    free(window);
    return result;
}

int main() {
    FILE* fp;
    fopen_s(&fp, "test.wav", "rb");
    if (fp == NULL) {
        printf("Failed to open file.\n");
        return -1;
    }
    WAV_Header header = read_wav_header(fp);
    uint32_t num_samples = header.data.Subchunk2Size / (header.fmt.BitsPerSample / 8);
    double* data = read_wav_data(fp, num_samples);
    fclose(fp);
    uint32_t window_size = 256;
    uint32_t hop_size = 128;
    double* WVD = wvd(data, num_samples, window_size, hop_size);
    free(data);
    free(WVD);
    return 0;
}

其中，`wvd` 函数用于计算 Wigner-Ville 分布，采用了固定长度的窗口和固定步长的重叠方式。对于每个时间点，我们计算窗口内所有时域点之间的 WVD 值，得到一个大小为 `window_size * window_size` 的矩阵，并将该矩阵展平成一个大小为 `window_size * window_size` 的向量。最终，我们得到了一个大小为 `num_samples * window_size * window_size` 的 WVD 矩阵。注意，这个程序只计算了幅度平方，如果需要计算相位信息，需要将 `real` 和 `imag` 记录下来。