Introduction

If you are into scientific computing, you probably have already heard about Julia, the magical language that aims to be (almost) as fast as C and as easy as MATLAB and Python to write. I have been playing with Julia for more than a year now, and I really like it and recommend checking it out.

The objective of this post is not to teach you Julia. Hopefully, everything used here is so similar to MATLAB or Python you will not need any explanation on what is happening :) The main idea here is to show the basic tools we already have available in DSP.jl (disclaimer: I am one of the developers!) and give you something to play with.

Setup

To be able to run the code in this post (also distributed as an IJulia notebook), you will need to install Julia (of course!) and the following packages:

  • DSP.jl: signal processing functions (filtering, spectrograms/periodograms, window functions...).
  • WAV.jl: loading/writing WAV files.
  • PyPlot.jl: a Julia interface to the Matplotlib plotting library. This is actually a Python library, but there is very little overhead (for examples, arrays are passed to it without making copies).
  • IJulia.jl: where the magic happens :) IJulia notebooks work exactly like IPython notebooks. In fact, they all share the same infrastructure (more on this here).

Installing packages in Julia is easy. Just run the following commands in the Julia REPL:

Pkg.init() # only needed if this is the first time you are dealing with packages
Pkg.update()
Pkg.add("PackageName")

where "PackageName" is the name of the package you want to install (without the .jl suffix).

Loading, visualizing, and processing an audio file


In [263]:
using DSP, WAV, PyPlot # "using" makes the listed modules available for the user, like "import" in other languages

# Loading and plotting an audio signal
s, fs = wavread("test.wav") # Note: multiple outputs do not have to be inside brackets, as in MATLAB

plot(0:1/fs:(length(s)-1)/fs, s)
xlabel("Time [s]")


Out[263]:
PyObject <matplotlib.text.Text object at 0x131433f50>

In [237]:
# PyPlot includes a specgram function, but let's use the native implementation in DSP.jl
# The function below takes a spectrogram with standard parameters for speech (25ms Hanning windows/10ms overlap),
# plots it and returns the spectrogram in case we want to do something with it.

function plot_spectrogram(s, fs)
    S = spectrogram(s[:,1], convert(Int, 25e-3*fs), convert(Int, 10e-3*fs); window=hanning)
    t = time(S)
    f = freq(S)
    imshow(flipud(log10(power(S))), extent=[first(t), last(t), fs*first(f), fs*last(f)], aspect="auto")
    S
end


Out[237]:
plot_spectrogram (generic function with 2 methods)

In [243]:
plot_spectrogram(s, fs)


Out[243]:
Spectrogram{Float64,Frequencies}(201x44 Array{Float64,2}:
 2.68723e-5   1.64834e-7   7.71032e-7   …  1.82563e-7   5.42877e-7 
 0.000627319  0.000113299  8.99384e-6      0.000276147  1.13923e-5 
 0.045069     0.0541711    0.0474623       0.000184796  1.11023e-5 
 0.0970084    0.133594     0.190475        2.16691e-5   4.4943e-6  
 0.00983235   0.0505738    0.043008        4.53554e-6   1.18623e-5 
 0.353044     0.345391     0.25022      …  2.42976e-5   5.33648e-7 
 0.320523     0.623649     0.873488        0.000212918  2.91675e-7 
 0.0111505    0.0616879    0.193093        0.000207865  3.15984e-6 
 0.0011965    0.0043068    0.00121565      2.95437e-6   2.32871e-6 
 0.00115196   0.00162184   0.00619124      2.46006e-5   5.25426e-6 
 0.00245089   0.00131054   0.00192368   …  2.33533e-7   1.06054e-5 
 0.0126481    0.00646363   0.00205837      1.23474e-5   2.57016e-5 
 0.00420449   0.00683724   0.00796684      1.15064e-6   1.0408e-5  
 ⋮                                      ⋱                          
 2.00642e-7   2.1111e-7    1.91358e-7      1.50266e-8   3.09954e-8 
 1.76311e-8   5.29343e-7   1.31138e-7   …  4.86625e-10  1.60534e-8 
 1.05053e-7   4.25714e-7   1.6439e-9       6.39737e-9   1.09893e-8 
 1.0099e-7    1.17022e-7   2.52101e-8      4.55025e-8   1.10105e-8 
 8.65721e-9   1.31638e-7   3.76413e-8      1.43908e-8   4.02047e-9 
 2.55758e-8   5.0852e-8    9.31235e-9      2.11052e-10  8.35175e-11
 8.62341e-9   7.3386e-10   2.20286e-9   …  2.9082e-10   1.20007e-9 
 1.45955e-10  7.599e-10    1.46898e-9      2.36385e-10  8.51917e-10
 1.82281e-10  1.10711e-9   1.87129e-10     3.16577e-11  6.32915e-10
 4.78642e-10  7.6351e-10   4.06273e-10     3.46594e-10  1.12455e-9 
 3.50298e-10  9.89562e-12  8.80213e-10     2.16358e-9   2.98944e-10
 1.96799e-10  5.10707e-11  3.86691e-10  …  1.69894e-9   4.65227e-12,[0.0,0.0025,0.005,0.0075,0.01,0.0125,0.015,0.0175,0.02,0.0225  …  0.4775,0.48,0.4825,0.485,0.4875,0.49,0.4925,0.495,0.4975,0.5],200.0:240.0:10520.0)

