#include <stdio.h>
#include <sys/stat.h> // man 2 stat
#include <sys/types.h>
#include <sys/mman.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#define C const
typedef double _Complex c;
typedef long double _Complex cx;
typedef unsigned char uchar;

static void pc(c a){printf("%12.2f + %12.2fi  ", __real__ a, __imag__ a);}
static double sq(double x){return x*x;}

static void fft(c a[], C int n){
  {int j=0, i; for(i=0; i<n-1; ++i) {
    if(i<j) {c t=a[j]; a[j]=a[i]; a[i]=t;}
    {int m=n/2; while(m<=j) {j -= m; m /= 2;} j += m;}}} 
  {int j=1, m; cx b=-1; while(j<n) {cx w = 1;
    for(m=0; m<j; ++m) {int i;
      for(i=m; i<n; i+=j+j)
        {c t = w*a[i+j]; a[i+j]=a[i]-t; a[i] += t;}
      w *= b;}
    j *= 2; b = csqrtl(b);}}}

typedef struct{char n[4];} nm;
void pID(void C * n, char C * C s){int j;
   for(j=0; j<4; ++j) putchar(*((char *)(n++)));
   printf(" %s\n", s);}

// routines gi and cn below claim to be ambidextrous
// but I have no big endian machine to test them.
int gi(void C * n, int l){int q=0; int j; int w=(l<<3)-8;
  for(j=0; j<l; ++j) q |= (*(uchar *)(n+j))<<(w-j*8);
  return q;}

int cn(char C * C c, void C * C n){int j=4; int C w = *(int*)n;
   while(j--) if(*(c+j) != ((w>>(j<<3))&255)) return 0; return 1;}

static void ex(int C d){exit(-3);}
static void bitch(char C * C s){ex(fprintf(stderr, "%s: Bye.\n", s));}

int main(){
    if(!stdin) bitch("No input file!");
    {struct stat Ss; if(fstat(fileno(stdin), &Ss)) bitch("fstat failed");
    if (!Ss.st_size) bitch("Empty file.");
    {void C * C pp = mmap(0, Ss.st_size, PROT_READ, MAP_SHARED,
       fileno(stdin), 0);
     if(0) printf("file starts: %08x %08x\n", *(int *)pp, *(int *)(pp+4));
     if(!cn("FORM", pp)) bitch("Doesn't begin 'FORM'.");
     printf("Size from stat: 0x%x, size from FormAIFCChunk: 0x%08x\n",
         (int)Ss.st_size, gi(pp+4, 4));
     if(!cn("AIFC", pp+8)) bitch("Expect AIFC chunk; not there.");
     printf("Size from AFIC chunk: 0x%x\n",
         gi(pp+12, 4));
         
         
     {void C * q = pp+12; void C * C pend = pp + gi(pp+4, 4);
      void C * fverp=0, * commp=0, * ssndp=0;
        while(q<pend) {int len = gi(q+4, 4);
           pID(q, "nxt:");
           printf("%08x\n", len);
           if(cn("FVER", q)) fverp = q;
           else if(cn("COMM", q)) commp = q;
           else if(cn("SSND", q)) ssndp = q;
           else bitch("What is this?");
           q += ((len+1)&~1)+8;}
        printf("end q=%08x pend=%08x\n", (int)q, (int)pend);
     if(!(fverp && commp && ssndp)) bitch("Not a full house");
     if(gi(fverp+8, 4) != 0xA2805140) bitch("wrong version");
     {int numChannels = gi(commp+8, 2), numSampleFrames = gi(commp+10, 4);
      int sampleSize = gi(commp+14, 2),
         offset = 0?gi(ssndp+8, 4):0xa3e0;
      int sampleRateXp = gi(commp+16, 2), sampleRateFr = gi(commp+18, 4);
      printf("numChannels=%d, sampleSize=%d, numSampleFrames=%d\n",
               numChannels, sampleSize, numSampleFrames);
      pID(commp+26, "Compression Type");
      if(!cn("sowt", commp+26)) bitch("Write some new code.");
      {int q; for(q=0; q< *(char*)(commp+30); ++q) putchar(*(char*)commp+31+q);}
      if(0) printf("rate exponent = %x, high frac = %08x\n", sampleRateXp, sampleRateFr);
      printf("Blocksize = 0x%x offset = 0x%x\n", gi(ssndp+12, 4), offset);
      if(sampleRateXp != 0x400e || sampleRateFr != 0xac440000)
           bitch("Sample rate not 44100 Hz");
      if(sampleSize != 16) bitch("Wrong sampe size for this program");
      if(numChannels != 2) bitch("This program only for 2 channel sound");
      {short C * C sd = 0?(ssndp + offset):(pp+0xaad0);
      {short C * w = sd;
       printf("sd = %08x, w at %08x in file\n", (int)sd, (int)w - (int)pp);
        int j=400; while(j--) {printf("%5d %5d  ", *w, *(w+1));
              if(!(j&3)) printf("\n");
              w += 2;}}
{int C n = 16; int C sz = 1<<n; c s[sz];
  void pa(char * m, int z){int k=z; printf("a: %s\n", m);
  while(k--) {pc(s[z-k]); if(!(k&3)) printf("\n");}}
  {int k = 1<<n; while(k--) s[k] = *(sd+2*k);}
  pa("before", 100);
  fft(s, 1<<n);
  pa("after", 100);
  {int j; for(j=0; j < (1<<n); ++j) {double m = sq(__real__ s[j]) + sq(__imag__ s[j]);
    if(m > 1.e13) {pc(s[j]); printf("%0.6f %8.0f\n", (float)j/(1<<n), m*1.e-8);}}}}

}}}}}

return 0;}