title: C# .wav file data analysis & draw Freq.Domain graph from it
body: <p>my first question in stackoverflow...so pls inform me, if i violate the forum rules</p>
  
  <p>i try to analyze the audio .wav data (basically a sinusoid with/without Noise) in FFT &amp; draw the result as a Freq. Domain graph using C#. I need to analyze the sine wave &amp; find whether any noise is present or not for 10 to 20seconds audio data. </p>
  
  <p>Sine wave was generated using Audacity with the following properties.</p>
  
  <ul>
  <li>1KHz Sine Tone .wav data for 15 to 20s.</li>
  <li>16Bits per sample</li>
  <li>44100Hz</li>
  <li>Mono Channel</li>
  </ul>
  
  <p>Here is the code part to separate the Audio data from the Header data.</p>
  
  <pre><code>        // --- Open &amp; Read .wav File from the Path
          System.IO.FileStream WaveFile = System.IO.File.OpenRead(waveFilePath);
          System.IO.BinaryReader br = new System.IO.BinaryReader(WaveFile);
  
  
          // Convert the length of the file into Byte value.
          waveLength = WaveFile.Length;
          wave = new byte[waveLength];
  
          // --- Subtract the .wav File Header Info (First 44 Bytes) to get the real Audio Data.  Source -&gt; http://www.topherlee.com/software/pcm-tut-wavformat.html 
          wavenew = new byte[(waveLength - 44)];
          for (long i = 0; i &lt; waveLength; i++)
          {
              wave[i] = br.ReadByte();
              //Get only the Audio Data byte values
              if (i &gt;= 44)
              {
                  wavenew[i - 44] = wave[i];
              }
          }
  
          // Bits per Sample(whether 8 Bit or 16 Bit.)- Header Information35, 36 is for Bits per Sample. 
          bitsPerSample = BitConverter.ToInt16(wave, 34); //16Bit
  
          //Number of Channels - Header Information 23, 24 is for No of Channels.  (Mono or Stereo)
          noOfChannels = BitConverter.ToInt16(wave, 22); //Mono
  
  
          // Calculate the over all Audio Data length, here for FFT only values of exponential 2 can be given. So take 2 power 17 = 131072 for this calculation. Remaining Bytes will be neglected. 
          // Divide by 4 since every audio data contains 4 bytes. (Mono,16 Bit or Stereo,8 bit - 4 Bytes. Mono, 8 bit - 2 Bytes)
          dataLength = (131116-44)/4;  // /4;
          data = new double[dataLength]; 
  
          //SampleRate - Header Information 25 to 28 (4 Bytes) are for SampleRates
          SampleRate = BitConverter.ToInt32(wave, 24);     // 44100;
  
          // The every 4 Bytes contain the data of Sinus wave. 
          for (int i = 0; i &lt; dataLength; i++)
          {
              data[i] = BitConverter.ToInt32(wavenew, (1 + i) * 4);
          }
  </code></pre>
  
  <p>Here are example raw Audio data after conversion.</p>
  
  <p>890381306
  1403864333
  1804165263
  2058843210
  2147319456
  2063433407
  1812231337
  1414947118
  903955929
  320087260
  -289603352
  -875963166
  -1391936590
  -1795514361
  -2054321130
  -2147384959</p>
  
  <p>For FFT,  I use  <a href="http://www.lomont.org/Software/Misc/FFT/LomontFFT.html" rel="nofollow">http://www.lomont.org/Software/Misc/FFT/LomontFFT.html</a> as a reference, which returns a Double -Real &amp; Imag Value. </p>
  
  <pre><code>    public void RealFFT(double[] data, bool forward)                                                         
      {                                                                                                        
  
          var n = data.Length; // # of real inputs, 1/2 the complex length                                     
          // checks n is a power of 2 in 2's complement format                                                 
          if ((n &amp; (n - 1)) != 0)                                                                              
              throw new ArgumentException(                                                                     
                  "data length " + n + " in FFT is not a power of 2"                                           
                  );                                                                                           
  
          var sign = -1.0; // assume inverse FFT, this controls how algebra below works                        
          if (forward)                                                                                         
          { // do packed FFT. This can be changed to FFT to save memory                                        
              TableFFT(data, true);                                                                            
              sign = 1.0;                                                                                      
              // scaling - divide by scaling for N/2, then mult by scaling for N                               
              if (A != 1)                                                                                      
              {                                                                                                
                  var scale = Math.Pow(2.0, (A - 1) / 2.0);                                                    
                  for (var i = 0; i &lt; data.Length; ++i)                                                        
                      data[i] *= scale;                                                                        
              }                                                                                                
          }                                                                                                    
  
          var theta = B * sign * 2 * Math.PI / n;                                                              
          var wpr = Math.Cos(theta);                                                                           
          var wpi = Math.Sin(theta);                                                                           
          var wjr = wpr;                                                                                       
          var wji = wpi;                                                                                       
  
          for (var j = 1; j &lt;= n/4; ++j)                                                                       
          {                                                                                                    
              var k = n / 2 - j;                                                                               
              var tkr = data[2 * k];    // real and imaginary parts of t_k  = t_(n/2 - j)                      
              var tki = data[2 * k + 1];                                                                       
              var tjr = data[2 * j];    // real and imaginary parts of t_j                                     
              var tji = data[2 * j + 1];                                                                       
  
              var a = (tjr - tkr) * wji;                                                                       
              var b = (tji + tki) * wjr;                                                                       
              var c = (tjr - tkr) * wjr;                                                                       
              var d = (tji + tki) * wji;                                                                       
              var e = (tjr + tkr);                                                                             
              var f = (tji - tki);                                                                             
  
              // compute entry y[j]                                                                            
              data[2 * j] = 0.5 * (e + sign * (a + b));                                                        
              data[2 * j + 1] = 0.5 * (f + sign * (d - c));                                                    
  
              // compute entry y[k]                                                                            
              data[2 * k] = 0.5 * (e - sign * (b + a));                                                        
              data[2 * k + 1] = 0.5 * (sign * (d - c) - f);                                                    
  
              var temp = wjr;                                                                                  
              // todo - allow more accurate version here? make option?                                         
              wjr = wjr * wpr - wji * wpi;                                                                     
              wji = temp * wpi + wji * wpr;                                                                    
          }                                                                                                    
  
          if (forward)                                                                                         
          {                                                                                                    
              // compute final y0 and y_{N/2}, store in data[0], data[1]                                       
              var temp = data[0];                                                                              
              data[0] += data[1];                                                                              
              data[1] = temp - data[1];                                                                        
          }                                                                                                    
          else                                                                                                 
          {                                                                                                    
              var temp = data[0]; // unpack the y0 and y_{N/2}, then invert FFT                                
              data[0] = 0.5 * (temp + data[1]);                                                                
              data[1] = 0.5 * (temp - data[1]);                                                                
              // do packed inverse (table based) FFT. This can be changed to regular inverse FFT to save memory
              TableFFT(data, false);                                                                           
              // scaling - divide by scaling for N, then mult by scaling for N/2                               
              //if (A != -1) // todo - off by factor of 2? this works, but something seems weird               
              {                                                                                                
                  var scale = Math.Pow(2.0, -(A + 1) / 2.0)*2;                                                 
                  for (var i = 0; i &lt; data.Length; ++i)                                                        
                      data[i] *= scale;                                                                        
              }                                                                                                
          }                                                                                                    
      }                                                                                                                void Scale(double[] data, int n, bool forward)                                                           
      {                                                                                                        
          // forward scaling if needed                                                                         
          if ((forward) &amp;&amp; (A != 1))                                                                           
          {                                                                                                    
              var scale = Math.Pow(n, (A - 1) / 2.0);                                                          
              for (var i = 0; i &lt; data.Length; ++i)                                                            
                  data[i] *= scale;                                                                            
          }                                                                                                    
  
          // inverse scaling if needed                                                                         
          if ((!forward) &amp;&amp; (A != -1))                                                                         
          {                                                                                                    
              var scale = Math.Pow(n, -(A + 1) / 2.0);                                                         
              for (var i = 0; i &lt; data.Length; ++i)                                                            
                  data[i] *= scale;                                                                            
          }                                                                                                    
      }    
  
      public void TableFFT(double[] data, bool forward)                                                        
      {                                                                                                        
          var n = data.Length;                                                                                 
          // checks n is a power of 2 in 2's complement format                                                 
          if ((n &amp; (n - 1)) != 0)                                                                              
              throw new ArgumentException(                                                                     
                  "data length " + n + " in FFT is not a power of 2"                                           
                  );                                                                                           
          n /= 2;    // n is the number of samples                                                             
  
          Reverse(data, n); // bit index data reversal                                                         
  
          // make table if needed                                                                              
          if ((cosTable == null) || (cosTable.Length != n))                                                    
              Initialize(n);                                                                                   
  
          // do transform: so single point transforms, then doubles, etc.                                      
          double sign = forward ? B : -B;                                                                      
          var mmax = 1;                                                                                        
          var tptr = 0;                                                                                        
          while (n &gt; mmax)                                                                                     
          {                                                                                                    
              var istep = 2 * mmax;                                                                            
              for (var m = 0; m &lt; istep; m += 2)                                                               
              {                                                                                                
                  var wr = cosTable[tptr];                                                                     
                  var wi = sign * sinTable[tptr++];                                                            
                  for (var k = m; k &lt; 2 * n; k += 2 * istep)                                                   
                  {                                                                                            
                      var j = k + istep;                                                                       
                      var tempr = wr * data[j] - wi * data[j + 1];                                             
                      var tempi = wi * data[j] + wr * data[j + 1];                                             
                      data[j] = data[k] - tempr;                                                               
                      data[j + 1] = data[k + 1] - tempi;                                                       
                      data[k] = data[k] + tempr;                                                               
                      data[k + 1] = data[k + 1] + tempi;                                                       
                  }                                                                                            
              }                                                                                                
              mmax = istep;                                                                                    
          }                                                                                                    
  
  
          // perform data scaling as needed                                                                    
          Scale(data, n, forward);                                                                             
      }                                                                                                        
  </code></pre>
  
  <p>Amplitude is derived from the received Real &amp; Imag. value as below</p>
  
  <pre><code>     //Converting RealFFT return value in to an amplitude for Graph     
       (SqurRoot of Real*Real + Imag * Imag)
              realPart = new double[(dataLength / 2)];
              imagPart = new double[(dataLength / 2)];
              dataNew = new double[(dataLength / 2)];
              for (int i = 0; i &lt; dataNew.Length; i++)
              {
                  if (i == 0) //Ignore the first two values, since it contain the special values. Then calculate the magnitude for real &amp; imaginary part. 
                  {           
                  }
                  else
                  {
                      realPart[i] = data[(i * 2)] * data[(i * 2)];
                      imagPart[i] = data[((i * 2) + 1)] * data[((i * 2) + 1)];
                      dataNew[i] = Math.Sqrt(realPart[i] + imagPart[i]);
                  }
              }
              // Get these values for Y Axis as Amplitude
              Values.Add(dataNew);
              int N = dataNew.Length;
  
              double[] frequencies = new double[N];
              for (int i = 0; i &lt; N; i++ )
              {
                  if (i &lt; (N / 2))
                  {
                      frequencies[i] = i * SampleRate / N;
                  }
                  else if (i &gt;= (N / 2))
                  {
                      frequencies[i] = (N - i) * SampleRate / N;
                  }
  
                  //frequencies[i] = (i / N); //*10;
  
              }
  
              // Get the Frequency values for X Axis
              Values.Add(frequencies);
  </code></pre>
  
  <p>Processed Amplitude values are looking like..</p>
  
  <p>0
  9030724,08220743
  78204971,4566076
  11562334,8871099
  10855402,9273669
  9213306,99124749
  39816810,42806
  9154491,10446211
  10747893,2800474
  9744695,14696198
  140738122,900694</p>
  
  <p>Derived Frequency from the sample rate are as follows</p>
  
  <p>0
  2
  5
  8
  10
  13
  16
  18
  21
  24
  26
  29</p>
  
  <p>Finally the graph does not show the correct Freq range &amp; peaks are observed through out the  chart, even though the .wav file is a pure sine wave &amp; does not contain any noise.  (I´m not able to attach the graph image due to missing reputation, but i´m ready to post it.)</p>
  
  <p>As per Sinus Tone, there should be only one peak around 1KHz, but I´m seeing again sinus like signal.   I´ve checked with the different frequencies, but the result does not match with the expectation. </p>
  
  <p>I could not clearly figure out where I´m doing the mistake, whether</p>
  
  <ul>
  <li>Processing Audio Data (Splitting 4 Bytes after Header Info to Double)</li>
  <li>Processing data within FFT</li>
  <li>or Amplitude/Freq Graph calculation.</li>
  </ul>
  
  <p>I´m ready to provide additional information, if required &amp; any help would be appreciated.</p>