#include <stdio.h>
#include <math.h>
#include "kiss_fft.h"
#include "kiss_fftr.h"

void analyze_frequency(const float* signal, int N, float Fs) {
    // 分配FFT配置
    kiss_fft_cfg cfg = kiss_fft_alloc(N, 0, NULL, NULL);

    // 准备输入输出缓冲区
    kiss_fft_cpx* in = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * N);
    kiss_fft_cpx* out = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * N);

    // 填充输入（实数信号，虚部为0）
    for (int i = 0; i < N; i++) {
      float window = 0.5 * (1 - cos(2 * M_PI * i / (N - 1)));
      in[i].r = signal[i] * window;
      in[i].i = 0;
    }
  // 执行FFT
    kiss_fft(cfg, in, out);

    // 计算频率轴和幅度谱
    float df = Fs / N;  // 频率分辨率

    printf("频率分量分析（Fs = %.1f Hz, N = %d, Δf = %.3f Hz）:\n", Fs, N, df);
    printf("k\t频率(Hz)\t幅度\t相位(度)\n");

    // 通常只分析前N/2+1个点（正频率部分）
    for (int k = 0; k <= N/2; k++) {
        // 计算幅度
        float magnitude = sqrt(out[k].r * out[k].r + out[k].i * out[k].i);

        // 计算相位（弧度转角度）
        float phase = atan2(out[k].i, out[k].r) * 180.0 / M_PI;

        // 实际频率
        float frequency = k * df;

        // 归一化幅度（直流和奈奎斯特频率除外）
        float normalized_magnitude;
        if (k == 0 || k == N/2) {
            normalized_magnitude = magnitude / N;  // 直流和奈奎斯特频率
        } else {
            normalized_magnitude = magnitude / (N/2);  // 其他频率
        }

        // 只显示明显的频率分量（可选阈值）
        if (magnitude > 0.1 * (N/2)) {  // 阈值示例
            printf("%d\t%.2f\t\t%.4f\t%.1f\n", k, frequency, magnitude, phase);
        }
    }

    // 查找幅度峰值
    float threshold = 0.1 * (N/2);  // 阈值
    for (int k = 1; k < N/2; k++) {  // 跳过直流
        float magnitude = sqrt(out[k].r * out[k].r + out[k].i * out[k].i);
        if (magnitude > threshold) {
            printf("发现频率分量: %.2f Hz (幅度: %.4f)\n", k * df, magnitude);
        }
    }


    // 释放资源
    free(cfg);
    free(in);
    free(out);
}

void spectrum(float *signal, int n) {
    kiss_fftr_cfg cfg = kiss_fftr_alloc(n, 0, NULL, NULL);
    kiss_fft_cpx *spectrum = malloc((n/2 + 1) * sizeof(kiss_fft_cpx));

    kiss_fftr(cfg, signal, spectrum);

    printf("幅度谱:\n");
    for (int i = 0; i <= n/2; i++) {
        float mag = sqrt(spectrum[i].r * spectrum[i].r +
                        spectrum[i].i * spectrum[i].i);
        printf("%.4f\n", mag);
    }

    free(spectrum);
    free(cfg);
}

// int main() {
//     int n = 1024;
//     float Fs = 1000.0;    // 采样频率 1000 Hz
//     float *signal = malloc(n * sizeof(float));

//     // 50Hz + 120Hz 信号
//     for (int i = 0; i < n; i++) {
//         float t = (float)i / Fs;
//         signal[i] = 0.5 * cos(2 * M_PI * 50 * t) +
//                     0.3 * sin(2 * M_PI * 120 * t) +
//                     1.3 * sin(2 * M_PI * 300 * t) +
//                     0.9 * sin(2 * M_PI * 40 * t);
//     }

//     spectrum(signal, n);
//     analyze_frequency(signal, n, Fs);
//     free(signal);
//     return 0;
// }

int main() {
    // 示例参数
    float Fs = 1000.0;    // 采样频率 1000 Hz
    int N = 1024;         // FFT点数

    // 创建测试信号：50Hz + 120Hz正弦波
    float signal[N];
    for (int i = 0; i < N; i++) {
        float t = (float)i / Fs;
        signal[i] = 0.7 * sin(2 * M_PI * 50 * t) +  // 50Hz
                    0.3 * sin(2 * M_PI * 120 * t);  // 120Hz
    }

    // 分析频率
    analyze_frequency(signal, N, Fs);

    return 0;
}