DSP 0.01: A tutorial

[zer0w1ng]

IIR Filter

Same with our previous lesson, the FIR filter, IIR filter also provides filtering our signal to our desired response.
IIR filter provides less computational complexity and time (MIPS) than its counterpart the FIR filter.
It has some disadvantages though, like stability, where it may cause oscillations if the coefficients are wrongly made.

IIR employs feedback where the previous output samples are used to somewhat influence the next output sample.

--

The basic response of the IIR filter, H(z), is computed as:

Code: [Select]

       Y(z)   b0*z^2 + b1*z + b2
H(z) = ---- = ------------------   (eq. 1)
       X(z)   a0*z^2 + a1*z + a2

Where Y(z) and X(z) are the output and input samples, respectively.

Simplifying:

Code: [Select]

Y(z)   b0 + b1*z^-1 + b2*z^-2
---- = ----------------------
X(z)   a0 + a1*z^-1 + a2*z^-2

Y(z)   (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
---- = -------------------------------------
X(z)   1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

Y(z) * ( 1 + (a1/a0)*z^-1 + (a2/a0)*z^-2) = X(z) * ((b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2)  (eq. 2)

z raised to the negative power theoritically is just the delayed sample.
Meaning, z^-1 is the previous sample and z^-2 is the 2nd to the last sample.

Thus our last equation (eq. 2) will become:

Code: [Select]

y[n] + (a1/a0)*y[n-1] + (a2/a0)*y[n-2] = (b0/a0)x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]

y[n] = (b0/a0)x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2]   (eq. 3)

Equation 3 is now the basic IIR filter.

(b0/a0), (b1/a0) and (b2/a0) are the INPUT coefficients "x(n)".
While (a1/a0) and (a2/a0) are the OUTPUT coefficents "y[n]".

So manually computation for IIR is as follows:

Code: [Select]

x_coef = [11 22 33]
y_coef = [88 99]
signal = [1 2 3 4 5]

output[1] = 0
output[2] = 0
output[3] = 11*3 + 22*2 + 33*1 - 0*88 - 0*99           = 110
output[4] = 11*4 + 22*3 + 33*2 - 110*88 - 0*99         = -9504
output[5] = 11*5 + 22*4 + 33*3 - (-9504*88) - 110*99   = 825704


In SCILAB it is also just convolution (multiply and add operation only):

Code: [Select]

//DSP Tutorial
//IIR Filter
//By: Regulus Berdin

clear;

function output = compute_iir(signal, x_coef, y_coef)
   N = length(signal);
   YN = length(y_coef);
   XN = length(x_coef);
   
   output = zeros(1:XN);
   
   for i = YN+1:N
      xframe = signal(i:-1:i-(XN-1));
      yframe = output(i-1:-1:i-YN);

      x = xframe .* x_coef;
      y = yframe .* y_coef;
     
      output(i) = sum(x) - sum(y);
   end
endfunction   

//--------

x_coef = [11 22 33];
y_coef = [88 99];
signal = [1 2 3 4 5];

output = compute_iir(signal, x_coef, y_coef);

output


--

There are many ways in computing for the 2 groups of coefficents.
One is using Winfilter or using formulas.

The best reference I had for formulas for IIR filters is this site:
http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

It covers all the kinds of filters and can specify the bandwidth or the slope of the filter.

--

For our next examples. We will be dealing with HiFi audio.
We need to save a sample of music and save it in WAV.
WAV is uncompressed so this is quite large.
We will just limit our sample music to 10 seconds or less for SCILAB may also choke up with the large memory requirement.

But in reality, using DSP chips, techniques are employed, such as double buffering, to ensure continous and real-time results are experienced for its output. 

DSP chips also employs optimized computaion techniques to minimize delay.
Such as the single cycle MAC instruction, no-overhead loops, reverse addressing, simultaneous execute, fetch and put and other optimizations on its hardware thus it is considered a specialized chip compared to an ordinary MCU or processor.

--

For our next example, I had cut the first 10 seconds of Billy Jean from MJ as our reference signal.
Download it here: http://www.fileden.com/files/2009/11/25/2664828/music.wav
You could also use other music files but just limit it to 10 seconds to lessen the cruelty to your PC.

