使用VS2019 fopen_s读取一个wav文件并对其进行短时傅里叶变换、小波变换和WVD变换c语言代码

以下是使用VS2019对wav文件进行短时傅里叶变换、小波变换和WVD变换的c语言代码示例：

短时傅里叶变换：

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

#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342
#define nfft 256
#define shift 128
#define winlen 256

typedef struct {
    short int format;
    short int channels;
    int sample_rate;
    int byte_rate;
    short int block_align;
    short int bits_per_sample;
    char data_chunk_header[4];
    int data_size;
} wav_header;

void stft(double *x, double *y, double *re, double *im) {
    int i, j, k;
    double w[nfft], xr[nfft], xi[nfft];

    for (i = 0; i < nfft; i++) {
        w[i] = 0.54 - 0.46 * cos(2 * PI * i / (nfft - 1));
    }

    for (i = 0; i < winlen; i++) {
        xr[i] = x[i] * w[i];
        xi[i] = y[i] * w[i];
    }

    for (i = winlen; i < nfft; i++) {
        xr[i] = 0;
        xi[i] = 0;
    }

    for (i = 0; i < nfft; i++) {
        re[i] = 0;
        im[i] = 0;
        for (j = 0; j < nfft; j++) {
            k = (i + j) % nfft;
            re[i] += xr[j] * cos(2 * PI * k / nfft) + xi[j] * sin(2 * PI * k / nfft);
            im[i] += -xr[j] * sin(2 * PI * k / nfft) + xi[j] * cos(2 * PI * k / nfft);
        }
    }
}

int main() {
    FILE *fp;
    wav_header header;
    char *filename = "test.wav";
    double x[winlen], y[winlen], re[nfft], im[nfft];
    short int buffer;
    int i, j, readsize;

    fp = fopen(filename, "rb");
    if (fp == NULL) {
        printf("Could not open file %s\n", filename);
        return -1;
    }

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

    while (strncmp(header.data_chunk_header, "data", 4) != 0) {
        fseek(fp, header.data_size, SEEK_CUR);
        fread(&header.data_chunk_header, sizeof(header.data_chunk_header), 1, fp);
    }

    printf("Channels: %d, Sample rate: %d, Bits per sample: %d\n", header.channels, header.sample_rate, header.bits_per_sample);

    readsize = winlen * header.channels * header.bits_per_sample / 8;

    while (fread(&buffer, sizeof(buffer), 1, fp) == 1) {
        for (i = 0; i < winlen; i++) {
            fread(&buffer, sizeof(buffer), 1, fp);
            if (header.channels == 2) {
                if (i % 2 == 0) {
                    x[i] = (double)buffer / 32768;
                } else {
                    y[i] = (double)buffer / 32768;
                }
            } else {
                x[i] = (double)buffer / 32768;
            }
        }

        if (i < winlen) {
            break;
        }

        stft(x, y, re, im);

        for (i = 0; i < nfft; i++) {
            printf("%f ", sqrt(re[i] * re[i] + im[i] * im[i]));
        }
        printf("\n");

        fseek(fp, -shift * header.channels * header.bits_per_sample / 8, SEEK_CUR);
    }

    fclose(fp);

    return 0;
}

小波变换：

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

#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342
#define nfft 256
#define shift 128
#define winlen 256

typedef struct {
    short int format;
    short int channels;
    int sample_rate;
    int byte_rate;
    short int block_align;
    short int bits_per_sample;
    char data_chunk_header[4];
    int data_size;
} wav_header;

void haar(double *x, double *y, double *output) {
    int i, j;

    for (i = 0; i < nfft; i++) {
        output[i] = 0;
        for (j = 0; j < nfft; j++) {
            output[i] += x[j] * pow(2, (i + 1) / 2.0 - 1) * pow(y[j], (i % 2 == 0 ? (i + 1) / 2.0 - 1 : i / 2.0));
        }
    }
}

