Tuesday, December 31, 2013

MD5 Revisited

Years ago I implemented MD5 in  C language as a college assignment. It was fun to implement something complex yet very straight forward. You just go over the RFC and  implement each step accurately. Only difficult part was debugging because even if you make a small mistake any where in the code the result is completely different and does not give any clue as to where the error is.

Years later I gave another look at the algorithm to get better understanding of the it. MD5 hash algorithm has three basic construction blocks.

1. Encryption Function: Converts a fixed sized input to fixed sized output such that output is a random permutation of the input. Encryption function also takes a key which can alter the permutations generated by the function. Encryption function by definition is reversible i.e. there exists a Decryption function such that given the output in can produce original input. In case of MD5, the encryption function takes a 32bit key and 128 bit input and 128 bit output. MD5 uses following encyptions function

def encrypt(f:(List[Int])=>Int,T:Int,S:Int)(data:List[Int],key:Int) =  
      data(3) :: 
      data(1) + rotateLeft((data(0)+f(data)+key+T),S) :: 
      data(1) :: 
      data(2) :: Nil
as you can see encryption function is parameterized by f,T and S. MD5 uses 64 different variations of encrypt function.

2. Compression Function: Compression Function takes a fixed sized input and converts it to a smaller fixed sized output. This function is one way i.e. there exists no function which can generate the input for a given output.Compression Functions are created using Encryption Function by applying it to a larger block again and again.

  In case of MD5 the compression function takes 512 bits as input and generates 128 bits of output.

 The 512 bit input is divided into 16 32 bit words. Each of these 32 bit words is used as key to the encryption function. This is repeated four times in four phases. In each phase a different permutation of 32 bit words is used. Thus, phases differ in two ways, first, the permutation of input and second, the encryption function used. Following function creates a four permutations, one for each phases and concatenates them. This results in 64 32 bit words.

  def permute(f:(Int)=>Int)(data:List[Int]):List[Int]={
      List.range(0,data.length).map((i)=>data(f(i)%data.length))
  }

  def permutations(chunk:List[Int]):List[Int]={
      permute((x:Int)=>x)(chunk):::
      permute((x:Int)=>(5*x+1))(chunk):::
      permute((x:Int)=>(3*x+5))(chunk):::
      permute((x:Int)=>(7*x))(chunk)
  }

The permute function takes a function and generates a permutation of data List according to it. permutations function returns four permutations one for each phase of compression function.

Once we have input data ready, we need to generate 64 different encryption functions, 16 for each phase. For this we need values for f,T and S. MD5 Compression function uses four different values of f, one for each phase.

  
  def F(x:Int,y:Int,z:Int) = x&y | ~(x)&z
  def G(x:Int,y:Int,z:Int) = x&z | y&(~z )
  def H(x:Int,y:Int,z:Int) = x ^ y ^ z
  def I(x:Int,y:Int,z:Int) = y ^ ( x | ~z )
  val fList = List(F _,G _,H _,I _).map(listyFy(_) _)

  def listyFy(f:(Int,Int,Int)=>Int)(d:List[Int]) = f(d(1),d(2),d(3))


  Above defines four different values of f. We also create a list these functions. listyFy function converts a function which takes three ints to a function which takes a list ints.

Following a list of functions which generate values for S for each phase.

val shiftfn=List( 
    (i:Int)=>5*i + 7,  (i:Int)=>(i*i + 7*i +10)/2,
    (i:Int)=>4 + 7*((i+1)/2) + 5*(i/2),  (i:Int)=>(i*i + 7*i +10)/2 + 1
  ).map((f)=>(x:Int)=>f(x%4));


Now below is the code which generates 64 encryption functions, which are used in each step of compression function .

def Ts(i:Int)=
      (math.floor(math.abs(math.sin(i+1))*(1L<<32)).toLong).toInt
val EList = 
      List.range(0,64).map((i:Int)=>encrypt(fList(i/16),
                                            Ts(i),
                       shiftfn(i/16)(i)) _)

Here we create EList, a list of 64 encryption function, by passing appropriate value of f,T and S to the encrypt function, which returns an encryption function corresponding to these values.
Now we can use this list of encryption function and permutation to define our compression Function.
  def compressMD5(state:List[Int],chunk:List[Byte]) = {
      val newstate =   permutations(asWords(chunk)).zip(EList).foldLeft(state)((x,y)=>y._2(x,y._1))
      newstate.zip(state).map({ case (a,b)=>(a+b)})
  }

compressMD5 takes a state, a list of four 32 bit words, and chunk, a list of 64 bytes. It converts 64 bytes in chunk to an array of 16, 32 bit words, using asWorkds utility function. It then creates permutations which results in 64 32 bit words, four permutation of 16 32 bit words. Then each of these words is applied as key to the corresponding encryption function in EList. The state variable is uses as initialization vector for this step which encrypted and then fed as data in to next encryption function. Finally, IV and resulting state and added word by word.