['yus]

yeah.. nakadownload na rin ako nung portable SciLab 5..
habol muna ako sa discussion... pag turn ko ng page 3, nagulat ako, biglang bumilis ang discussion.

  gawa din ako ng portable python(Eric) with scipy

 thanks sir zer0w1ng

[zer0w1ng]

You're welcome.
Pause muna ako. Baka sa lunes na ang susunod na example.

[♪ ♫ ♫ ♪ ♪ ♪ ♫ ♫ ♫]

akala ko dsp sim itong thread na to, dsp pala talaga ang dinidiscuss, thanks sir zero1ng, hindi ako makabwelo sa mcu dahil malabo pa sa akin ang dsp, ang ganda ng timing nyo a!

charge kita palagi para hindi maputol 

 

[0b00000111]

pag turn ko ng page 3, nagulat ako, biglang bumilis ang discussion.

nice thread!

[e.novacek]

This is a good start zer0w1ng.. it is useful and timely relevant especially to engineering students, professionals, hobbyist and learners as well..

but since the usage of Scilab was emphasize then I would like to request to change the title if ever.. my suggestion is "Practical DSP Tutorial using Scilab" ..

I would like to request to admin electron (as well as with our moderators) to make the tutorial "sticky" so as newbies, guests and others may find this thread easily..

We are looking forward to continue doing the tutorial and finishing it as well.. I know sometimes this is difficult but it's one way to serve and educating other's..

+1..

I hope "yus" (..with his python) and other's as well may do start doing technical tutorial..

[waterbender]

I so wished na sa undergrad ko ito na experience.. 

[0b00000111]

me too... sa batch namin, kami ang pinaka unang nakatikim ng dsp na subject. kaya lahat nangangapa pati prof.

[waterbender]

me too... sa batch namin, kami ang pinaka unang nakatikim ng dsp na subject. kaya lahat nangangapa pati prof.

buti nga sa inyo sir 7 may DSP subject.

[Google]

sa amin naman alam na ng prof yung DSP pero ang problema sa software naman

[marcelino]

sa amin naman alam na ng prof yung DSP pero ang problema sa software naman

normally, sa mga universities, matlab ang ginagamit... kaso mahal talga yan. eto nga scilab or octave pwede naman...

[ef - el - ay - pee]

kabibili lang namin ng matlab sa lab for 5 concurrent users, and swear ang mahal niya ... pero sa dsp grad subject required na octave or scilab ang gamit so i think ok naman na ang open-source languages na ito ... tsaka active naman ang community na gumawa ng functions sa octave at scilab ...

cant wait for the next lesson.

[zer0w1ng]

PRACTICAL EXAMPLE FOR IIR FILTER

We had discussed IIR filtering on my last post.
Due to its minimal MAC operation compared to FIR, this filter is commonly used in new HiFi audio.
Mainly for canned equalization we see such as ROCK, POP, CLASSIC, SOFT, JAZZ, etc.
Nevertheless, this is commonly found in digital parametric equalizers where one can adjust the Q and the frequency of the band to set.

--

For this example, I will give an example for bass boosting using DSP.
This is pretty simple.  It is just using a shelving filter.

A shelving filter either increases or decreases one part of the band, either low or high but does not change the other band.
In our case, we will use a low shelving filter.

We will set our gain to positive to increase the bass and set our cut-off frequency to 400Hz.
A gain of 9 dB should give us a HUGE bass boost.

We will use the "Cookbook formulae for audio EQ biquad filter coefficients" by Robert Bristow-Johnson for our IIR coefficients computation:
http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Using MathCad or any formula solver (even MsExcel could be used), compute the following:

Code: [Select]

fo = 400
fs = 44100
S = 1

A = 10^(9/20)
w0 = 2 * 3.1416 * fo / fs

alpha = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 )

b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   )
b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
a1 =   -2*( (A-1) + (A+1)*cos(w0)                   )
a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha

x0 = b0/a0
x1 = b1/a0
x2 = b2/a0

y0 = a1/a0
y1 = a2/a0

The value of S (shelf-slope) equal to "1" is chosen to be so the filter's response is as steep as it can be.

The resulting values x0, X1 and X2 are the coefficients for the input.
  input coefficients = (+1.04457887924806,-1.94808417101059,+0.91244684520764)

WHile y0 and y1 are for the coefficients for the output.
  output coefficients = (-1.95199211029050,+0.95311778517579)

---

Now that we have our coefficients for our filter, we will use our previous hifi sound sample.
We will only be using 2 seconds of it.  
My computer, even a core2quad with 4Gb memory, freezed when I process the whole 10 seconds.
Scilab sucks on high memory requirement programs. Python with scipy/pylab should be nicer.
Anyway, c'est la vie, we will continue and only use 2 seconds of it.

--

To run our filter, we will use the our iir computation function discussed on the last topic.
We will cut the clip at between 5 and 7 seconds for our 2 second sample.
Then we separate the 2 channels (left and right) and process them separately.

After the 2 processed outputs are combined to become stereo again.

A demo matrix is then created by having the original clip and the processed clip concatenated.

The following is our resulting bass boost Scilab program:

Code: [Select]

//DSP Tutorial
//IIR Filter
//By: Regulus Berdin

clear;

function output = compute_iir(signal, x_coef, y_coef)
   N = length(signal);
   YN = length(y_coef);
   XN = length(x_coef);
  
   temp = zeros(1:XN);
   win1 = progressionbar("Working.");
  
   for i = YN+1:N
      xframe = signal(i:-1:i-(XN-1));
      yframe = temp(i-1:-1:i-YN);
      temp(i) = sum(xframe * x_coef') - sum(yframe * y_coef');
   end
   output = temp;
   winclose(win1)
  
endfunction  

//--------

