#define _USE_MATH_DEFINES

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

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

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

/* プロトタイプ宣言 */
Complex getComplex(double r, double i);
void bit_reversal(Complex *x, int N);
void fft(Complex *x, int N, int inverse);
void genTwiddleTable(Complex*, int);
Complex calcTwiddlePower(Complex*, int, int);
void dft(Complex*, Complex*, int);
void applyWindow(Complex *x, int N);

Complex cexptheta(double theta);
Complex addComplex(Complex a, Complex b);
Complex subComplex(Complex a, Complex b);
Complex timeComplex(Complex a, Complex b);
Complex divComplexByReal(Complex z, double div);
Complex powComplex(Complex, int);

int main() {
    int N = LENGTH; /* 信号長 */
    double amp; /* 振幅 */

    int i;

    double start_sec = 2.454; /* 読み込む音声の開始位置 [秒] */

    long fs = 16000; /* 標本化周波数 */
    long bytes_per_sample = 2; /* 16bitなので2バイト */
    long start_byte = (long)(start_sec * fs * bytes_per_sample);    /* 指定時刻の開始位置バイト */

    /* 入力信号配列 */
    short input[LENGTH];
    Complex x[LENGTH];
    /* N点DFT係数 */
    Complex c[LENGTH];


    FILE *fp = fopen ( "j22330.sw", "rb" );
    if ( fp == NULL ) {
        fprintf ( stderr, "File open error (j22330.sw)\n" );
        exit ( 1 );
    }

    /* DFT係数の初期化 */
    for(i = 0; i < N; i++) {
        c[i] = getComplex(0.0, 0.0);
    }

    /* 読み込む位置にファイルポインタを移動 */
    fseek(fp, start_byte, SEEK_SET);

    /* 音声データを読み込み */
    if (fread(input, sizeof(short), N, fp) != (size_t)N) {
        fprintf(stderr, "File read error (j22330.sw)\n");
        fclose(fp);
        exit(1);
    }

    /* 音声データを複素数型に変換 */
    for(int i = 0; i < N; i++) {
        x[i] = getComplex(input[i], 0);
    }

    /* 窓関数の適用 */
    applyWindow(x, N);

    /* FFTの実行 */
    fft(x, N, 0);
    
    /* 振幅スペクトルの出力 */
    for(int i = 0; i < N; i++) {
        amp = sqrt(x[i].r * x[i].r + x[i].i * x[i].i);
        printf("%f\n", amp);
    }
    

    fclose(fp);

    return 0;
}

/**
 * @brief 配列の要素をビット逆順に並べ替える
 *
 * @param x 信号配列
 * @param N 信号長 ただし，2の整数乗
 */
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番目のビットが 1 だったら */
            j -= k; /* そのビットを 0 にする */
            k >>= 1;    /* 次のビットに移る */
        }
        j += k; /* 0 を見つけたら 1 にする */
    }

}

/**
 * @brief 高速フーリエ変換
 *
 * @param x 入力信号配列
 * @param N 信号長 ただし，2の整数乗
 * @param inverse 0ならDFT，1ならIDFT
 */
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];   /* X_0[k] */
                Complex v = timeComplex(w, x[k + j + block2]);  /* W^k * X_1[k] */

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

                w = timeComplex(w, wm);
            }
        }
    }

    /* 逆変換の場合、最後に N で割る */
    if (inverse) {
        for (int i = 0; i < N; i++) {
            x[i] = divComplexByReal(x[i], N);
        }
    }
}

/* 回転子W_Nの0乗からN-1乗を計算 */
void genTwiddleTable(Complex *W, int N) {
    int i;
    for (i = 0; i < N; i++) {
        W[i] = getComplex(cos(i * 2.0 * M_PI / N), -sin(i * 2.0 * M_PI / N));
    }
}

/* 回転子W_Nのn乗のk乗を計算 */
Complex calcTwiddlePower(Complex *W, int n, int k) {
    return powComplex(W[n], k);
}

/* 離散フーリエ変換 */
void dft(Complex *x, Complex *c, int N) {
    int k, n;
    Complex W[LENGTH];

    genTwiddleTable(W, N);

    /* 信号長で反復 */
    for(k = 0; k < N; k++) {
        /* シグマループ */
        for(n = 0; n < N; n++) {
            c[k] = addComplex(
                c[k],
                timeComplex(
                    x[n],
                    calcTwiddlePower(
                        W, k, n
                    )
                )
            );
        }
    }
}

/* 配列にハニング窓を適用 */
void applyWindow(Complex *x, int N) {
    for (int i = 0; i < N; i++) {
        double w = 0.5 - 0.5 * cos(2 * M_PI * i / (N - 1));
        x[i].r *= w;
        x[i].i *= w;
    }
}

Complex getComplex(double r, double i) {
    Complex c;
    c.r = r;
    c.i = i;
    return c;
}

Complex cexptheta(double theta) {
    return getComplex(cos(theta), sin(theta));
}

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);
}

/* 第一引数の第二引数乗を返す
 * 第二引数が負の場合は，第一引数の大きさが1のときのみ正しく計算される
 */
Complex powComplex(Complex z, int pow) {
    if(pow == 0) return getComplex(1.0, 0.0);

    int i;
    int p = (pow > 0) ? pow : -pow;
    Complex ret = z;
    for(i = 0; i < p-1; i++) {
        ret = timeComplex(ret, z);
    }

    /* |z|==1のとき，zの-1乗はzの共役複素数 */
    if(pow < 0)
        ret.i = -ret.i;

    return ret;
}