In [255]:
# Now let's bandpass the signal to simulate telephone bandwidth and plot its spectrogram again
responsetype = Bandpass(300, 3400; fs=fs)
prototype = Butterworth(8)
telephone_filter = digitalfilter(responsetype, prototype)

# Let's take a look at the filter response now
ω = 0:0.01:pi # variables can have Unicode names! This is typed in the notebook as \omega + tab.
H = freqz(telephone_filter, ω)

plot(fs/2*ω/pi, 20*log10(abs(H)))
xlabel("Frequency [Hz]")
ylabel("Gain [dB]")

# Filtering our signal with the filter
sf = filt(telephone_filter, s)


Out[255]:
10891x1 Array{Float64,2}:
 -6.8193e-5  
 -0.000725265
 -0.00346986 
 -0.00971748 
 -0.0170789  
 -0.0176696  
 -0.00493417 
  0.0162909  
  0.0314648  
  0.0312873  
  0.0212183  
  0.0129548  
  0.0101269  
  ⋮          
  0.000166085
  0.000217562
  0.000477398
  0.000685943
  0.000629586
  0.000374882
  0.000200836
  0.000299072
  0.00056069 
  0.000698629
  0.000568364
  0.000302171

In [252]:
# Just to be sure the filter has done the right thing, let's see the spectrogram of the filtered signal
plot_spectrogram(sf, fs)


Out[252]:
Spectrogram{Float64,Frequencies}(201x44 Array{Float64,2}:
 3.90714e-9   1.11785e-8   1.00194e-8   …  6.89713e-12  1.04341e-12
 1.31113e-8   3.61911e-8   2.96357e-8      1.29152e-9   6.52004e-12
 6.86915e-8   1.19946e-7   4.45702e-8      7.91152e-9   3.45687e-11
 2.98211e-7   1.10841e-6   6.73752e-7      3.04682e-8   8.11075e-11
 4.088e-5     3.66121e-5   4.72577e-6      1.23609e-7   3.11901e-9 
 0.00118734   0.00223387   0.00352084   …  2.30944e-6   9.68379e-8 
 0.00187092   0.00342642   0.0103947       2.75304e-5   2.03604e-6 
 0.00200812   0.000550608  0.00209669      1.24224e-5   7.90119e-6 
 0.00251585   0.00124552   0.00189303      6.11721e-5   1.89558e-5 
 7.80461e-5   0.00169565   0.00542889      6.92665e-5   2.47635e-5 
 0.00444074   0.00047806   0.00173713   …  1.11671e-5   2.77241e-5 
 0.0113532    0.0084156    0.00157983      5.68708e-6   3.24141e-5 
 0.00335498   0.00640711   0.0074975       2.77168e-6   8.22006e-6 
 ⋮                                      ⋱                          
 7.79231e-18  1.22182e-17  4.46617e-17     1.01478e-17  7.68596e-19
 7.55378e-18  1.14103e-17  4.52902e-17  …  8.4437e-18   7.37437e-19
 7.3406e-18   1.06961e-17  4.58588e-17     6.91511e-18  7.09632e-19
 7.15102e-18  1.00715e-17  4.63662e-17     5.55765e-18  6.85018e-19
 6.9848e-18   9.5309e-18   4.68132e-17     4.36791e-18  6.63537e-19
 6.84153e-18  9.07052e-18  4.71997e-17     3.34281e-18  6.4508e-19 
 6.72087e-18  8.68679e-18  4.75262e-17  …  2.47977e-18  6.29583e-19
 6.62254e-18  8.37677e-18  4.7793e-17      1.77658e-18  6.16984e-19
 6.5463e-18   8.13808e-18  4.80002e-17     1.23149e-18  6.07235e-19
 6.49196e-18  7.96889e-18  4.8148e-17      8.43117e-19  6.00299e-19
 6.45942e-18  7.8679e-18   4.82367e-17     6.10481e-19  5.96147e-19
 3.22429e-18  3.91716e-18  2.41331e-17  …  2.665e-19    2.97383e-19,[0.0,0.0025,0.005,0.0075,0.01,0.0125,0.015,0.0175,0.02,0.0225  …  0.4775,0.48,0.4825,0.485,0.4875,0.49,0.4925,0.495,0.4975,0.5],200.0:240.0:10520.0)

In [260]:
# IJulia still does not support playing audio, so we include some code to do it.
include("AudioDisplay.jl")

using AudioDisplay

audioplayer(s, fs)
audioplayer(sf, fs)


Warning: replacing module AudioDisplay
Warning: using AudioDisplay.audioplayer in module Main conflicts with an existing identifier.

(Just a sidenote: the AudioIO.jl module has more advanced functions for loading and playing audio, but in this example I wanted to use the notebook for all multimedia rendering.)

Where to go from here?

I know this is all really basic, but Julia has much more to offer. There are currently 429 packages (and counting) available to install via the package manager and much more in development on GitHub. Package development has been especially strong in the statistics/machine learning and optimization domains.

If you want to learn more about the language, check the official documentation and the learning resources at the website. There is also a very active mailing list for users.

I sincerely hope you find Julia interesting and worth of your time!