3. Hash Function: This is the final function which given input of any arbitrary length generates an output of fixed length. It usually uses merkel-damgard construction to apply a compression function to arbitrary large data.
   val initialState = List(0x67452301, 0xEFCDAB89, 0x98BADCFE, 0x10325476)
   def merkelDamgard(bytes:List[Byte]) = 
    (initialState /: packAs(64)( bytes++padding(bytes.length) ))(compressMD5)
This is the most straight forward step. It appends padding as required by merkelDamgard construction and divides resulting bytes in to chunks of 64 bytes which is equal to size of input of our compression function.

Sunday, November 17, 2013

Karplus–Strong string synthesis in JavaScript


Recently I came across this awesome algorithm which can generate sound of a plucked string, kind of like when picking guitar strings. Idea of generating sound of a guitar string fascinated me and I wanted to understand and finally implement the algorithm to see it working.

Before I go ahead and start explaining, you can try it out by click the button below.

Frequency: Hz.

There are plenty of places on the web that give details of  Karplus-Strong string synthesis algorithm. Basically, there are three main components that are required for this algorithm.
    1. Noise burst:  A source of white noise for a brief time. It is used to inject energy into the string, which happens when you pluck the string.
    2. Delay Line: A delay line which would re-inject the noise again into the circuit after certain time delay. The length of the time after which it should re-inject noise depends on the fundamental frequency of the string we want to simulate.
    3. Filter: A low pass filter to dampen out energy out of the string.

The algorithm by itself generates a stream of audio samples. These audio samples must be converted into an audio signal which will be sent to the speakers to produce sound. This can be done using Web Audio API which provides interfaces that take audio samples and send them to computer's hardware which can play them though the speakers. So, let us first get the setup ready which can send audio samples to the speaker.

The Setup

  if(window.webkitAudioContext){
       // Google Chrome
      var audioContext = new window.webkitAudioContext();
  } else {
       // Firefox
      var audioContext = new window.AudioContext();
  }

Here we create a global instance of AudioContext. It provides all the APIs require to interact with hardware.  There must be only one instance of AudioContext, which is why I have a created a global object for it.

if(audioContext.createJavaScriptNode){
         // Older API, Chrome
       var jsNode = audioContext.createJavaScriptNode(4096,1,1);
  } else {
        // Newer API, Firefox 
      var jsNode = audioContext.createScriptProcessor(4096,1,1);
  }


This creates a JavaScriptNode. JavaScriptNode can be used for sending an endless stream of audio samples to the hardware. AudioContext's createJavaScriptNode takes three parameters. First parameter is the size of the buffer which will be used to send data to the hardware. Think of it like the buffer you would normally use for reading files from file system. Value of this should be a power of two like 256, 512, 1024, 2048, 4096 etc. Second and third parameters are the number of input and output channels respectively. One thing to note here is that jsNode variable must also be global otherwise chrome's garbage collector will clean it up and jsNode will stop working once it is out of scope. This is probably a bug in chrome.


 
 jsNode.onaudioprocess = function(evt){
    var buffer = evt.outputBuffer.getChannelData(0);
    var n = buffer.length;
    for(var i = 0;i<n;i++){
         // write next sample to output
        buffer[i]=getNextSample();
    }
}



Next we define the callback that jsNode will call when it needs next buffersize number of samples. In this callback we iterate over the buffer array and fill it with new samples. The getNextSample() always returns the next sample in the sequence of samples that we want to send. All the magic now happens inside getNextSample() and we do not need to worry about buffersize, all that is taken care by this callback.

jsNode.connect(audioContext.destination);

And finally we connect our JavaScript node to the hardware. With this we are done with code that can send arbitrary audio samples to the hardware.

1. Noise Burst

// version 1.0
function getNextSample(){
      // return a random number between -1 and 1
    return 2*Math.random()-1;
} 


Noise source is nothing but a stream of samples with random values. Math.random() returns random number between 0 and 1. The above function returns a random number between -1 and 1. Now if you copy paste all of the above code in script tag of an html and load the page, you should hear some noise. In order to stop that noise you would have to close the tab since we are sending endless number of samples. For noise burst we should be able to send a fixed number of noise samples. Following version of getNextSample() does exactly that.


var noise_samples=0;
//version 2.0
function getNextSample(){
    var nextSample;
    if(noise_samples>0){
        nextSample= 2*Math.random() - 1;
        noise_samples--;
    } else {
        nextSample = 0;
    }
    return nextSample;
}

Now in some other function which can be called on click of a button we will set noise_samples to a finite value, so that after sometime it stops sending noise to audio device.


//version 1.0
function activate(){
   noise_samples = 10000;
}

With this we have implemented first component, a noise burst, of karplus strong algorithm. It can inject a finite number of samples of white noise. Try this in your browser to see if it works.

2. Delay Line

 Next we need to create a delay line. A delay line is used to re-introduce the samples back after a delay of time T. We can implement a delay line by storing the samples in an array. But how do we calculate number of samples the array should contain? In order to calculate number of samples we need to know the sampling rate of the audio card. Audio card plays a certain number of samples every second, called the Sample Rate. So to create a delay of time T we need T times he Sampling Rate of the Audio Card. sampleRate property in AudioContext tells us precisely that.

