Seminar4

#include <cmath>
#include <cstdio>
#include <iostream>
#include <vector>

const double Pi = acos(-1);

struct Complex {
  double Re = 0.;
  double Im = 0.;

  Complex operator+(const Complex& right) const { return {Re + right.Re, Im + right.Im}; }
  Complex operator*(const Complex& right) const { return {Re * right.Re - Im * right.Im, Re * right.Im + Im * right.Re}; }
  Complex& operator+=(const Complex& right) { Re += right.Re; Im += right.Im; return *this; }
};

Complex operator*(double left, const Complex& right) { return {left * right.Re, left * right.Im}; }

std::vector<Complex> DFT(short* data, size_t n) {
  std::vector<Complex> wr(n);
  for (size_t i = 0; i < n; ++i) {
    wr[i].Re = cos(-2 * Pi * i / n);
    wr[i].Im = sin(-2 * Pi * i / n);
  }

  std::vector<Complex> dft(n);
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < n; ++j) {
      dft[j] += data[i] * wr[(i * j) % n];
    }
  }
  return dft;
}

std::vector<Complex> Zip(short* data, size_t size) {
  auto result = DFT(data, size);
  size_t cutted = size * 50 / 100;
  for (size_t i = cutted / 2; i < size - cutted / 2; ++i) result[i] = {0., 0.};
  return result;
}

std::vector<short> DFTReverse(const std::vector<Complex>& zipped) {
  size_t n = zipped.size();
  std::vector<Complex> w(n);
  for (size_t i = 0; i < n; ++i) {
    w[i].Re = cos(2 * Pi * i / n);
    w[i].Im = sin(2 * Pi * i / n);
  }

  std::vector<Complex> dft_rev(n);
  double coef = 1. / n;
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = 0; j < n; ++j) {
      dft_rev[j] += coef * zipped[i] * w[(i * j) % n];
    }
  }
  std::vector<short> restored(n);
  for (size_t i = 0; i < n; ++i) {
    restored[i] = static_cast<short>(dft_rev[i].Re);
  }
  return restored;

  return restored;
}

// Структура, описывающая заголовок WAV файла.
struct WAVHEADER {
  // WAV-формат начинается с RIFF-заголовка:

  // Содержит символы "RIFF" в ASCII кодировке
  // (0x52494646 в big-endian представлении)
  char chunkId[4];
  // 36 + subchunk2Size, или более точно:
  // 4 + (8 + subchunk1Size) + (8 + subchunk2Size)
  // Это оставшийся размер цепочки, начиная с этой позиции.
  // Иначе говоря, это размер файла - 8, то есть,
  // исключены поля chunkId и chunkSize.
  unsigned int chunkSize;
  // Содержит символы "WAVE"
  // (0x57415645 в big-endian представлении)
  char format[4];

  // Формат "WAVE" состоит из двух подцепочек: "fmt " и "data":
  // Подцепочка "fmt " описывает формат звуковых данных:

  // Содержит символы "fmt "
  // (0x666d7420 в big-endian представлении)
  char subchunk1Id[4];
  // 16 для формата PCM.
  // Это оставшийся размер подцепочки, начиная с этой позиции.
  unsigned int subchunk1Size;
  // Аудио формат, полный список можно получить здесь http://audiocoding.ru/wav_formats.txt
  // Для PCM = 1 (то есть, Линейное квантование).
  // Значения, отличающиеся от 1, обозначают некоторый формат сжатия.
  unsigned short audioFormat;
  // Количество каналов. Моно = 1, Стерео = 2 и т.д.
  unsigned short numChannels;
  // Частота дискретизации. 8000 Гц, 44100 Гц и т.д.
  unsigned int sampleRate;
  // sampleRate * numChannels * bitsPerSample/8
  unsigned int byteRate;
  // numChannels * bitsPerSample/8
  // Количество байт для одного сэмпла, включая все каналы.
  unsigned short blockAlign;
  // Так называемая "глубиная" или точность звучания. 8 бит, 16 бит и т.д.
  unsigned short bitsPerSample;
  // Подцепочка "data" содержит аудио-данные и их размер.
  // Содержит символы "data"
  // (0x64617461 в big-endian представлении)
  char subchunk2Id[4];
  // numSamples * numChannels * bitsPerSample/8
  // Количество байт в области данных.
  unsigned int subchunk2Size;

  // Далее следуют непосредственно Wav данные.
};

int main() {
  FILE *file = fopen("/Users/mstepa/tasks/FFT/speech.wav", "r");
  if (!file) {
    std::cout << "Failed open file";
    return 0;
  }

  WAVHEADER header;
  fread(&header, sizeof(WAVHEADER), 1, file);

  // Выводим полученные данные
  std::cout << header.chunkId[0] << header.chunkId[1] << header.chunkId[2] << header.chunkId[3] << std::endl;
  printf("Chunk size: %d\n", header.chunkSize);
  std::cout << header.format[0] << header.format[1] << header.format[2] << header.format[3] << std::endl;
  std::cout << header.subchunk1Id[0] << header.subchunk1Id[1] << header.subchunk1Id[2] << header.subchunk1Id[3] << std::endl;
  printf("Subchunk1Size: %d\n", header.subchunk1Size);
  printf("Audio format: %d\n", header.audioFormat);
  printf("Channels: %d\n", header.numChannels);
  printf("Sample rate: %d\n", header.sampleRate);
  printf("Bits per sample: %d\n", header.bitsPerSample);
  std::cout << header.subchunk2Id[0] << header.subchunk2Id[1] << header.subchunk2Id[2] << header.subchunk2Id[3] << std::endl;
  printf("Subchunk2Size: %d\n", header.subchunk2Size);

  // Посчитаем длительность воспроизведения в секундах
  double fDurationSeconds = 1. * header.subchunk2Size / (header.bitsPerSample / 8) / header.numChannels / header.sampleRate;
  int iDurationMinutes = (int)floor(fDurationSeconds) / 60;
  fDurationSeconds = fDurationSeconds - (iDurationMinutes * 60);
  printf("Duration: %02d:%02.3f\n", iDurationMinutes, fDurationSeconds);

  char* data = new char[header.subchunk2Size];
  fread(data, header.subchunk2Size, 1, file);

  std::cout << "Data is successfully loaded." << std::endl;

  std::cout << "Pi = " << Pi << std::endl;

  size_t n = header.subchunk2Size / 2;
  short* source = reinterpret_cast<short*>(data);
  auto zipped = Zip(source, n);
  std::cout << "Zipped first 100:" << std::endl;
  for (int i = 0; i < 100; ++i) std::cout << zipped[i].Re << " ";
  std::cout << "\nOther every 100:" << std::endl;
  for (int i = 100; i < zipped.size(); i += 100) std::cout << zipped[i].Re << " ";
  std::cout << std::endl;

  std::vector<short> restored = DFTReverse(zipped);

  for (int i = 0; i < 100; ++i) std::cout << source[i] << " ";
  std::cout << std::endl;
  for (int i = 0; i < 100; ++i) std::cout << restored[i] << " ";

  for (int i = 0; i < n; ++i) source[i] = restored[i];
  FILE* output_file = fopen("/Users/mstepa/tasks/FFT/speech_out.wav", "w");
  if (!output_file) {
    std::cout << "Failed to open output file" << std::endl;
    return 0;
  }
  fwrite(&header, sizeof(header), 1, output_file);
  fwrite(data, n * 2, 1, output_file);

  delete[] data;
  fclose(file);
  return 0;
}