#define _USE_MATH_DEFINES

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

#define LENGTH 512 /* 信号長 ただし，2の整数乗 */

typedef struct {
    double r; /* 実部 */
    double i; /* 虚部 */
} Complex;  /* 複素数 */

Complex getComplex(double r, double i) {
    return (Complex){r, i};
}
Complex addComplex(Complex a, Complex b) {
    return getComplex(a.r + b.r, a.i + b.i);
}
Complex subComplex(Complex a, Complex b) {
    return getComplex(a.r - b.r, a.i - b.i);
}
Complex timeComplex(Complex a, Complex b) {
    return getComplex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}
Complex divComplexByReal(Complex z, double div) {
    return getComplex(z.r / div, z.i / div);
}
Complex cexptheta(double theta) {
    return getComplex(cos(theta), sin(theta));
}

void applyWindow(Complex *x, int N) {
    for (int i = 0; i < N; i++) {
        double w = 0.5 - 0.5 * cos(2.0 * M_PI * i / (N - 1));
        x[i].r *= w;
        x[i].i *= w;
    }
}

void bit_reversal(Complex *x, int N) {
    int i, j, k;
    Complex temp;
    j = 0;
    for (i = 0; i < N; i++) {
        if (i < j) {
            temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
        k = N >> 1;
        while (k <= j && k > 0) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

void fft(Complex *x, int N, int inverse) {
    bit_reversal(x, N);
    for (int stage = 1; (1 << stage) <= N; stage++) {
        int block = 1 << stage;
        int block2 = block >> 1;
        double theta = (inverse ? 2.0 : -2.0) * M_PI / block;
        Complex wm = cexptheta(theta);

        for (int k = 0; k < N; k += block) {
            Complex w = getComplex(1.0, 0);
            for (int j = 0; j < block2; j++) {
                Complex u = x[k + j];
                Complex v = timeComplex(w, x[k + j + block2]);

                x[k + j] = addComplex(u, v);
                x[k + j + block2] = subComplex(u, v);

                w = timeComplex(w, wm);
            }
        }
    }
    if (inverse) {
        for (int i = 0; i < N; i++)
            x[i] = divComplexByReal(x[i], N);
    }
}

void dft(Complex *x, Complex *c, int N) {
    for (int k = 0; k < N; k++) {
        c[k] = getComplex(0.0, 0.0);
        for (int n = 0; n < N; n++) {
            double theta = -2.0 * M_PI * k * n / N;
            Complex w = cexptheta(theta);
            c[k] = addComplex(c[k], timeComplex(x[n], w));
        }
    }
}

/* 振幅スペクトルの最大値を0dBとしてテキストファイルに出力する関数 */
void outputSpectrum(const char *fname, Complex *x, int N, double fs) {
    double max_amp = 0.0;
    double *amps = malloc((N / 2 + 1) * sizeof(double));

    // 全データの中で最大振幅を探す
    for (int i = 0; i <= N / 2; i++) {
        amps[i] = sqrt(x[i].r * x[i].r + x[i].i * x[i].i);
        if (amps[i] > max_amp)
            max_amp = amps[i];
    }

    FILE *fp = fopen(fname, "w");
    if (!fp) {
        free(amps);
        return;
    }
    for (int i = 0; i <= N / 2; i++) {
        double f = (double)i * fs / N;
        // 最大振幅を基準(X_b)としてdB値を計算
        double db = (amps[i] > 0 && max_amp > 0) ? 20.0 * log10(amps[i] / max_amp) : -100.0;
        fprintf(fp, "%f %f %f\n", f, amps[i], db);
    }
    fclose(fp);
    free(amps);
}

/* 複素数配列の平均パワーを求める関数 */
double calcAmpPower(Complex *x, int N) {
    double amp, power = 0.0;
    for (int i = 0; i < N; i++) {
        amp = x[i].r * x[i].r + x[i].i * x[i].i;
        power += amp;
    }
    return power / N;
}

int main(void) {
    int N = LENGTH;
    short input[LENGTH];
    Complex x_dft[LENGTH], x_fft[LENGTH], c[LENGTH];
    FILE *fp;

    double fs = 16000.0;

    printf("Processing Kadai 1...\n");
    fp = fopen("sample01.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open sample01.sw\n");
        return 1;
    }

    FILE *out1 = fopen("kadai1.txt", "w");
    short val;
    for (int i = 0; i < 480; i++) {
        if (fread(&val, sizeof(short), 1, fp) != 1)
           break;
        double time_ms = (double)i / fs * 1000.0;
        double norm_val = (double)val / 32768.0;
        fprintf(out1, "%f %f\n", time_ms, norm_val);
    }
    fclose(out1);
    fclose(fp);

    printf("Processing Kadai 2...\n");
    fp = fopen("sample01.sw", "rb");
    fread(input, sizeof(short), N, fp);
    fclose(fp);

    for (int i = 0; i < N; i++) {
        x_dft[i] = getComplex(input[i], 0);
        x_fft[i] = getComplex(input[i], 0);
    }

    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);

    QueryPerformanceCounter(&start);
    dft(x_dft, c, N);
    QueryPerformanceCounter(&end);
    printf("Kadai 2: DFT Elapsed Time: %f s\n", (double)(end.QuadPart - start.QuadPart) / freq.QuadPart);

    QueryPerformanceCounter(&start);
    fft(x_fft, N, 0);
    QueryPerformanceCounter(&end);
    printf("Kadai 2: FFT Elapsed Time: %f s\n", (double)(end.QuadPart - start.QuadPart) / freq.QuadPart);

    outputSpectrum("kadai2_dft.txt", c, N, fs);
    outputSpectrum("kadai2_fft.txt", x_fft, N, fs);

    printf("Processing Kadai 3...\n");
    Complex x_win[LENGTH];
    for (int i = 0; i < N; i++)
        x_win[i] = getComplex(input[i], 0);
    applyWindow(x_win, N);
    fft(x_win, N, 0);

    outputSpectrum("kadai3_window.txt", x_win, N, fs);

    printf("Processing Kadai 4...\n");
    typedef struct {
        double time_sec;
        const char *vowel;
    } SeekTarget;
    // 各母音の音声の開始時刻（秒）
    SeekTarget targets[] = {
        {0.361, "a"}, {0.969, "i"}, {1.415, "u"}, {1.936, "e"}, {2.431, "o"}
    };

    fp = fopen("j22330.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open j22330.sw\n");
        return 1;
    }

    for (int t = 0; t < 5; t++) {
        // 開始時刻からファイル上のバイト位置を計算 (16kHz, 16bit=2byte)
        long start_byte = (long)(targets[t].time_sec * fs * 2);
        fseek(fp, start_byte, SEEK_SET);
        if (fread(input, sizeof(short), N, fp) != (size_t)N)
            continue;

        Complex x_k4[LENGTH];
        for (int i = 0; i < N; i++)
            x_k4[i] = getComplex(input[i], 0);
        applyWindow(x_k4, N);
        fft(x_k4, N, 0);

        char fname[256];
        sprintf(fname, "kadai4_%s.txt", targets[t].vowel);
        outputSpectrum(fname, x_k4, N, fs);
    }
    fclose(fp);

    printf("Processing Kosatsu 1...\n");
    fp = fopen("sample01.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open sample01.sw\n");
        return 1;
    }
    fread(input, sizeof(short), N, fp);
    fclose(fp);

    Complex x_power[LENGTH];
    for (int i = 0; i < N; i++)
        x_power[i] = getComplex(input[i], 0);

    double beforePower, afterPower;
    beforePower = calcAmpPower(x_power, N);
    applyWindow(x_power, N);
    afterPower = calcAmpPower(x_power, N);

    printf("Kosatsu 1: BEFORE : %16f\n", beforePower);
    printf("Kosatsu 1: AFTER  : %16f\n", afterPower);
    printf("Kosatsu 1: DIFF   : %16f\n", beforePower - afterPower);
    printf("Kosatsu 1: RATIO  : %16f\n", afterPower / beforePower);

    printf("Processing Hatten Kadai 1...\n");
    fp = fopen("j22330.sw", "rb");
    if (!fp) {
        fprintf(stderr, "Cannot open j22330.sw\n");
        return 1;
    }
    fseek(fp, 0, SEEK_END);
    long file_size = ftell(fp);
    rewind(fp);
    long total_samples = file_size / sizeof(short);
    short *all_input = malloc(file_size);
    if (!all_input) {
        fprintf(stderr, "Memory allocation failed\n");
        fclose(fp);
        return 1;
    }
    fread(all_input, sizeof(short), total_samples, fp);
    fclose(fp);

    int shift = 160; // フレームをずらす間隔 (10ms)
    double global_max = 0.0;
    Complex x_hatten[LENGTH];

    // 0dBの基準となる全体の最大振幅を探索
    for (long i = 0; i <= total_samples - N; i += shift) {
        for (int j = 0; j < N; j++) {
            x_hatten[j] = getComplex(all_input[i + j], 0);
        }
        applyWindow(x_hatten, N);
        fft(x_hatten, N, 0);

        for (int j = 0; j <= N / 2; j++) {
            double amp = sqrt(x_hatten[j].r * x_hatten[j].r + x_hatten[j].i * x_hatten[j].i);
            if (amp > global_max) {
                global_max = amp;
            }
        }
    }

    // 各フレームのFFTを実行し、最大振幅を基準としたdB値を出力
    FILE *out_sp = fopen("spectrogram.txt", "w");
    if (out_sp && global_max > 0) {
        for (long i = 0; i <= total_samples - N; i += shift) {
            double time_sec = (i + N / 2.0) / fs;
            for (int j = 0; j < N; j++) {
                x_hatten[j] = getComplex(all_input[i + j], 0);
            }
            applyWindow(x_hatten, N);
            fft(x_hatten, N, 0);

            for (int j = 0; j <= N / 2; j++) {
                double f = (double)j * fs / N;
                // 振幅の計算とdB変換 (-90dBで下限クリッピング)
                double amp = sqrt(x_hatten[j].r * x_hatten[j].r + x_hatten[j].i * x_hatten[j].i);
                double db = (amp > 0) ? 20.0 * log10(amp / global_max) : -100.0;
                if (db < -90.0) db = -90.0;
                fprintf(out_sp, "%f %f %f\n", time_sec, f, db);
            }
            fprintf(out_sp, "\n"); // gnuplotのpm3d（面グラフ）描画用に空行を挿入
        }
        fclose(out_sp);
    }
    free(all_input);

    printf("All Kadai processed.\n");
    return 0;
}