#include "stdafx.h"
#include <stdio.h>
#include <math.h>

int main()
{
    double data_in1[128],data_in2[128],data_in3[128];
    double data_out1[128],data_out2[128],data_out3[128],data_out4[128],data_out5[128],data_out6[128];
    double data_out7[128],data_out8[128],data_out9[128],data_out10[128],data_out11[128],data_out12[128];
    double data_out13[128],data_out14[128],data_out15[128],data_out16[128],data_out17[128],data_out18[128];
    int k,N=128,n,an,bn;
    double sum1=0,sum2=0;
    double sita,pi;
    double data_powerspectrum[128],sum_real[128],sum_imag[128];
    int i=128;
    FILE *fp1;

    pi=3.14159265358979;


    if((fp1=fopen("data_sin_noise.txt","r"))==NULL){
        printf("Cannot open the file %s \n");
    }
    for(i=0;i<128;i++){
        fscanf(fp1,"%lf\n",&data_in1[i]);
        printf("[%d]=%f\n",i,data_in1[i]);
    }
    fclose(fp1);
    //data_inへsinの数値（データ）を入力


    for(k=0;k<N;k++){
        sum1=0;sum2=0;
        for(n=0;n<128;n++){
            sita=-1*(2*pi*k*n)/N;

            an=data_in1[n];
            bn=0;//フーリエの時はbnは0でいい。しかし逆フーリエの時は0ではない。jは無視していい。

            sum1+=an*cos(sita)-bn*sin(sita);
            sum2+=an*sin(sita)+bn*cos(sita);
        }
        data_out1[k]=sum1;
        data_out2[k]=sum2;
    }

    for(i=0;i<N;i++){
        data_out3[i]=sqrt(data_out1[i]*data_out1[i]+data_out2[i]*data_out2[i]);//正弦波パワースペクトル
    }

    for (k=0;k<N;k++){
        if(20<k&&k<128-20){
            data_out1[k]=0;           //フィルタ処理
            data_out2[k]=0;
        }
    }

    for(k=0;k<N;k++){
        data_out4[k]=sqrt(data_out1[k]*data_out1[k]+data_out2[k]*data_out2[k]);//処理後のパワースペクトル
    }

    for(i=0;i<N;i++){
        sum_real[i]=0;sum_imag[i]=0;
        for(n=0;n<N;n++){
            sita=(2*pi*k*n)/N;

            sum_real[i]+=(data_out1[n]*cos(sita))+(-1*data_out2[n]*sin(sita));
            sum_imag[i]+=(data_out2[n]*cos(sita))+(data_out1[n]*sin(sita));
        }
        data_out5[i]=sum_real[i]/N;//逆離散フーリエ変換
        data_out6[i]=sum_imag[i]/N;
    }

    if((fp1=fopen("data_out3.txt","w")) == NULL){  //正弦波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out3[i]);
        printf("[%d]=%f\n",i,data_out3[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out4.txt","w")) == NULL){  //フィルタ処理後の正弦波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out4[i]);
        printf("[%d]=%f\n",i,data_out4[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out5.txt","w")) == NULL){  //処理後に出力した正弦波信号の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out5[i]);
        printf("[%d]=%f\n",i,data_out3[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out6.txt","w")) == NULL){  //処理後に出力した正弦波出力の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out6[i]);
        printf("[%d]=%f\n",i,data_out6[i]);
    }
    fclose(fp1);



    if( ( fp1=fopen("data_square_noise.txt","r") ) == NULL ){
        printf ( "Cannot open the file %s \n" );
    }
    for(i=0; i<128; i++){
        fscanf(fp1,"%lf\n", &data_in2[i]);
        printf("[%d]=%f\n", i, data_in2[i]);
    }
    fclose( fp1 );
    //data_inへsquareの数値（データ）を入力



    for(k=0;k<N;k++){
        sum1=0;sum2=0;
        for(n=0;n<128;n++){
            sita=-1*(2*pi*k*n)/N;

            an=data_in2[n];
            bn=0;//フーリエの時はbnは0でいい。しかし逆フーリエの時は0ではない。jは無視していい。

            sum1+=an*cos(sita)-bn*sin(sita);
            sum2+=an*sin(sita)+bn*cos(sita);
        }
        data_out7[k]=sum1;
        data_out8[k]=sum2;
    }

    for(i=0;i<N;i++){
        data_out9[i]=sqrt(data_out7[i]*data_out7[i]+data_out8[i]*data_out8[i]);
    }
    for (k=0;k<N;k++){
        if(20<k&&k<128-20){
            data_out7[k]=0;
            data_out8[k]=0;
        }
    }
    for(k=0;k<N;k++){
        data_out10[k]=sqrt(data_out7[k]*data_out7[k]+data_out8[k]*data_out8[k]);
    }

    for(i=0;i<N;i++){
        sum_real[i]=0;sum_imag[i]=0;
        for(n=0;n<N;n++){
            sita=(2*pi*k*n)/N;

            sum_real[i]+=(data_out7[n]*cos(sita))+(-1*data_out7[n]*sin(sita));
            sum_imag[i]+=(data_out8[n]*cos(sita))+(data_out8[n]*sin(sita));
        }
        data_out11[i]=sum_real[i]/N;
        data_out12[i]=sum_imag[i]/N;
    }

    if((fp1=fopen("data_out9.txt","w")) == NULL){  //矩形波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out3[i]);
        printf("[%d]=%f\n",i,data_out3[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out10.txt","w")) == NULL){  //フィルタ処理後の矩形波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out3[i]);
        printf("[%d]=%f\n",i,data_out3[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out11.txt","w")) == NULL){  //処理後に出力した矩形波信号の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out9[i]);
        printf("[%d]=%f\n",i,data_out9[i]);
    }

    if((fp1=fopen("data_out12.txt","w")) == NULL){  //処理後に出力した矩形波信号の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out9[i]);
        printf("[%d]=%f\n",i,data_out9[i]);
    }

    fclose(fp1);



    if( ( fp1=fopen("data_triangle_noise.txt","r") ) == NULL ){
        printf ( "Cannot open the file %s \n" );
    }
    for(i=0; i<128; i++){
        fscanf(fp1,"%lf\n", &data_in3[i]);
        printf("[%d]=%f\n", i, data_in3[i]);
    }
    fclose( fp1 );
    //data_inへtriangleの数値（データ）を入力




    for(k=0;k<N;k++){
        sum1=0;sum2=0;
        for(n=0;n<128;n++){
            sita=-1*(2*pi*k*n)/N;

            an=data_in1[n];
            bn=0;//フーリエの時はbnは0でいい。しかし逆フーリエの時は0ではない。jは無視していい。

            sum1+=an*cos(sita)-bn*sin(sita);
            sum2+=an*sin(sita)+bn*cos(sita);
        }
        data_out13[k]=sum1;
        data_out14[k]=sum2;
    }


    for(i=0;i<N;i++){
        data_out15[i]=sqrt(data_out13[i]*data_out13[i]+data_out14[i]*data_out14[i]);
    }


    for (k=0;k<N;k++){
        if(20<k&&k<128-20){
            data_out1[k]=0;
            data_out2[i]=0;
        }
    }

    for(k=0;k<N;k++){
        data_out16[k]=sqrt(data_out13[k]*data_out13[k]+data_out14[k]*data_out14[k]);//√(real*real+img*img)
    }


    for(i=0;i<N;i++){
        sum_real[i]=0;sum_imag[i]=0;
        for(n=0;n<N;n++){
            sita=(2*pi*k*n)/N;

            sum_real[i]+=(data_out1[n]*cos(sita))+(-1*data_out2[n]*sin(sita));
            sum_imag[i]+=(data_out2[n]*cos(sita))+(data_out1[n]*sin(sita));
        }
        data_out17[i]=sum_real[i]/N;
        data_out18[i]=sum_imag[i]/N;
    }

    if((fp1=fopen("data_out15.txt","w")) == NULL){  //三角波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out15[i]);
        printf("[%d]=%f\n",i,data_out15[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out16.txt","w")) == NULL){  //フィルタ処理後の三角波のパワースペクトルの出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out16[i]);
        printf("[%d]=%f\n",i,data_out16[i]);
    }
    fclose(fp1);

    if((fp1=fopen("data_out17.txt","w")) == NULL){  //処理後に出力した三角波信号の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out17[i]);
        printf("[%d]=%f\n",i,data_out17[i]);
    }

    if((fp1=fopen("data_out18.txt","w")) == NULL){  //処理後に出力した三角波信号の出力
        printf("Cannot open the file %s \n");
    }
    for(i=0; i<128; i++){
        fprintf(fp1,"%f\n",data_out18[i]);
        printf("[%d]=%f\n",i,data_out18[i]);
    }

    fclose(fp1);





    return 0;
}