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.
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:
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:
n is using type level numeral and so you need ghc 7.6 with data type promotion. So n has kind Nat
s has the kind Saturation and can be the type Saturated or Unsaturated.
So, an example of fixed point type could be:
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
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 = 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 = 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:
, 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.
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:
(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.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:
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) + (2 :: Nat) = 17
It is ugly and I hope there is a better solution.
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.
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.