x_coef = [+1.04457887924806 -1.94808417101059 +0.91244684520764];
y_coef = [-1.95199211029050 +0.95311778517579];

clip_start = 5;
clip_end = 7;
[signal,fs,bits] = wavread('c:/music.wav',[44100*clip_start, 44100*clip_end]);

//scale so not to clip
signal = signal * 0.2;

left_ch = compute_iir(signal(1,:), x_coef, y_coef);
right_ch = compute_iir(signal(2,:), x_coef, y_coef);

stereo_out = [left_ch; right_ch];

demo = [signal stereo_out];

printf("\nPlaying original and the processed sound clip....\n");
playsnd(demo,fs);

You should hear a very big difference in the bass on the first 2 seconds with the next 2 seconds but not changing the higher band.

--

Try experimenting with different filters, such as the high shelf.
Combine it with the low by passing the processed output into another IIR filter.

You can try also using the IIR filter on our previous FIR examples.

Please show your programs and describe the results so you could share with others what you have learned and experienced.

--

Now you have learned all the basic time-domain filtering.
We shall proceed next with the the frequency domain and creating spectrographs.

But before that, I will give another time-domain process example.
This is the STEREO WIDENING or 3D audio effect.
This is not using filters but an just an enhancement.

This is pretty much easy once you know its theory.
http://en.wikipedia.org/wiki/3D_audio_effect

I will give +points to those who could do this in advance.

[anthonydp]

Here's a try for the stereo widening exercise.  I can't here much effect with the 30% attenuation but with 50% the center beat lessens while the cymbal sounds at the right speaker remains loud.

Stereo Widening:

Another way of looking at this same effect, without extrapolating a center and side signal from the left and right signals, is to simply add the left signal, slightly attenuated and phase inverted, into the right channel and vice-versa.

Code: [Select]

//prevent stack memory error:

stacksize('max');

//load the 10-second Billy Jean intro:

[signal,fs,bits] = wavread('c:/music.wav');

//separate the left and right channel:

left_ch = signal(1,:);
right_ch = signal(2,:);

//define our attenuation factor:
//try:
// atten = 0.50;
 
atten = 0.30;

//widen the left and right channel:

wide_left_ch = left_ch - atten * right_ch;
wide_right_ch = right_ch - atten * left_ch;

//build the new stereo signal from the widened channels

wide_stereo = [wide_left_ch; wide_right_ch];

//let's hear it:

playsnd(wide_stereo,fs);



[anthonydp]

omg! There's a much wider stereo effect if we delay the right channel for 10msec upon combining the left and right signal.

Code: [Select]

//define a delay function
//fs -- the sampling rate from wavread
//tsec -- number of seconds of delay (30msec = 0.030sec)

function output = delay(signal, fs, tsec)
  n = length(signal);
  timedelay = fs * tsec; //number of samples to skip at the beginning
  temp = zeros(1:n);
  for i = timedelay:n
    temp(i) = signal(i-timedelay+1);
  end
  output = temp;
endfunction

//prevent stack memory error:

stacksize('max');

//load the 10-second Billy Jean intro:

[signal,fs,bits] = wavread('c:/music.wav');

//separate the left and right channel:

left_ch = signal(1,:);
right_ch = signal(2,:);

//define our attenuation factor:
//try:
// atten = 0.50;
 
atten = 0.30;

//widen the left and right channel:

wide_left_ch = left_ch - atten * right_ch;
wide_right_ch = right_ch - atten * left_ch;

//build the new stereo signal from the widened channels
//but now we delay the right channel for 10msec

wide_stereo = [wide_left_ch; delay(wide_right_ch, fs, 0.01)];

//let's hear it:

playsnd([signal wide_stereo],fs);



[zer0w1ng]

+1 to anthony.

[zer0w1ng]

BTW the delay could just be:
delay_samples = int(0.01 * fs);
delayed_right_ch = [ zeros(1:delay_samples)  right_ch(1:length(right_ch) - delay_samples) ];

It is faster to slice the array than to make one by one.

[anthonydp]

BTW the delay could just be:
delay_samples = int(0.01 * fs);
delayed_right_ch = [ zeros(1:delay_samples)  right_ch(1:length(right_ch) - delay_samples) ];

It is faster to slice the array than to make one by one.

oo nga, very nice! thanks.
and thanks for the point.

[rdpzycho]

normally, sa mga universities, matlab ang ginagamit... kaso mahal talga yan. eto nga scilab or octave pwede naman...

alam ko sis sa SLU Scilab ang gamit nila.

mahal nga ang Matlab. naalala ko nun 'yung dongle na may key for installation ng Matlab sa lab namin naiwan sa isang computer. mukha lang siyang extender ng parallel port, muntik na namin maitapon.

[tiktak]

scilab nga po gamit namin sa SLU dati dati matlab pero licensing issues kaya scilab na sayang lang ayaw nila ipagamit mga ADI na blackfin dun kinakapa padin kasi ng prof sira na din kasi mga TI na board