Discrete Fourier Transformation (DFT) based Spectrum Tracer

DFT can be used to convert audio samples from time domain to frequency domain. That is the dancing spectrum in graphic equalizers. This piece of code (C++ with MFC) traces the spectrum of a wave file.

const double PI = 3.14159265;
const int N = 1024; // DFT point count

void TraceSpectrum(CDC* pDC)
{
    char* pzFile ="D:\\in.wav"; 
    HMMIO h = mmioOpen((LPTSTR)pzFile, NULL, MMIO_READ | MMIO_ALLOCBUF | MMIO_DENYWRITE);
    for (int i = 0; i < 50000; ++i)
    {
        double realInput[N] = { 0 };
        for (int j = 0; j < N; ++j)
        {
            SHORT iLeft;
            SHORT iRight;
            mmioRead(h, (HPSTR)&iLeft, 2); // 16 bit samples
            mmioRead(h, (HPSTR)&iRight, 2); // 16 bit samples
            realInput[j] = iLeft / 1000.0;
        }

        TraceDFT(realInput, pDC);
        Sleep(10);
    }
    mmioClose(h, 0);
}

void TraceDFT(double* realInput, CDC* pDC) 

{
   double realOutput[N] = { 0 };
   double imagOutput[N] = { 0 }; 

   for (int k = 0; k < N / 2; ++k) // Output is symmetric, hence N/2
   {
      for (int n = 0; n < N; ++n)
      {
          // Theory: F(k) = sigma(T(n) * (cos(x) - j*sin(x))) Where x = 2PIkn/N double x = 2 * PI * k * n / N;
          double re = realInput[n] * cos(x);
          double im = realInput[n] * -sin(x);
          realOutput[k] += re;
          imagOutput[k] += im;
      }
   }

   // To Decibels
   double spectrum[N / 2];
   for (int i = 0; i < N / 2; ++i)
   {
      double d = sqrt(realOutput[i] * realOutput[i] + imagOutput[i] * imagOutput[i]);
      d = d < 1 ? 0 : 20 * log(d);
      spectrum[i] = min(d, 200);
   }

   // Draw the graph
   pDC->FillSolidRect(0, 0, 1000, 500, 0xffffff);
   for (int i = 0; i < N / 2; ++i)
      pDC->FillSolidRect(20 + i, 200 - spectrum[i], 1, spectrum[i], 0x0);
}