#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;
typedef unsigned int uint;

static void pc(c C a){printf("%10.0f + %10.0fi  ", __real__ a, __imag__ a);}
static void pcx(c C a){printf("%e + %ei\n", __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 C 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 = gi(ssndp+8, 4);
      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 = (ssndp + offset);

if(0) { //Print beginning samples in fixed point.
  short C * w = sd+2*numSampleFrames; // two shorts per sample frame
  printf("sd = %08x, w at %08x in file\n", (int)sd, (int)((void*)w - pp));
      {int j=100; while(j--) {printf("%5d %5d  ", *w, *(w+1));
              if(!(j&3)) printf("\n");
              w += 2;}}}
if(0) // Compute Fourier transform for 64K samples of left channel
 {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 = sz; while(k--) s[k] = *(sd+2*k);}
  pa("before", 100);
  fft(s, sz);
  pa("after", 100);
  {int j; for(j=0; j < sz; ++j) {double m = sq(__real__ s[j]) + sq(__imag__ s[j]);
    if(m > 3.e13) {pc(s[j]); printf("%0.6f %8.0f\n", (float)j/sz, m*1.e-8);}}}}
if (0) // Compute histogram of samples of right channel.
 {int i; short max=0, min=0; int mni = 0, mxi = 0;
    int qt[200]; {int j=200; while(j--) qt[j]=0;}
   for(i=0; i< numSampleFrames; ++i) {short q = *(sd+2*i+1);
     if(q<0) {if(q< min) {min = q; mni = i;}} else {if(q > max) {max = q; mxi = i;}}
     {int Q = ((int)q+100000)/200-500;
      if(abs(Q)>= 100) {printf("%08x \n", q); bitch("Too loud.");}
     ++qt[Q+100];}}
     printf("min = %d at %d, max = %d at %d\n", min, mni, max, mxi);
     {int j=200; while(j--) {printf("%d ", qt[j]); if(!(j&7)) printf("\n");}}}
if(0) { // Compute energy in each of many adjacent intervals.
   int is = 8820; // samples per interval
   short C * p = sd; int ic = 0;// sample cursor; interval counter
   double C f = (float)1/44100;
   while((void*)p+2*is<pend) {double s = 0; int k = is; while(k--) {
   s += sq(*p);
   p += 2;}
       printf("%4.1f %7.1f\n", is*f*(ic++), s*1.e-8);}}
if(1) { // Convolution of left and right channels
   int C n = 16; int C sz = 1<<n; c l[sz], r[sz];
   void ps(int k, int n){int j;
     for(j = k; j<n; ++j) {printf("%4d ", j); pcx(l[(j+sz)&(sz-1)]);}}
   {int k = sz; while(k--) {l[k] = *(sd+2*k); r[k] = *(sd+2*k+1);}}
  fft(l, sz); fft(r, sz);
  {int j = sz; while(j--) l[j] = conj(l[j]*r[j])/sz;}
  fft(l, sz);
  ps(-25, 25); ps(3000, 3030); ps(8000, 8030);
  {int j=sz; float q = 0; while(j--) q = fmax(q, fabs(__imag__ l[j]));
    printf("max imaginary = %e\n", q);}}

}}}}}

return 0;}