Fourier transform a wav file.

dft_wav.c

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <fftw3.h>

typedef struct __attribute__((__packed__)) {
   char     ChunkID[4];
   uint32_t ChunkSize;
   char     Format[4];
   char     Subchunk1ID[4];
   uint32_t Subchunk1Size;
   uint16_t AudioFormat;
   uint16_t NumChannels;
   uint32_t SampleRate;
   uint32_t ByteRate;
   uint16_t BlockAlign;
   uint16_t BitsPerSample;

} wav_header_struct;

/* prints the structure of a wav_struct */
void print_wav_struct(wav_header_struct wav_file);

/* num of samples */
uint32_t NumOfSamples(wav_header_struct wav_file, FILE *fp);

/* subchunk2 size */
uint32_t Subchunk2Size(FILE *fp);

int main(){
   int16_t *raw_PCM;

   FILE *fp = fopen("write.wav", "rb");
   if (!fp) return 1;
   wav_header_struct wav_file1;

   fread(&wav_file1,sizeof(wav_file1),1,fp);
   if (strncmp(wav_file1.ChunkID, "RIFF", 4)) return 2; // not a riff file
   if (wav_file1.NumChannels == 2) return 3; // must be a mono file!
   // print_wav_struct(wav_file1);

   // printf("NumOfSamples %d\n", NumOfSamples(wav_file1, fp));
   uint32_t wav_samples = NumOfSamples(wav_file1, fp);
   fftw_complex *in, *out;
   in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*wav_samples);

   raw_PCM = malloc(sizeof(int16_t)*wav_samples);
   fread(raw_PCM, sizeof(int16_t), wav_samples, fp);
/*   for (int i = 0; i < 44100; i++){
      printf("%.3f   %d\n", (double) (i/44100.0), raw_PCM[i]);
   }
*/
   fclose(fp);
   for (int i = 0; i < wav_file1.SampleRate; i++){
      in[i][0] = (double) raw_PCM[i];
      in[i][1] = (double) raw_PCM[i];
   }
   free(raw_PCM);
   out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*wav_samples);

   fftw_plan my_plan = fftw_plan_dft_1d(wav_file1.SampleRate, in, out, FFTW_FORWARD, FFTW_PATIENT);

   fftw_execute(my_plan);
   fftw_execute(my_plan);
   for (int i = 0; i < 500; i++){
      printf("%d   %f\n", i, sqrt((out[i][0]*out[i][0])+(out[i][1]*out[i][1])));
   }
   fftw_destroy_plan(my_plan);
   fftw_free(in); fftw_free(out);
   return 0;
}

void print_wav_struct(wav_header_struct wav_file){
   printf("ChunkID:         %.*s\n", 4, wav_file.ChunkID);
   printf("ChunkSize:       %d\n", wav_file.ChunkSize);
   printf("Format:          %.*s\n", 4, wav_file.Format);
   printf("Subchunk1ID:     %.*s\n", 4, wav_file.Subchunk1ID);
   printf("Subchunk1Size:   %d\n", wav_file.Subchunk1Size);
   printf("AudioFormat:     %hu\n", wav_file.AudioFormat);
   printf("NumChannels:     %hu\n", wav_file.NumChannels);
   printf("SampleRate:      %d\n", wav_file.SampleRate);
   printf("ByteRate:        %d\n", wav_file.ByteRate);
   printf("BlockAlign:      %hu\n", wav_file.BlockAlign);
   printf("BitsPerSample:   %hu\n", wav_file.BitsPerSample);

}

uint32_t NumOfSamples(wav_header_struct wav_file, FILE *fp){

   uint32_t dummy_NumOfSamples =  (uint32_t) (8*((double) Subchunk2Size(fp)/wav_file.BitsPerSample));
   return dummy_NumOfSamples;
}

uint32_t Subchunk2Size(FILE *fp){
   char Subchunk2ID[4];
   uint32_t wav_size;

   fseek(fp, 36, SEEK_SET);

   while (strncmp(Subchunk2ID, "data", 4)){
      fread(Subchunk2ID, 4, 1, fp);
      fseek(fp, -3, SEEK_CUR);
   }
   fseek(fp, 3, SEEK_CUR);
   fread(&wav_size, 4, 1, fp);

   return wav_size;
}