短时傅里叶变换谱图 C语言

短时傅里叶变换（Short-Time Fourier Transform, STFT）是一种时频分析方法，它将信号分成若干个短时窗口，对每个窗口内的信号进行傅里叶变换，然后将变换结果拼接在一起形成一个二维谱图，其中横轴表示时间，纵轴表示频率，颜色表示信号的强度。

以下是一个使用C语言实现的短时傅里叶变换谱图的示例代码：

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

#define PI 3.14159265358979323846

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

// 计算复数加法
complex_t add(complex_t a, complex_t b) {
    complex_t c;
    c.real = a.real + b.real;
    c.imag = a.imag + b.imag;
    return c;
}

// 计算复数乘法
complex_t mul(complex_t a, complex_t b) {
    complex_t c;
    c.real = a.real * b.real - a.imag * b.imag;
    c.imag = a.real * b.imag + a.imag * b.real;
    return c;
}

// 计算离散傅里叶变换
void dft(complex_t *x, int N) {
    complex_t *X = (complex_t *)malloc(N * sizeof(complex_t));
    for (int k = 0; k < N; k++) {
        X[k].real = 0;
        X[k].imag = 0;
        for (int n = 0; n < N; n++) {
            complex_t W;
            W.real = cos(2 * PI * k * n / N);
            W.imag = -sin(2 * PI * k * n / N);
            X[k] = add(X[k], mul(x[n], W));
        }
    }
    for (int k = 0; k < N; k++) {
        x[k] = X[k];
    }
    free(X);
}

// 计算短时傅里叶变换
void stft(double *x, int N, int M, int L, complex_t **S) {
    complex_t *window = (complex_t *)malloc(L * sizeof(complex_t));
    for (int n = 0; n < L; n++) {
        window[n].real = cos(2 * PI * n / L);
        window[n].imag = -sin(2 * PI * n / L);
    }
    for (int m = 0; m < M; m++) {
        complex_t *X = (complex_t *)malloc(L * sizeof(complex_t));
        for (int n = 0; n < L; n++) {
            X[n].real = x[m * L + n] * window[n].real;
            X[n].imag = x[m * L + n] * window[n].imag;
        }
        dft(X, L);
        for (int k = 0; k < L; k++) {
            S[m][k] = X[k];
        }
        free(X);
    }
    free(window);
}

// 绘制二维谱图
void plot(complex_t **S, int M, int L) {
    for (int k = 0; k < L / 2 + 1; k++) {
        for (int m = 0; m < M; m++) {
            double P = sqrt(S[m][k].real * S[m][k].real + S[m][k].imag * S[m][k].imag);
            printf("%lf ", P);
        }
        printf("\n");
    }
}

int main() {
    int N = 1024; // 信号长度
    int M = 32; // 窗口数量
    int L = 32; // 窗口长度

    double *x = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n++) {
        x[n] = sin(2 * PI * 50 * n / N) + sin(2 * PI * 100 * n / N);
    }

    complex_t **S = (complex_t **)malloc((L / 2 + 1) * sizeof(complex_t *));
    for (int k = 0; k < L / 2 + 1; k++) {
        S[k] = (complex_t *)malloc(M * sizeof(complex_t));
    }

    stft(x, N / L, M, L, S);
    plot(S, M, L);

    for (int k = 0; k < L / 2 + 1; k++) {
        free(S[k]);
    }
    free(S);
    free(x);

    return 0;
}

该示例代码生成一个由两个正弦波叠加而成的信号，并使用短时傅里叶变换将其分成32个窗口，每个窗口长度为32。最终输出一个32行、16列的二维谱图，其中每个元素表示对应时间和频率上的信号强度。