Signal processing with Haskell : First steps

As said in my previous post, I have now started working on an Haskell package for prototyping signal processing algorithms : mainly speech and acoustic. None of the packages I have seen so far were matching my needs. So, I decided to start working on something. And, this project is also a good opportunity for me to learn the latest ghc stuff like data type promotion and list fusion.

This work is available on github but warning : I am in an exploratory phase. Anything can change at any time. Also, don't trust the code unless I have written some unit tests. In this first version I only have the foundations and it is only partially tested.

So, let's see what's in the package so far.

1. Signal

The core of the package is the Signal type. It is just a wrapper around lists (using the stream-fusion package). But, I want Signal to be infinite lists. An finite lists will be vectors in this project. So that's why I am not using the list type directly.

To avoid conflict with the Prelude, the standard list manipulation functions are followed with a S when used on a signal.

So, for instance, you can write:

 mapS (+2) mySignal 

to add a DC offset to the signal.

2. Fixed point arithmetic

Implementing a signal processing algorithm on a true DSP is tricky. That's why professional working on speech / acoustic are using the ITU G.191 specification : a set of C fixed point reference implementations with definition of basic operators (multiply, addition ...) with the fixed point behavior (precision, saturation).

At some point when studying an algorithm, and before porting it to a DSP, you need to study the behavior when passing to a fixed point representation (unless you're rich and using floating point arithmetic on your DSP :-). It is not very fun to do it using G.191 which is very low level and close to the DSP instruction sets.

So, one goal of this project is to support enough fixed point arithmetic to be able to easily study the behavior of an algorithm and keep a high level description.

For this, I am using the type:

Fixed a n s 

with:

a is Int16 or Int32 or Int40 or Int64

n is using type level numeral and so you need ghc 7.6 with data type promotion. So n has kind Nat

 n :: Nat 

s has the kind Saturation and can be the type Saturated or Unsaturated.

So, an example of fixed point type could be:

Fixed Int16 8 Saturated

It is the type of a signed number on 16 bits with 8 bits for the fractional part and which is using saturated arithmetic. So for instance

 0x7FFE + 2 -> 0x7FFF and not 0x8000

Saturation is introducing distortion but it is not as bad as a sign change. So, the tradeoff is between quantization noise if you don't have enough bits for the fractional part or saturation. Haskell will use the type level natural to automatically introduce the right shift in the computations needed to position the fixed point at the right place.

Here is an example generated with this library:

On the above picture you can see a signal computed with Double numbers and with different fixed point formats.

One is showing the quantization noise, another the saturation and another the sign change when saturation is disabled and there is not enough available dynamics.

All the curves above are generated with the same code :

genericSignal :: forall a. Sample a => Signal a
genericSignal = mapS (\t ->  2*(fromDouble $ 1.5*sin (2*pi*getT t) 
                       + 0.4*sin(2*pi*20*getT t))) theTimes

fromDouble is a simple way to convert a double into a fixed point. That's the reason for HasDoubleRepresentation. In a future version I may have a fixed point algorithm to compute the sine.

Sample is a constraint kind.

I have typed the time and frequency hence the getT.

Then, I can just write:

mySignalA :: Signal Double 
mySignalA = genericSignal 

mySignalC :: Signal (Fixed Int16 14 Saturated) 
mySignalC = genericSignal

mySignalD :: Signal (Fixed Int16 3 Saturated) 
mySignalD = genericSignal

mySignalE:: Signal (Fixed Int16 14 Unsaturated) 
mySignalE = genericSignal

for different versions.

The display (using my OS X viewer) is generated with:

display $ discreteSignalsWithStyle (takeWhileS (<= duration) theTimes) [ AS mySignalA
                                                                       , AS mySignalC
                                                                       , AS mySignalD 
                                                                       ]

theTimes being a Signal, it is an infinite sequence of time instants and we need to select a finite sequence for plotting.

The plot function can display signals using different fixed point formats (and even Double). So, an existential type is required.

data AnySignal = forall b. HasDoubleRepresentation b => AS (Signal b) 

The plotting is using Double at the end (since it is based on PDF).

Once you have displayed some signals, you may want to display some spectrums:

It is done with the code:

spectruma :: Signal Double
(freqR,spectruma) = spectrum samplingFrequency duration (noWindow) mySignalA

freqR is the frequency resolution. It is needed to draw the spectrum.

Note in the previous graphic, only one spectrum is computed with a fixed point FFT (the one with Fixed Int16 3 Saturated). Others are using a double FFT but applied on the fixed point signals.

3. Problems

3.1. Type level arithmetic

For some fixed point computations, you may need to temporarily do the accumulation with higher accuracy. So, I have introduced a specific multiply which is keeping the full accuracy in output:

amul :: (SingI na, SingI nb) => Fixed a na s -> Fixed a nb s -> Fixed b (na + nb) s

It is a complex type. First, I am using two type level naturals : na and nb. But, I need their value during the computation to change the fixed point position. So, I am using the Singleton classes.

Then, the type is specifying that I have more fractional bits at the output using type level arithmetic. Now, the problem is that there is not yet any solver in ghc 7.6.

So, na + nb is never reduced to a natural type. It is causing lots of trouble. In this first version, I only need to compute the type 15 + n and I know n is inferior to 16. So, I have hacked a solution using type family instances:

type instance (15 :: Nat) + (1 :: Nat)  = 16
type instance (15 :: Nat) + (2 :: Nat)  = 17
...

It is ugly and I hope there is a better solution.

3.2. Unbox

I am using Unboxed vectors for the FFT. And I need to use Complex. Unfortunately, Unbox (Complex a) is requiring RealFloat. Which means that to be able to use unboxed vector on Complex of fixed point numbers, I have to define instances RealFloat Int16 etc ... It is not good. I am defining those instances using undefined !

I understand where this constraint is coming from : some Complex functions (like magnitude) are requiring RealFloat. But, I don't like this.

4. Conclusion

So far, things are working like I had imagined. Now, I need to focus on testing the fixed point arithmetic. Then, I'll be able to code some standard signal processing agorithms.

blog comments powered by Disqus