DFT

#include <stdio.h>
#include <stdlib.h>
#define _USE_MATH_DEFINES
#include <math.h>
static const int N = 128;
static FILE *open_for_read(const char *file_name) {
  FILE *in = fopen(file_name, "r");
  if (!in) {
    fprintf(stderr, "open err: %s\n", file_name);
    exit(1);
  }
  return in;
}
static void read_data(FILE *in, double *data) {
  int i;
  for (i = 0; i < N; i++)
    fscanf(in, "%lf", &data[i]);
}
static void dft(double *data, double *real, double *image) {
  int k, n;
  double theta;
  for (k = 0; k < N; k++) {
    real[k] = 0.;
    image[k] = 0.;
    for (n = 0; n < N; n++) {
      theta = 2 * M_PI * k * n / N;
      real[k] += data[n] * cos(theta);
      image[k] += - data[n] * sin(theta);
    }
  }
}
static void power_spectrum(double *real, double *image, double *result) {
  int i;
  for (i = 0; i < N; i++)
    result[i] = sqrt(real[i] * real[i] + image[i] * image[i]);
}
static void filter(double *real, double *image) {
  int i;
  for (i = 0; i < N; i++)
    if (i < 20 || N - 20 < i) {
      real[i] = 0.;
      image[i] = 0.;
    }
}
static void idft(double *real, double *image, double *idft_real, double *idft_image) {
  int k, n;
  double re, im, theta;
  for (k = 0; k < N; k++) {
    re = 0.;
    im = 0.;
    for (n = 0; n < N; n++) {
      theta = 2 * M_PI * k * n / N;
      re += real[n] * cos(theta) - image[n] * sin(theta);
      im += real[n] * sin(theta) + image[n] * cos(theta);
    }
    idft_real[k] = re / N;
    idft_image[k] = im / N;
  }
}
FILE *open_for_write(const char *file_name) {
  FILE *out = fopen(file_name, "w");
  if (!out) {
    fprintf(stderr, "open error: %s\n", file_name);
    exit(1);
  }
  return out;
}
static void print_data(FILE *out, double *data) {
  int i;
  for (i = 0; i < N; i++)
    fprintf(out, "%lf\n", data[i]);
}
static void process(const char *noize, const char *power, const char *power_filter,
    const char *idft_amp, const char *idft_wav) {
  FILE *in, *out;
  double data[N], real[N], image[N], power_spec[N], power_spec_filter[N], idft_real[N], idft_image[N];
  in = open_for_read(noize);
  read_data(in, data);
  fclose(in);
  dft(data, real, image);
  power_spectrum(real, image, power_spec);
  filter(real, image);
  power_spectrum(real, image, power_spec_filter);
  idft(real, image, idft_real, idft_image);
  out = open_for_write(power);
  print_data(out, power_spec);
  fclose(out);
  out = open_for_write(power_filter);
  print_data(out, power_spec_filter);
  fclose(out);
  out = open_for_write(idft_amp);
  print_data(out, idft_real);
  fclose(out);
  out = open_for_write(idft_wav);
  print_data(out, idft_image);
  fclose(out);
}
int main(void) {
  const char *sine_noize = "data_sin_noise.txt";
  const char *sine_power1 = "data_out3.txt";
  const char *sine_power2 = "data_out4.txt";
  const char *sine_idft_amp = "data_out5.txt";
  const char *sine_idft_wav = "data_out6.txt";
  const char *square_noize = "data_square_noise.txt";
  const char *square_power1 = "data_out9.txt";
  const char *square_power2 = "data_out10.txt";
  const char *square_idft_amp = "data_out11.txt";
  const char *square_idft_wav = "data_out12.txt";
  const char *triangle_noise = "data_triangle_noise.txt";
  const char *triangle_power1 = "data_out15.txt";
  const char *triangle_power2 = "data_out16.txt";
  const char *triangle_idft_amp = "data_out17.txt";
  const char *triangle_idft_wav = "data_out18.txt";
  process(sine_noize, sine_power1, sine_power2, sine_idft_amp, sine_idft_wav);
  process(square_noize, square_power1, square_power2, square_idft_amp, square_idft_wav);
  process(triangle_noise, triangle_power1, triangle_power2, triangle_idft_amp, triangle_idft_wav);
  return 0;
}