But wait, what is time T? The length of the delay line depends on the fundamental frequency of the string we are simulating, it is actually the length of the one complete cycle or time period of the wave with frequency same as the fundamental frequency. We know that time period of any wave is inverse of its frequency. So, let freq be the fundamental frequency of the wave so
      T = 1/freq
      Number of samples in delay line = T x SamplingRate
which means
       Number of samples in delay line = SamplingRate / freq
based on above equations we define following variables



var freq = 440; // 440Hz
var delayLineLen = Math.round(audioContext.sampleRate/freq); // number of sample has to be an integer, round out
var delayLine=[]; // array used creating delay line
var delayLinePos = 0; // marks the position at which the sample will be stored in delay line.


Now let us plug this delay line into our getNextSample() function.
var noise_samples=0;
//version 3.0
function getNextSample(){
    var nextSample;
    if(noise_samples>0){
        nextSample = 2*Math.random() - 1;
        noise_samples --;
    } else {
        nextSample = delayLine[delayLinePos];
    }
    delayLine[delayLinePos]=nextSample;
    delayLinePos = (delayLinePos + 1)%delayLineLen;
    return nextSample;
}

That's it. We have created a delay line and plugged that into out sample generator. Currently, delay line does nothing but stores the samples when noise samples are being sent and then plays them later when there are no noise samples. Since there is no damping or losses this creates an endless stream of same noise samples repeating at a time interval of 1/freq. You can try this out, you should  hear a clean note around 440Hz, and you may have to close the tab to shut it down.

3. Filter

 Now lets add our third component which is a filter or to be accurate, a low pass filter. A discrete low pass filter can be implemented using the formula. y[i] = y[i-1] + alpha*(x[i] - y[i-1]) where, alpha is the gain factor, usually ranges between 0 to 1,   y[i] is the output sample to be calculated y[i-1] is last calculated value x[i] is the current input value. Let us add this to our sample generator (getNextSample)

var noise_samples=0;
var alpha = 0.5;// you can try different values of alpha
//version 4.0
function getNextSample(){
    var nextSample;
    if(noise_samples>0){
        nextSample = 2*Math.random() - 1;
        noise_samples --;
    } else {
        var x_i = delayLine[delayLinePos];
        var y_i_1 = delayLine[(delayLinePos-1+delayLineLen)%delayLineLen]; // last sample in delayLine is yi minus 1 because it was created in last time.
        nextSample = y_i_1 + alpha*(x_i - y_i_1); // here nextSample is the y[i]
    }
    delayLine[delayLinePos]=nextSample; 
    delayLinePos = (delayLinePos + 1)%delayLineLen;
    return nextSample;
}

With this our implementation of Karplus-Strong algorithm is almost complete, there is one simple change we need to make to activate method.

// version 2.0
function activate(){
   noise_samples = delayLineLen; // fill entire delayLine and the stop.
}

If you now call activate method should hear the sound of a plucked string. You can play around with different frequencies and alpha values. Closer the alpha value to 1 longer the string would vibrate. Below is the complete source code listing.

/* *********************************************
 * 
 * Source Code Listing
 * Karplus-Strong String Synthesis
 * =============================================
 */
 // Global Variables
 var audioContext,
     jsNode,
     freq=440, // 440Hz
     delayLine=[], // array used creating delay line
     delayLineLen,
     noise_samples=0, // number of noise sample to play
     alpha = 0.5,
     delayLinePos = 0; // marks the position at which the sample will be stored in delay line.

 if(window.webkitAudioContext){
       // Google Chrome
     audioContext = new window.webkitAudioContext();
 } else {
       // Firefox
     audioContext = new window.AudioContext();
 }


 if(audioContext.createJavaScriptNode){
         // Older API, Chrome
     jsNode = audioContext.createJavaScriptNode(4096,1,1);
 } else {
        // Newer API, Firefox 
     jsNode = audioContext.createScriptProcessor(4096,1,1);
 }
 jsNode.onaudioprocess = function(evt){
     var buffer = evt.outputBuffer.getChannelData(0);
     var n = buffer.length;
     for(var i = 0;i<n;i++){
        // write next sample to output
         buffer[i]=getNextSample();
     }
 }

 jsNode.connect(audioContext.destination);

 delayLineLen = Math.round(audioContext.sampleRate/freq); // number of sample has to be an integer, round out

 function getNextSample(){
     var nextSample;
     if(noise_samples>0){
         nextSample = 2*Math.random() - 1;
         noise_samples--;
     } else {
         var x_i = delayLine[delayLinePos];
         
         // last sample in delayLine is yi minus 1 
         // because it was created in previous call to this function
         var y_i_1 = delayLine[(delayLinePos -1 + delayLineLen) % delayLineLen]; 
         nextSample = y_i_1 + alpha * (x_i - y_i_1); // here nextSample is the y[i]
     }
     delayLine[delayLinePos] = nextSample; 
     delayLinePos = (delayLinePos + 1) % delayLineLen;
     return nextSample;
 }
 // this function should be called on some
 // event, like onclick of a button
 function activate(){
     noise_samples = delayLineLen; // fill entire delayLine and the stop.
 }

/*************************************************************/