int main() {
    FILE *fp;
    wav_header header;
    char *filename = "test.wav";
    double x[winlen], y[winlen], output[nfft];
    short int buffer;
    int i, j, readsize;

    fp = fopen(filename, "rb");
    if (fp == NULL) {
        printf("Could not open file %s\n", filename);
        return -1;
    }

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

    while (strncmp(header.data_chunk_header, "data", 4) != 0) {
        fseek(fp, header.data_size, SEEK_CUR);
        fread(&header.data_chunk_header, sizeof(header.data_chunk_header), 1, fp);
    }

    printf("Channels: %d, Sample rate: %d, Bits per sample: %d\n", header.channels, header.sample_rate, header.bits_per_sample);

    readsize = winlen * header.channels * header.bits_per_sample / 8;

    while (fread(&buffer, sizeof(buffer), 1, fp) == 1) {
        for (i = 0; i < winlen; i++) {
            fread(&buffer, sizeof(buffer), 1, fp);
            if (header.channels == 2) {
                if (i % 2 == 0) {
                    x[i] = (double)buffer / 32768;
                } else {
                    y[i] = (double)buffer / 32768;
                }
            } else {
                x[i] = (double)buffer / 32768;
            }
        }

        if (i < winlen) {
            break;
        }

        haar(x, y, output);

        for (i = 0; i < nfft; i++) {
            printf("%f ", output[i]);
        }
        printf("\n");

        fseek(fp, -shift * header.channels * header.bits_per_sample / 8, SEEK_CUR);
    }

    fclose(fp);

    return 0;
}

WVD变换：

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

#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342
#define nfft 256
#define shift 128
#define winlen 256

typedef struct {
    short int format;
    short int channels;
    int sample_rate;
    int byte_rate;
    short int block_align;
    short int bits_per_sample;
    char data_chunk_header[4];
    int data_size;
} wav_header;

void wvd(double *x, double *y, double *output) {
    int i, j;
    double xr[nfft], xi[nfft], yr[nfft], yi[nfft], zr[nfft], zi[nfft];

    for (i = 0; i < nfft; i++) {
        xr[i] = 0;
        xi[i] = 0;
        yr[i] = 0;
        yi[i] = 0;
        zr[i] = 0;
        zi[i] = 0;
    }

    for (i = 0; i < winlen; i++) {
        xr[i] = x[i];
        xi[i] = 0;
        yr[i] = y[i];
        yi[i] = 0;
    }

    for (i = 0; i < nfft; i++) {
        output[i] = 0;
        for (j = 0; j < nfft; j++) {
            zr[j] = xr[(i + j) % nfft] * yr[(nfft + i - j) % nfft];
            zi[j] = xi[(i + j) % nfft] * yi[(nfft + i - j) % nfft];
            output[i] += zr[j] - zi[j];
        }
    }
}

int main() {
    FILE *fp;
    wav_header header;
    char *filename = "test.wav";
    double x[winlen], y[winlen], output[nfft];
    short int buffer;
    int i, j, readsize;

    fp = fopen(filename, "rb");
    if (fp == NULL) {
        printf("Could not open file %s\n", filename);
        return -1;
    }

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

    while (strncmp(header.data_chunk_header, "data", 4) != 0) {
        fseek(fp, header.data_size, SEEK_CUR);
        fread(&header.data_chunk_header, sizeof(header.data_chunk_header), 1, fp);
    }

    printf("Channels: %d, Sample rate: %d, Bits per sample: %d\n", header.channels, header.sample_rate, header.bits_per_sample);

    readsize = winlen * header.channels * header.bits_per_sample / 8;

    while (fread(&buffer, sizeof(buffer), 1, fp) == 1) {
        for (i = 0; i < winlen; i++) {
            fread(&buffer, sizeof(buffer), 1, fp);
            if (header.channels == 2) {
                if (i % 2 == 0) {
                    x[i] = (double)buffer / 32768;
                } else {
                    y[i] = (double)buffer / 32768;
                }
            } else {
                x[i] = (double)buffer / 32768;
            }
        }

        if (i < winlen) {
            break;
        }

        wvd(x, y, output);

        for (i = 0; i < nfft; i++) {
            printf("%f ", output[i]);
        }
        printf("\n");

        fseek(fp, -shift * header.channels * header.bits_per_sample / 8, SEEK_CUR);
    }

    fclose(fp);

    return 0;
}