1.       Write a program (in C/C++) to design a FIR low pass filter. The program should ask for the number of taps, sampling frequency, the cut-off frequency and whether windowing is to be used. Choose a suitable window.

The C++ programming language was chosen to complete the objectives outlined for this assignment. The FIR low pass filter was created through research of similar implementations in books, notes and web based examples. The program was compiled using Visual Studio.

F.I.R. filters (Finite Impulse Response) are implemented using a finite number “n” delay taps on a delay line and “n” computation coefficients to compute the algorithm (filter) function.  They can offer shape, factor accuracy and stability equivalent to very high-order linear active filters that cannot be achieved in the analogue domain.

FIR filters can create transfer functions that have no equivalent in linear circuit technology. FIR filters are formed with only the equivalent of zeros in the linear domain. This means that the taps depress or push down the amplitude of the transfer function. The amount of depression for each tap depends upon the value of the multiplier coefficient. Hence, the total number of taps determines the “steepness” of the slope. More taps increase the steepness of the filter roll-off while increasing calculation time (delay) and for high order filters, limiting bandwidth. There is tradeoff between phase delay and filter precision when designing FIR filters.

Filter Design:

For my filter, I have utilized the Windowing and the Hanning Technique.

Windowed Filter:

The simplest technique is known as “Windowed” filters. This technique is based on designing a filter using well-known frequency domain transition functions called “windows”. The use of windows often involves a choice of the lesser of two evils. Some windows, such as the Rectangular, yield fast roll-off in the frequency domain, but have limited attenuation in the stop-band along with poor group delay characteristics. Other windows like the Blackman, have better stop-band attenuation and group delay, but have a wide transition-band (the band-width between the corner frequency and the frequency attenuation floor). Windowed filters are easy to use, are scalable (give the same results no matter what the corner frequency is) and can be computed on-the-fly by the DSP.

Rectangular and Hanning Windowing

Rectangular windowing is a simple multiplication by 0 or 1, the impulse response is cut to zero at the point the window stops. In this example the window is the size of the entire filter effectively making no difference to the output.

w(n) = 1

Hanning Window:

The “Raised Cosine” window uses a similar technique although the coefficient is calculated differently.

Calculate Coeefficient

To determine the coefficients and fill the coefficients array, requires the “Nyquist Frequency” of the filter.

nyquist frequency

H(Θ) is the frequency response as a function of the Nyquist frequency. Using a discrete inverse Fourier transform the coefficients for a low pass filter can be described as:

Low Pass Filter Coefficients

The code for my filter is contained in Appendix A.

filter.cpp is the main application source file.

filter.h is the header file.

main.cpp has the entry point for the program. In that object of filter, a class is created

and filtering is called to start the filtering.

2.       Use your program to design a filter (i)with and (ii)without windowing, whose number of taps is equal to the last two digits of your id number. The sampling frequency is to be 8000Hz and cut off frequency is 200Hz.

For this part, I used 101 as my number of taps, my id ends in 00 so I added 101.The graphs are plotted using Microsoft Excel. The function CalCoeff(), calculates the coefficients and these are displayed using the display() function. The code for my filter is contained in Appendix A.

(i) With Hanning Window

Phase Response:

Phase Response with windowing

Figure 1 – Phase Response with windowing

Appendix B gives the values used to calculate this graph.

Amplitude Response:

Purple indicates the use of Hanning windowing.

Blue indicates without the use of windowing.

Amplitude Response with and without windowing

Figure 2 – Amplitude Response with and without windowing

Appendix B gives the values used to calculate this graph.

(ii) Without Windowing:

Phase Response:

Phase Response without Windowing

Figure 3 – Phase Response without Windowing

Appendix B gives the values used to calculate this graph.

Amplitude Response:

Please see Figure 2.

Appendix B gives the values used to calculate this graph.

3.       Write a program (in C/C++) to calculate the frequency response (Amplitude & Phase) of a digital filter and use it to plot the frequency response of the filters you designed in 2. Compare the responses (Amplitude & Phase).

For this part, I added FreqResponse() function to the program which calculates the Amplitude & Phase response. The calculated phase response and amplitude response values are stored in phase.txt and amp.txt respectively. The code for my filter is contained in Appendix A.

Phase Response of Filter designed in Part 2

Figure 4 – Phase Response of Filter designed in Part 2

Amplitude Response of Filter designed in Part 2

Figure 5 – Amplitude Response of Filter designed in Part 2

4.       Make a file of about 400 samples [at 8kHz] of your own speech waveform.[try to get the transition between a fricative (like "sh") and a vowel (like "a")]. Plot the original waveform, process it through your filter and plot the resultant output. Compare the input and output signals. What is the time shift between the input and output? Why?

I created a file at 8kHz using the word “shark” to create the waveform beneath. The results used to calculate Figure 6 and 7 are contained in Appendix C. There is a time shift between input & output signal, as the speech signal passes through the filter. The filtered output is longer than the input signal by no of taps chosen – in our case 100 (101 – 1). So, time shift = 100.

I changed the CalCoeff() function of my previous program to write coefficients to a file.

void Filter::CalCoeff()

{

//Allocate memory for coefficients

p_mcoeff = new double[m_tap];

ofstream coeffFile;

coeffFile.open (“C:\coeff.txt”);

int i=1,count;

//To keep our filter causal, shift the centre of it to the middle of taps.

if(m_tap%2 == 0)

{

count = (m_tap – 1)>>1;

}

else

{

count = m_tap >> 1;

}

//Calculate centre value of coefficient.

p_mcoeff[count] = m_nqFreq/PI;

//Calculate remaining values of coefficient.

while(i<=count)

{

p_mcoeff[count-i]=sin(i*m_nqFreq)/PI/i;

p_mcoeff[count+i]=p_mcoeff[count-i];

i++;

}

coeffFile<<m_tap<<” “;

for(i=0;i<m_tap;i++)

{

Then I processed the coefficients of my filter and my waveform through the new filter. This is contained in Appendix B. The same coefficients are used in Part 5. It utilizes the same header file. I edited the main file to prompt the user to input files. This is contained in Appendix B. This filter performs 4 main tasks:

1.         Reads the coefficients of filter from the specified .txt file.

2.         Reads the samples of the speech signal from the specified .wav file.

3.         Convolution of coefficients & samples.

4.         Writes the output of the filtered signal in a file.

Speech Waveform

Figure 6 –  Speech Waveform

Output Waveform

Figure 7 –  Output Waveform

5.       Using Matlab, plot the spectrum of the input and output [of the vowel part of the] speech signals. Comment.

From the graphs, it can be seen that both signals have different levels of power. The filtered output can be seen in Figure 8.

The Matlab code used to complete this:

sound_data = wavread(‘C:\temp\part_4_voice.wav’);

filter_coeff = [0.00636938;0.00641538;0.00630221;0.00602589;0.00558612;0.00498633;0.00423384;0.00333978;0.00231904;0.00119006;-2.53607e-005;-0.00130248;-0.00261392;-0.00393023;-0.00522029;-0.006452;-0.00759278;-0.00861027;-0.00947294;-0.0101507;-0.0106157;-0.0108426;-0.0108094;-0.0104981;-0.00989464;-0.00898978;-0.00777917;-0.00626361;-0.00444919;-0.00234732;2.53607e-005;0.00264715;0.00549151;0.00852747;0.0117201;0.0150308;0.0184184;0.0218392;0.0252482;0.0285994;0.0318471;0.0349461;0.0378527;0.0405256;0.0429265;0.0450207;0.0467777;0.0481717;0.0491824;0.0497948;0.05;0.0497948;0.0491824;0.0481717;0.0467777;0.0450207;0.0429265;0.0405256;0.0378527;0.0349461;0.0318471;0.0285994;0.0252482;0.0218392;0.0184184;0.0150308;0.0117201;0.00852747;0.00549151;0.00264715;2.53607e-005;-0.00234732;-0.00444919;-0.00626361;-0.00777917;-0.00898978;-0.00989464;-0.0104981;-0.0108094;-0.0108426;-0.0106157;-0.0101507;-0.00947294;-0.00861027;-0.00759278;-0.006452;-0.00522029;-0.00393023;-0.00261392;-0.00130248;-2.53607e-005;0.00119006;0.00231904;0.00333978;0.00423384;0.00498633;0.00558612;0.00602589;0.00630221;0.00641538;0.00636938;];

filtered_output = conv(filter_coeff,sound_data);

fft_output = abs(fft(filtered_output));

fft_input = abs(fft(sound_data));

subplot(2,1,1), plot(fft_input,’red’);

xlabel(‘frequency’);

ylabel(‘Amplitude’);

title(‘Specturm of input speech signal’);

subplot(2,1,2),plot(fft_output,’blue’);

xlabel(‘frequency’);

ylabel(‘Amplitude’);

title(‘Specturm of output speech signal’);

Input and Output Waveform

Figure 8 –  Input and Output Waveform

6.       Construct a real – time implementation of your filter on the C31 DSK. Record a few seconds o your speech.

To complete this part, I had to write a C31 Assembler version of my filter. I used RFIR.ASM as a guide.

There were 3 main changes I made to this:

  • TA        .set     48.828

This had to be changed to 48.828 as the original sampling frequency was set too high.

  • The Data buffer was changed to 128 bits(2^8) as I used 101 taps for my filter.
  • The coefficients were changed to the values of my filter.

;———————————————————

; FIR.ASM

; Keith Larson

; TMS320 DSP Applications

; (C) Copyright 1995,1996

; Texas Instruments Incorporated

;

; This is unsupported freeware with no implied warranties or

; liabilities.  See the disclaimer document for details

;

; The FIR filter code used in this example is taken

; from the TMS320C3x Users Guide.  The AIC setup and

; control is designed to work with the TMS320C31 DSK

;

; This example can either be loaded and run from the debugger

; or by directly loading and running from DSK3LOAD

;———————————————————

;

;     How to change the sampling frequency.

;

;     [Ronan Scaife, DCU, October 17th 2001]

;

;     By changing either TA/TB and RA/RB in the AIC,

;     OR by changing the figures for the C31 timer

;     see the section of code labelled “ST-STUB”

;     below, you can change the sampling frequency.

;     The default settings in Keith Larsen’s

;     original version give fS ~ 48,828 s/s.

;     For more interesting effects, I have modified

;     TA and RA to 16 to reduce fS to 24,414 s/s.

;

;     Edited By Sameer Kumar

;———————————————————

;

; Define constants used by program    ;

;TA        .set     8                  ; AIC timing register values

TA        .set     48.828                 ; AIC timing register values/this value must be changed as the original sampling frequency is too high

TB      .set         16             ;

;RA        .set     8                  ;

RA        .set     16                  ;

RB      .set         16             ;

GIE       .set     0×2000             ; This bit in ST turns on interrupts

.include “C3XMMRS.ASM”      ;

;

;     change start addr for WinDSK [RS 8-x-01]

;

;          .start   “AICTEST”,0×809820 ; Start assembling here

.start   “AICTEST”,0×809900 ; Start assembling here

.sect    “AICTEST”          ;

;- – - – - – - – - – - – - – - – - – - – - – - – - – - – - – -

; 128 bit float data buffer of incoming data from the AIC

; The location in memory must be on a 2^N boundary and the

; size of the coefficeint table must be the same as the data

; Taps for my filter = 101, therefore 128 bit data buffer must be utilised

;- – - – - – - – - – - – - – - – - – - – - – - – - – - – - – -

ADC_recv     .float    0.0          ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0           ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0           ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

.float    0.0            ;

;- – - – - – - – - – - – - – - – - – - – - – - – -

; FIR filter coefficients

;- – - – - – - – - – - – - – - – - – - – - – - – -

FIR_coef    .float    0.00636938

.float       0.00641538

.float       0.00630221

.float       0.00602589

.float       0.00558612

.float      0.00498633

.float       0.00423384

.float       0.00333978

.float       0.00231904

.float       0.00119006

.float       -2.53607e-005

.float       -0.00130248

.float       -0.00261392

.float       -0.00393023

.float       -0.00522029

.float       -0.006452

.float       -0.00759278

.float       -0.00861027

.float       -0.00947294

.float       -0.0101507

.float       -0.0106157

.float       -0.0108426

.float       -0.0108094

.float       -0.0104981

.float       -0.00989464

.float       -0.00898978

.float       -0.00777917

.float       -0.00626361

.float       -0.00444919

.float       -0.00234732

.float       2.53607e-005

.float       0.00264715

.float       0.00549151

.float       0.00852747

.float       0.0117201

.float       0.0150308

.float       0.0184184

.float       0.0218392

.float       0.0252482

.float       0.0285994

.float       0.0318471

.float       0.0349461

.float       0.0378527

.float       0.0405256

.float       0.0429265

.float       0.0450207

.float       0.0467777

.float       0.0481717

.float       0.0491824

.float      0.0497948

.float      0.05

.float      0.0497948

.float      0.0491824

.float      0.0481717

.float      0.0467777

.float      0.0450207

.float      0.0429265

.float      0.0405256

.float      0.0378527

.float      0.0349461

.float      0.0318471

.float      0.0285994

.float      0.0252482

.float      0.0218392

.float      0.0184184

.float      0.0150308

.float      0.0117201

.float      0.00852747

.float      0.00549151

.float      0.00264715

.float      2.53607e-005

.float      -0.00234732

.float      -0.00444919

.float      -0.00626361

.float      -0.00777917

.float      -0.00898978

.float      -0.00989464

.float      -0.0104981

.float      -0.0108094

.float      -0.0108426

.float      -0.0106157

.float      -0.0101507

.float      -0.00947294

.float      -0.00861027

.float      -0.00759278

.float      -0.006452

.float      -0.00522029

.float      -0.00393023

.float      -0.00261392

.float      -0.00130248

.float      -2.53607e-005

.float      0.00119006

.float      0.00231904

.float      0.00333978

.float      0.00423384

.float      0.00498633

.float      0.00558612

.float      0.00602589

.float      0.00630221

.float      0.00641538

.float      0.00636938       ;  FIR filter coefficients

END_coef

;- – - – - – - – - – - – - – - – - – -

BufSz      .set      END_coef – FIR_coef

;- – - – - – - – - – - – - – - – - – -

SIZE         .word    BufSz           ; Size of filter

ADC_first    .word    ADC_recv

ADC_end      .word    FIR_coef

ADC_last     .word    ADC_recv

FIR_coefx    .word    FIR_coef

;————————————

; Define some constant storage data

;————————————

A_REG        .word  (TA<<9)+(RA<<2)+0 ; A registers

B_REG        .word  (TB<<9)+(RB<<2)+2 ; B registers

C_REG      .word  00000011b         ; control

S0_gctrl_val .word  0x0E970300        ; Serial port control register values

S0_xctrl_val .word  0×00000111        ;

S0_rctrl_val .word  0×00000111        ;

;****************************************************

; Begin main code loop here

;****************************************************

main      or    GIE,ST          ; Turn on INTS

ldi   0xF4,IE         ; Enable XINT/RINT/INT2

b     main            ; Do it again!

;——————————-

DAC2      push  ST              ; DAC Interrupt service routine

push  R0              ;

pushf R0              ;

push  R2              ;

pushf R2              ;

push  AR0             ;

push  AR1             ;

ldi   @ADC_last,AR1   ;

ldi   @FIR_coefx,AR0  ;

ldi   @SIZE,BK

FIR       mpyf3 *AR0++,*AR1++(1)%,R0

ldf   0.0,R2

ldi   @SIZE,RC

subi  2,RC

rptb  FIR2

mpyf3 *AR0++,*AR1++(1)%,R0

FIR2  ||  addf3 R0,R2,R2

addf  R2,R0

fix   R0,R0

andn  3,R0            ;

sti   R0,@S0_xdata    ; Output the new DAC value

pop   AR1             ;

pop   AR0             ;

popf  R2              ;

pop   R2              ;

popf  R0              ;

pop   R0              ;

pop   ST              ;

reti                  ;

;——————————-

ADC2      push  ST              ;

push  R3              ;

pushf R3              ;

push  AR0             ;

ldi   @S0_rdata,R3    ;

lsh   16,R3

ash   -16,R3

ldi   @ADC_last,AR0

float R3,R3

stf   R3,*AR0++

cmpi  @ADC_end,AR0

ldige @ADC_first,AR0

sti   AR0,@ADC_last

pop   AR0

popf  R3              ;

pop   R3              ;

pop   ST              ;

reti                  ;

;*****************************************************;

; The startup stub is used during initialization only ;

; and can be safely overwritten by the stack or data  ;

;*****************************************************;

.entry   ST_STUB      ; Debugger starts here

ST_STUB   ldp   T0_ctrl         ; Use kernel data page and stack

ldi   @stack,SP

ldi   0,R0            ; Halt TIM0 & TIM1

sti   R0,@T0_ctrl     ;

sti   R0,@T0_count    ; Set counts to 0

ldi   1,R0            ; Set periods to 1

sti   R0,@T0_prd      ;

ldi   0x2C1,R0        ; Restart both timers

sti   R0,@T0_ctrl     ;

;———————

ldi   @S0_xctrl_val,R0;

sti   R0,@S0_xctrl    ; transmit control

ldi   @S0_rctrl_val,R0;

sti   R0,@S0_rctrl    ; receive control

ldi   0,R0            ;

sti   R0,@S0_xdata    ; DXR data value

ldi   @S0_gctrl_val,R0; Setup serial port

sti   R0,@S0_gctrl    ; global control

;======================================================;

; This section of code initializes the AIC             ;

;======================================================;

AIC_INIT  LDI   0×10,IE         ; Enable only XINT interrupt

andn  0×34,IF         ;

ldi   0,R0            ;

sti   R0,@S0_xdata    ;

RPTS  0×040           ;

LDI   2,IOF           ; XF0=0 resets AIC

rpts  0×40            ;

LDI   6,IOF           ; XF0=1 runs AIC

;———————

ldi   @C_REG,R0       ; Setup control register

call  prog_AIC        ;

ldi   0xfffc  ,R0     ; Program the AIC to be real slow

call  prog_AIC        ;

ldi   0xfffc|2,R0     ;

call  prog_AIC        ;

ldi   @B_REG,R0       ; Bump up the Fs to final rate

call  prog_AIC        ; (smallest divisor should be last)

ldi   @A_REG,R0       ;

call  prog_AIC        ;

b     main            ; the DRR before going to the main loop

;——————————-

prog_AIC  ldi   @S0_xdata,R1    ; Use original DXR data during 2 ndy

sti   R1,@S0_xdata    ;

idle

ldi   @S0_xdata,R1    ; Use original DXR data during 2 ndy

or    3,R1            ; Request 2 ndy XMIT

sti   R1,@S0_xdata    ;

idle                  ;

sti   R0,@S0_xdata    ; Send register value

idle                  ;

andn  3,R1            ;

sti   R1,@S0_xdata    ; Leave with original safe value in DXR

;———————

ldi   @S0_rdata,R0    ; Fix the receiver underrun by reading

rets                  ;

stack     .word   $             ; Put stack here

;****************************************************;

; Install the XINT/RINT ISR handler directly into    ;

; the vector RAM location it will be used for        ;

;****************************************************;

.start   “SP0VECTS”,0x809FC5

.sect    “SP0VECTS”

B        DAC2         ; XINT0

B        ADC2         ; RINT0

7.       Play back the original speech and the signal  filtered by your filter. Listen carefully to both and comment on what you hear.

The signal filtered by my filter is muffled in comparison to the original speech. The ‘s’ part of ‘shark’ is barely audible.

Appendix A.

The code used for part 1, 2 and 3.

filter.h

//Author: Sameer Kumar

//Date: 22/11/2006

#include <iostream>

#include <math.h>

#include <fstream>

using namespace std;

class Filter

{

//number of taps

int m_tap;

//Cut-of frequency

int m_fc;

//Sampling frequency

int m_fs;

//Nyquist Frequency

double m_nqFreq;

//Array of coefficients

double *p_mcoeff;

//Take input parameters from prompt

void TakeArguments();

//Calculate coefficients

void CalCoeff();

//Apply windowing function to coefficients

void Windowing();

//Display coefficients

void Display();

//Calculate frequecncy response of FIR filter

void FreqResponse();

public:

//main function to do filtering

void Filtering();

//Destructor to delete memory allocated for coefficeints.

~Filter()

{

delete [] p_mcoeff;

};

};

main.cpp

//Author: Sameer Kumar

//Date: 22/11/2006

// main.cpp : Defines the entry point for the console application.

//

#include “filter.h”

int main(int argc, char* argv[])

{

Filter f;

f.Filtering();

return 0;

}

filter.cpp

//Author: Sameer Kumar

//Date: 22/11/2006

#include “filter.h”

const double PI = 3.14;

//Take input parameters from prompt

void Filter::TakeArguments()

{

//Take taps as a input parameter.

cout<<”Please enter number of taps”<<endl;

cin>>m_tap;

//Take cut-off frequecy as a input parameter

cout<<”Please enter the cut-off frequency”<<endl;

cin>>m_fc;

//Take sampling frequency as a input parameter.

cout<<”Please enter the sampling frequency”<<endl;

cin>>m_fs;

//Calculate nyquist frequency

m_nqFreq = 2 * PI * m_fc/m_fs;

}

//Calculate FIR coefficeints.

void Filter::CalCoeff()

{

//Allocate memory for coefficients

p_mcoeff = new double[m_tap];

int i=1,count;

//To keep our filter causal, shift the centre of it to the middle of taps.

if(m_tap%2 == 0)

{

count = (m_tap – 1)>>1;

}

else

{

count = m_tap >> 1;

}

//Calculate centre value of coefficient.

p_mcoeff[count] = m_nqFreq/PI;

//Calculate remaining values of coefficient.

while(i<=count)

{

p_mcoeff[count-i]=sin(i*m_nqFreq)/PI/i;

p_mcoeff[count+i]=p_mcoeff[count-i];

i++;

}

}

//Apply hanning window function to coefficients.

void Filter::Windowing()

{

for(int i = 0 ; i < m_tap ; i++)

{

double win_val = 0.5 – 0.5 * cos( 2*PI*(i+ 1)/(m_tap+1));

p_mcoeff[i] *= win_val;

}

}

//Display the coefficients of FIR filter.

void Filter::Display()

{

for(int i = 0 ; i < m_tap ; i++)

{

cout<<p_mcoeff[i]<<endl;

}

}

//Calculate frequecy response (Amplitude & Phase)

void Filter::FreqResponse()

{

int degree = 180;

ofstream ampFile,phaseFile;

//Open files to write amplitude and phase values.

ampFile.open (“amp.txt”);

phaseFile.open (“phase.txt”);

double temp;

for(int i=0;i<=degree;i++ )

{

double real=0,imag = 0;

double nqfreq = ((double) 2*i/degree -1)*PI;

//Calculate real and imaginary frequecies.

for(int j=0;j<m_tap;j++ )

{

int time = j – (m_tap – 1)/2;

real += p_mcoeff[j]*cos(time*nqfreq);

imag -= p_mcoeff[j]*sin(time*nqfreq );

}

if(nqfreq > 0)

{

//Calculate amplitude

temp = sqrt(pow(real,2) + pow(imag,2));

ampFile<<temp<<endl;

//Calculate phase.

temp = atan(imag/real);

phaseFile<<temp<<endl;

}

}

//Close txt files.

ampFile.close();

phaseFile.close();

}

//Do filtering

void Filter::Filtering()

{

bool window;

//Take input parameter

TakeArguments();

//Calculate coefficients

CalCoeff();

//Take input parameter for window selection

cout<< “Please Enter the windowing option”<<endl;

cout<<”0: No Window”<<endl;

cout<<”1: With Hanning Window”<<endl;

cout<<”2: With rectangular Window”<<endl;

cin>>window;

//Depending upon window option, select windowing function.

//For rectangular and no window option, w(n) = 1

if(window == 1)

{

Windowing();

}

Display();

FreqResponse();

}

Appendix B.

The code used for Part 4.

main.cpp

//Author: Sameer Kumar

//Date: 22/11/2006

#include “filter.h”

int main(int argc, char* argv[])

{

Filter fp;

char inFile[20],coeffFile[20];

//Take the coefficient file as argument.

cout<<”Please enter the coefficient file”<<endl;

cin>>coeffFile;

//Take the input wave file as argument.

cout<<”Please enter the input wave file”<<endl;

cin>>inFile;

try

{

//Do filtering

fp.Filtering(coeffFile,inFile);

}

catch(Exception e)

{

cout<<e.GetMessage()<<endl;

}

return 0;

}

filter.cpp

//Author: Sameer Kumar

//Date: 22/11/2006

#include “filter.h”

//Check wave header

void WavHeader::ConfirmHeader()

{

//Compares rID with “RIFF”

if(!((rID[0]==’R')||(rID[1]==’I')||(rID[2]==’F')||(rID[3]==’F')))

{

throw Exception(“Wrong RIFF Format”);

}

//Compares wID with “WAVE”

if(!((wID[0]==’W')||(wID[1]==’A')||(wID[2]==’V')||(wID[3]==’E')))

{

throw(“Wrong WAVE Format”);

}

//Compares fID with “fmt”

if(!((fId[0]==’f')||(fId[1]==’m')||(fId[2]==’t')))

{

throw(“Wrong fmt format”);

}

//Checks formatTag

if(wFormatTag!=1)

{

throw(“Wrong Wave Format tag”);

}

//Checks bits per sample with 16 & 8

if( (nBitsPerSample != 16) && (nBitsPerSample != 8) )

{

throw(“Wrong Bits per sample”);

}

}

//Checks whether chunk is data chunk or not

int ChunkHeader::IsDataChunk()

{

//Checks dID with “data”

if((dId[0]==’d') && (dId[1]==’a') && (dId[2]==’t')&& (dId[3]==’a'))

{

return 1;

}

return 0;

}

//Reads coefficients from file.

void Filter::ReadCoeff(char* fileName)

{

ifstream myfile (fileName);

if (myfile.is_open())

{

//Read length of coeffcients

myfile>>m_coeffLen;

//Allocate memory for coeffcients

m_coeff = new double[m_coeffLen];

//Reads all coefficients

for(int i=0;((i<m_coeffLen) && (!myfile.eof()));i++)

{

myfile>> m_coeff[i];

}

myfile.close();

}

else

{

throw(“Can not open coefficient file”);

}

}

//Read Input speech samples

void Filter::ReadSpeechInput(char *fileName)

{

ifstream::pos_type size;

char* memblock;

int index=0;

ifstream file (fileName, ios::in|ios::binary|ios::ate);

WavHeader *hdr = new WavHeader;

ChunkHeader *hdr1 = new ChunkHeader;

int chunkhdrSize = sizeof(ChunkHeader);

ofstream inFile;

inFile.open(“input.txt”);

if (file.is_open())

{

//Get the size of wav file

size = file.tellg();

//Allocate memory for memblock

memblock = new char [size];

file.seekg (0, ios::beg);

//reads all data from file

file.read (memblock, size);

file.close();

}

else

{

throw(“Unable to open file”);

}

//Read wave header

memcpy(hdr,memblock,sizeof(WavHeader));

//Check wave header

hdr->ConfirmHeader();

index = sizeof(WavHeader);

memblock += index;

// read chunks until a ‘data’ chunk is found

int data_chunk = 1;

while(index<(int)size)

{

if(data_chunk > 10)

{

throw(“Too many chunks”);

}

if(((int)size – index) > chunkhdrSize)

{

//read the chunk header

data_chunk++;

index +=chunkhdrSize;

memcpy(hdr1,memblock,chunkhdrSize);

memblock+=chunkhdrSize;

}

else

{

throw(“Wrong data”);

}

//Breaks the loop if , got data chunk

if(hdr1->IsDataChunk()==1)

{

break;

}

}

//Calculate the length of samples

m_inputLen = hdr1->dLen * 8/hdr->nBitsPerSample;

//Allocate memory for input samples

m_input = new double[m_inputLen];

//Reads data in m_input samples

if(hdr->nBitsPerSample == 16)

{

short *temp = (short*)memblock;

for(int i=0;i<m_inputLen;i++)

{

m_input[i] = (double)temp[i];

//Each 64th sample of sound is added to file.

if((i%64)== 0)

{

inFile<<m_input[i]<<endl;

}

}

}

else

{

unsigned char* temp = (unsigned char*) memblock;

for(int i=0;i<m_inputLen;i++)

{

m_input[i] = (double)temp[i];

//Each 64th sample of sound is added to file.

if((i%64)== 0)

{

inFile<<m_input[i]<<endl;

}

}

}

inFile.close();

}

//Convolution of two arrays

void Filter::convol()

{

//Calculate the length of output samples.

int outLen = m_inputLen + m_coeffLen -1;

m_output = new double[outLen];

ofstream outFile;

outFile.open(“output.txt”);

for(int i=0;i<outLen;i++)

{

m_output[i] = 0;

for(int j=0;j<m_coeffLen;j++)

{

if(((i-j)>=0) && ((i-j)<m_inputLen))

{

m_output[i]+=(m_input[i-j]*m_coeff[j]);

}

}

//Each 64th sample of sound is added to file.

if((i%64)== 0)

{

outFile<<m_output[i]<<endl;

}

}

outFile.close();

}

//Filter the input samples

void Filter::Filtering(char* coeffFile,char*wavFile)

{

//Read coefficients

ReadCoeff(coeffFile);

//Read input speech sample

ReadSpeechInput(wavFile);

//Convolution of two arrays

convol();

}

VN:F [1.8.8_1072]
Rating: 0.0/10 (0 votes cast)
VN:F [1.8.8_1072]
Rating: -1 (from 1 vote)

Algorithm Description and Implementation :

The application is built using the C++ Language. The Huffman algorithm is one of the most well known algorithms worldwide and has been used for many years for compressing/encrypting data.

The algorithm’s purpose is to reduce the number of bits representing a byte. Generally, each byte is represented in 8 bits which allows 256 different characters to be represented.

Huffmans algorithm tries to reduce the number of bits required by the most frequently used characters in a certain file or stream. This varies from one file to another.

The implementation of the application has two main programs :


i) encoder.cpp

The program starts by reading the stream byte by byte and calculating the frequency of each byte. The most used characters have a shorter representation.

There may be some characters that the algorithm will give a longer representation to, but these are rare characters.

The corresponding Huffman code is obtained using a priority queue. The STL, standard template libraries for C was referred to in creating this. The file queue contains the priority queue template.

All the bytes are inserted such that the most used comes in the back and the least used goes in the front. The least used 2 bytes are popped from the queue and the sum of the 2 bytes is pushed.

This sum becomes the parent and the 2 bytes that were popped are its children. This process continues until only one node remains. This node becomes head of the Huffman Tree.

As the decoder will decode the file using the same tree then the encoder must save the tree data in the output file for decoding purposes.

This creates a problem regarding the size occupied by the tree data. A well known method for solving this problem is the Canonical Huffman method.

This method involves changing the codes generated for each byte in a systematic way that will give the feature of saving the bytes in the file without saving its corresponding code.

This method will make the header of the output file which holds the tree data to be a maximum of only 265 bytes.

As the original Huffman tree would contain the symbol, number of bits and the code, listing the symbols in sequential order allows the removal of the symbols themselves as they can be determined by the order.

It is also known that later code is higher in value than earlier code. This can be omitted. This leaves the number of bits for each symbol to be transmitted.

For example, the following case would utilise the Canonical Method. A text file that contains only a combination of the character a and b. Let a=97 and b=98.

Using the general Huffman method, it will be stored as the following:

(0,0,0,0,……,length of encoding of a,,0,0,0,0,……,length of encoding code of b)

Here the fields of a and b are filled with many useless zeros. Canonical Huffman stores this as

(length of a,length of b). As can be seen this is saved in only 4 bytes instead of 256 bytes.

The encoder then reads the input stream again byte by byte and inserts the corresponding code in the output file.

Another problem that needs to be solved is the last byte in the stream. As the length of the bytes is changed, the final bit-stream length will not be divisable by 8.

This can be solved by inserting a byte in the header that tells how many bytes will be extracted from the last byte by the decoder.

This saves space in the header and saves the decoder from decoding garbage data.


ii) decoder.cpp

As specified in the encoder section, the decoder will be doing the reverse of the encoder operation.

It starts by building the tree from the saved data. In order to build this tree, the decoder checks the byte value. If the value = 255, the decoder knows is the normal Huffman method.

Else if the byte value of 255 wasnt found, the decoder assumes the canonical Huffman method.

The decoder then reads the encoded file and using the built tree, it decodes it generating the output file.

Each program is a class which its constructor takes two arguments.

The first is the path of the first file (to be compressed file path in case of encoder and previously compressed file in case of decoder) and the second argument takes the path of the file to be stored.

Each of them has only one public method called start ( ) which do all the work. Starting the encoder program will prompt for a file path to be compressed and another one to be saved to.

If no errors happen (the only expected error is failure to open a file), the data will be read from the source file and compressed to the destination file.

The decoder program has the same usage method of the encoder. It will prompt for the path of the previously compressed file and will extract the data to the destination file.

The time used to build Huffman tree for each file and the space occupied (even if it is small) makes it not the best compression method although it works satisfactorily for text files.

The compression ratio for the multimedia file is not significant.

The compression ratios of the different file types is not the same. It depends on the type of file and the number of different characters in the file.

As this number increases, the compressibility may decrease. In fact, many compression programs fail to compress if the size of the file is already extremely small as the output file will even be larger than the original file.

This is due to the tree data and that each character must be represented once. The cost in time is directly proportional to the file size in both compression and decompression operations.

Huffman encoding is practical if the encoded string is large relative to the code table.

Test Procedure :

The following test procedure, for test1, was used to test the application during design. The test procedure was the same for all data. The file names, extensions or locations used are not mandatory.

Compress File:

1. Copy the test file to c:/test1.dat

2. Run the encoder and when prompted for the file to be compressed, enter c:/test1.dat

3. When asked for the file to store, enter c:/test1.huf

4. A file named test1.huf will be created in the c:/ which is the compressed one.

Decompress File :

1. Run the decoder and when prompted for the file to be decompressed, enter c:/test1.huf.

2. When asked for the file to store, enter c:/test2.dat.

3. test1.dat and test2.dat were then compared.

Test Results :

The application was developed in Microsoft Visual Studio 2005 and was tested in DevC++. The application executed successfully ensuring that it does not depend on particular libraries or environments.

It has also been tested on MS Windows XP and MS Windows Vista.

1. test1.dat :

Original Size : 60kb

Compressed Size : 37kb

2. test2.dat :

Original Size : 149kb

Compressed Size : 97kb

3. test3.dat :

Original Size : 23329kb

Compressed Size : 21419kb

Code:

encoder.cpp

//Encoder.cpp
//Sameer Kumar
//27/04/2007
//DME4
//CA438 Multimedia Technology
//Dr. Stephen Blott
#include <iostream>
#include <fstream>
#include <queue>
#include <vector>
#include <stdlib.h>
class encoder
{   //the node for the tree
class TreeNode{
public:
//frequency of the node
unsigned int Frequency;
//the character for this node
unsigned char c;
//length of the encoding code for this node
int length;
//2 pointers for left and right nodes in the tree
TreeNode *child0,*child1;
TreeNode(){Frequency=0;child0=child1=NULL;length=0;c=255;}
//constructor
TreeNode(unsigned char c,unsigned int Frequency)
{
child0=child1=NULL;
this->c=c;
this->Frequency=Frequency;
length=0;
c=255;
}
//constructor
TreeNode( TreeNode* n0, TreeNode* n1)
{
Frequency=n0->Frequency+n1->Frequency;
child0=n0;
child1=n1;
length=0;
}
//this operator is used for priority queue to determine the location in which to put the inserted node
bool operator <(const TreeNode& a)const{return a.Frequency
};
//the canonical sort used by the sort algorithm
struct CanonicalHuffman{
bool operator()( TreeNode n1,TreeNode n2 ){
if(n1.length
return 1;
else if(n1.length==n2.length)
if(n1.c
return 1;
return 0;
}
};
//used by sort algorithm to sort nodes alphabetically
struct AlphaSort{
bool operator()( TreeNode n1,TreeNode n2 ){
return(n1.c
}
};
std::ifstream in;
std::ofstream out;
//frequency of the 266 bytes
TreeNode SymbolsFrequency[256];
//number of symbols that exist in the file
int ActiveSymbols;
//the number of bytes that will be in the last byte
char LastByte;
//for bulding huffman tree
std::priority_queue < TreeNode,std::vector
> q;
//traverse huffman tree to get the encode for each byte
void traverse(const TreeNode *Tree,int length){
//go through the whole tree to get the encoding
if(!Tree->child0&&!Tree->child1)
{
SymbolsFrequency[Tree->c].length=length;
return;
}
length++;
if(Tree->child0)
traverse(Tree->child0,length);
if(Tree->child1)
traverse(Tree->child1,length);
}
//calculates the frequency of all the bytes storing them in the SymbolsFrequency array
void calculateFrequency()
{   //calculate all frequencies
while(true)
{
//calculate the frequency of the bytes of the input stream
unsigned char c;
c=in.get();
if(in.eof())break;
if(!SymbolsFrequency[c].Frequency)
{ActiveSymbols++;SymbolsFrequency[c].c=c;}
SymbolsFrequency[c].Frequency++;
}
//building the tree using priority queue
//push all the non-zero frequency nodes to the tree
for(int i=0;i<256;i++){ 			SymbolsFrequency[i].c=i; 			if(SymbolsFrequency[i].Frequency) 				q.push(SymbolsFrequency[i]); 			} 		while(q.size()>1)
{
TreeNode *child0=new TreeNode(q.top());
q.pop();
TreeNode *child1=new TreeNode(q.top());
q.pop();
q.push(TreeNode(child0,child1));
}
//traverse the tree to obtain pre-final codes
traverse(&q.top(),0);
//start canonical method and sort accordingly - changing codes
std::sort(SymbolsFrequency,SymbolsFrequency+256,CanonicalHuffman());
int code=-1;
int ii=0;
while(!SymbolsFrequency[ii].Frequency)
ii++;
int current_length=SymbolsFrequency[ii].length;
//canonical huffman algorithm
for(;ii<256;ii++)
{
code++;
if(SymbolsFrequency[ii].length!=current_length)
{
current_length=SymbolsFrequency[ii].length;
code=code<<1; 			} 			SymbolsFrequency[ii].Frequency=code; 		} 		//end of canonical method 		//re sort alphabetically 		std::sort(SymbolsFrequency,SymbolsFrequency+256,AlphaSort()); 		//check whether it is better to store associative huffman table or canonical table 		//for the associative table...127byte*2 bytes = 254 bytes...use this if it is greater than 127          if(ActiveSymbols>127)
for(int i=0;i<256;i++)
out.put((unsigned char)SymbolsFrequency[i].length);
//use the canonical table because the ActiveSymbols are less than or equal to 127
else
{   //canonical table indication mark
out.put((unsigned char)255);
out.put((unsigned char)ActiveSymbols);
for(int i=0;i<256;i++)
if(SymbolsFrequency[i].length)
{
     out.put(SymbolsFrequency[i].c);
     out.put((unsigned char)SymbolsFrequency[i].length);
}
}
//restore the input stream to the beginning
in.clear();
in.seekg(0,std::ios::beg);
}
public: 	encoder(char *in,char *out)
{
this->in.open(in,std::ios::binary);
this->out.open(out,std::ios::binary);
ActiveSymbols=LastByte=0;
}
~encoder()
{
in.close();
out.close();
}
void start()
{
if(in.is_open()&&out.is_open())
{
//reserve place for last byte length to be stored
out.put(LastByte);
std::cout<<"Building Huffman Tree..."<
//calculate the frequency and build the tree and save it to the output stream
calculateFrequency();
std::cout<<"Huffman Tree Building finished."<
unsigned char b=0,c,inserter=1<<7;
unsigned int checker;
std::cout<<"Compressing file..."<
while(true)
{
c=in.get();
if(in.eof())break;
checker=1<<(SymbolsFrequency[c].length-1); 				while(checker) 				{ 					if(checker&SymbolsFrequency[c].Frequency)                         //packing bytes in one byte 						b=b|inserter; 					inserter=inserter>>1;
//wait till full 8-bits are packed
if(!inserter)
{
inserter=1<<7; 						out.put(b); 						b=LastByte=0; 					} 					checker=checker>>1;
}
LastByte++;
}
//there may be a bit in the last byte that contains data and garbage. check if any data is present to save it
if(inserter)
out.put(b);
//setting last byte length
//modify the size of the last byte length as it may be invalid
out.seekp(std::ios::beg);
out.put(LastByte);
out.flush();
std::cout<<"Operation Completed."<
}
else
std::cout <<"Error Opening Files."<
}
};
int main(){
char FP[256],HP[256];
//prompt user
std::cout<<"Enter Source Path To Be Compressed: ";
std::cin.getline(FP,256);
std::cout<<"Enter Destination Path To Store: ";
std::cin.getline(HP,256);
encoder e(FP,HP);
e.start();
}

decoder.cpp

// Decoder.cpp
//Sameer Kumar
//27/04/2007
//DME4
//CA438 Multimedia Technology
//Dr. Stephen Blott
#include <iostream>
#include <fstream>
#include <queue>
#include <vector>
#include <stdlib.h>
class decoder{
//the node for the tree
class TreeNode{
public:
//frequency of the node
unsigned int Frequency;
//the character for this node
unsigned char c;
//length of the encoding code for this node
int length;
//2 pointers for left and right nodes in the tree
TreeNode *child0,*child1;
TreeNode(){Frequency=0;child0=child1=NULL;length=0;c=255;}
//constructor
TreeNode(unsigned char c,unsigned int Frequency)
{
child0=child1=NULL;
this->c=c;
this->Frequency=Frequency;
length=0;
c=255;
}
//constructor
TreeNode( TreeNode* n0, TreeNode* n1)
{
Frequency=n0->Frequency>>1;
length=n0->length-1;
if(n0->Frequency&1)
{child1=n0;child0=n1;}
else
{child1=n1;child0=n0;}
}
//this operator is used for priority queue to determine the location in which to find the inserted node
bool operator <(const TreeNode& a)const
{
if (length==a.length)
return Frequency<a.Frequency;
return length<a.length;
}
};
////the canonical sort used by the sort algorithm
struct CanonicalHuffman{
bool operator()( TreeNode n1,TreeNode n2 ){
if(n1.length<n2.length)
return 1;
else if(n1.length==n2.length)
if(n1.c<n2.c)
return 1;
return 0;
}
};
//used by sort algorithm to sort nodes alphabetically
struct AlphaSort{
bool operator()( TreeNode n1,TreeNode n2 ){
return(n1.c<n2.c);
}
//end of TreeNode class
};
//input stream
std::ifstream in;
//output stream
std::ofstream out;
//array for all the possible 256 characters. it will be converted to the tree from here using the queue
TreeNode SymbolsFrequency[256];
//the head of the tree
TreeNode* Tree;
//the number of bytes that will be extracted from the last byte
unsigned char LastByte;
//the queue
std::priority_queue < TreeNode,std::vector<TreeNode> > q;
//builds the huffman tree from the stored table determining whether it is from the general huffman table or canonical table
void build_tree()
{
unsigned char b;
b=in.get();
//canonical mark indicated
if(b==255)
{
unsigned char i=in.get();
while(i)
{
unsigned char symbol,len;
symbol=in.get();
len=in.get();
SymbolsFrequency[symbol].length=len;
SymbolsFrequency[symbol].c=symbol;
i--;
}
}
//huffman method
else
{
SymbolsFrequency[0].length=b;
SymbolsFrequency[0].c=0;
int i=1;
do
{
SymbolsFrequency[i].length=in.get();
SymbolsFrequency[i].c=i;
}while(++i<256);
}
//start of canonical method by sorting
std::sort(SymbolsFrequency,SymbolsFrequency+256,CanonicalHuffman());
int ii=0;
int code=-1;
while(!SymbolsFrequency[ii].length)
ii++;
int current_length=SymbolsFrequency[ii].length;
//generate huffman codes
for(int i=ii;i<256;i++)
{
code++;
if(SymbolsFrequency[i].length!=current_length)
{
current_length=SymbolsFrequency[i].length;
code=code<<1;
}
SymbolsFrequency[i].Frequency=code;
}
//canonical codes generated
//build the tree
TreeNode* t=Tree;
//start by filling the queue
for(;ii<256;ii++)
q.push(SymbolsFrequency[ii]);
//loop till only one node remains in the queue.
while(q.size()>1)
{
TreeNode *child0=new TreeNode(q.top());
q.pop();
TreeNode *child1=new TreeNode(q.top());
q.pop();
q.push(TreeNode(child0,child1));
}
//the remaining node is the head of the tree. tree has been constructed successfully
Tree=(TreeNode*)&q.top();
}
public:
//constructor
decoder(char *in,char *out)
{
this->in.open(in,std::ios::binary);
this->out.open(out,std::ios::binary);
Tree=new TreeNode;
}
//start decode
void start(){
if(in.is_open()&&out.is_open())
{
//catching last byte length
LastByte=in.get();
std::cout<<"Building Huffman Tree..."<<std::endl;
//call for building the tree first
build_tree();
std::cout<<"Huffman Tree Building finished."<<std::endl;
TreeNode* t=Tree;
unsigned char b,bb,checker=1<<7;
b=in.get();
bool isLastByte=false;
std::cout<<"Decompressing file..."<<std::endl;
//decompressing data
while(!in.eof()){
bb=in.get();
//indicates if it is decoding the last byte of the code. stops the decoder from decoding any garbage data
if(in.eof())
isLastByte=true;
for(int i=0;i<8;i++)
{
//read each byte and dig down through the tree till we reach the required byte
if(b&checker)
t=t->child1;
else t=t->child0;
//byte found
if(!(t->child0||t->child1)){
//put the byte to the output stream
out.put(t->c);
if(isLastByte){
LastByte--;
if(!LastByte)
break;
}
//reset the tree digger pointer to the top again to make it ready for decoding new bytes
t=Tree;
}
checker=checker>>1;
}
checker=1<<7;
b=bb;
}
//flushing output
out.flush();
std::cout<<"Operation Completed."<<std::endl;
}
else
std::cout <<"error opening files."<<std::endl;
}
};
int main(){
char FP[256],HP[256];
//prompt user
std::cout<<"Enter source path to be decompressed: ";
std::cin.getline(FP,256);
std::cout<<"Enter destination path to extract: ";
std::cin.getline(HP,256);
decoder d(FP,HP);
 d.start();
}

References:

01. Algorithms in C, Robert Sedgewick, Addison-Wesley Publishing, ISBN 0-201-51425-751425

02. http://en.wikipedia.org/wiki/Priority_queue

03. http://datacompression.info/Tutorials.shtml

04. http://www.dogma.net/markn/articles/pq_stl/priority.htm

05. http://en.literateprograms.org/Huffman_coding_(C_Plus_Plus)

06. http://en.wikipedia.org/wiki/Huffman_coding

07. http://www.compressconsult.com/huffman/

08. http://en.wikipedia.org/wiki/Canonical_Huffman_coding

09. http://www.daniweb.com/code/snippet5.html

10. http://planetmath.org/encyclopedia/HuffmansAlgorithm.html

11. http://www.arturocampos.com/cp_ch3-1.html

12. http://www.cprogramming.com/tutorial.html

13. http://cprogramming.com/tutorial/computersciencetheory/huffman.html

14. http://www.cprogramming.com/tutorial/lesson15.html

15. http://mitpress.mit.edu/sicp/full-text/sicp/book/node41.html

16. http://www.inference.phy.cam.ac.uk/mackay/itprnn/code/c/compress/

17. http://www.codeproject.com/cpp/HandyHuffmanCoding.asp

18. http://www.answers.com/topic/huffman-coding

19. http://michael.dipperstein.com/huffman/index.html

20. http://www.cs.duke.edu/csed/poop/huff/info/

21. http://www.datastructures.info/huffman-encoding-algorithm/

22. http://www.codeproject.com/cpp/Huffman_coding.asp

The files can be downloaded here.

VN:F [1.8.8_1072]
Rating: 0.0/10 (0 votes cast)
VN:F [1.8.8_1072]
Rating: 0 (from 0 votes)

Chapter 1 – Introduction

1.1 High Concept

The aim of this project is to incorporate fluid dynamics and particles systems to model natural phenomena in real time. Using the Navier-Stokes equations as a basis for modelling fluid flow, the resultant output will be feed into the particle system. All of the calculations will be performed using the GPU, displaying the robustness and benefits of designing such an application.

1.2 Project Overview

The animation of gases and fluids such as smoke, water and fire provide some of the most stunning visuals in computer games today. These can be achieved with a variety of methods and currently many projects are looking to expand on previous techniques. Particle systems are used in all modern action / shooter games to model a variety of phenomena such as explosions, smoke and clouds.
Volcanic activity is chosen as it represents a combination of gas, fluid and combustion. When these phenomena interact, the results should model the real world. Common particle systems use simple static forces to update individual pieces. Modelling these interactions is CPU intensive due to the sheer number of calculations involved. When utilising particles for achieving a realistic visualisation of any phenomena, the main factor that determines the quality is the number of particles used to represent the system. Increasing this number automatically increases the overhead on the system. Data transfer between the CPU and the GPU is limited and so a bottleneck will occur if systems containing tens of thousands of particles are rendered. Performing all of these calculations on the GPU erases this issue. It can also update particle positions freeing up more CPU power. The GPU can also calculate other environmental forces and apply those to the particle system as will be demonstrated by this project.

This project will focus on

  • Techniques for modelling interactions between natural phenomena
  • Performing these interactions on the GPU in real-time

1.3 Project Scope

This paper aims to investigate research and current techniques across the fields of fluid dynamics, particle systems and GPU programming. The following chapter summarises these technologies.
After considering the literature, a suitable method for designing the application is derived. The chosen method will be compatible for designing other types of similar phenomena.

Chapter 2 – Background and Literature Review

This chapter deals with the methodologies required to complete this project.

2.1 Fluid Dynamics

2.1.1 Fluid Physics

The field of fluid dynamics deals with the study of the motion of fluids and gases. The Navier-Stokes equations can be used to describe this motion. A fluid with constant density and temperature is described by a velocity and pressure field. The equations model the relationship between the velocity and the rates of change of the pressure in a given fluid.

THE NAVIER STOKES EQUATIONSFIGURE 1 – THE NAVIER STOKES EQUATIONS [ STAM03 ]
Equation 2.1 shows the velocity in a compact vector notation and 2.2 shows a density moving through the velocity field.

From this, the main properties are:

Advection:

The velocity of a fluid causes the fluid to transport objects, densities and other quantities along with the flow.

Pressure:

When force is applied to a fluid, it does not instantly affect the entire volume. Instead the pressure builds up, and any pressure leads to an increase in acceleration.

Diffusion:

Thick fluids such as lava have a high viscosity. Viscosity is a measure of how resistive a fluid is to flow. This resistance results in diffusion of the momentum.

External Forces:

These forces may be either local forces or body forces. Local forces may be applied to a specific region of the fluid whilst body forces are applied equally to the entire fluid.

2.1.2 Fluid Dynamics in Computer Games

There are different methods available to model fluid dynamics. The primary concern is the final aesthetic effect and the ability to run in real time so accuracy is off a lesser value. Modelling fluids with accuracy in mind would result in a computationally expensive system.

To use the Navier-Stokes equations, the space must be converted into a grid. The density of the object, velocity and pressure at the centre of each grid cell must be defined. At each time step, the equations are used to determine the new values of the fields. An extra layer of cells is provided to account for the boundary conditions.

GRID WITH BOUNDARY CELLSFIGURE 2 – GRID WITH BOUNDARY CELLS [ STAM 03 ]

Many systems have been created which can model accurate and realistic simulations. However, these are not designed to run in real-time. [ Fedkiw01 ] assumes that the fluid flow is of constant density and uses Euler equations. It uses an implicit semi-Lagrangian scheme which ensures stability over the time step. This algorithm builds a new grid of velocities from the ones already computed by tracing the midpoints of each voxel face through the velocity field. New velocities are then interpolated at these points and the values are transferred to the cell they originated from.

[ Stam03 ] describes a system which utilises an Eulerian grid approach. This uses a coarser resolution than that required for accurate simulations reducing the overhead. The simulation space is divided in a 3D voxel grid. The Navier-Stokes equations for the gas volume in each cell are calculated by applying numerical methods. It provides an implicit integration scheme, allowing the system to run at any time step. This is based on equation 2.2. The solver is used to evolve the smoke density over time in each cell given the rate of diffusion and a static velocity field. This moves the smoke and allows changing forces caused by the environment.

RESULTS FROM SOLVER INPLEMENTED IN [ STAM03 ]FIGURE 3 – RESULTS FROM SOLVER INPLEMENTED IN [ STAM03 ]

2.1.3 Rendering Fluid

Rendering fluids or gases is known as volume rendering. Mentioned in the previous section, the space must be divided into a 3D grid containing the gas densities of each cell. [ Stam03 ] uses a particle based approach which moves the particles through the generated velocity field instead of the gas densities. The densities can be modelled as particles at the centre of each cell.

2.2 Particle Systems

2.2.1 General Particle Systems

A particle is a very small object that is usually modelled as a point mathematically. Point sprites can have textures mapped to them and the sizes can vary. The conventional particle system contains an emitter which creates the particles. This also controls their life-span and applies initial forces to them. When the particle is instantiated, it moves to a random 3D point
within a specified area and constantly moves under the force applied to it.
Traditionally the particles were stored in a vertex buffer in the main system memory and these points were sent to the GPU. This method can handle up to 10000 particles before bottlenecks begin to occur. This is computationally expensive for large particle systems and so the current trend is to compute all the required information on the GPU.

2.2.2 GPU – Based Particle Systems

Parallel processing capabilities of computers today allow systems to perform calculations simultaneously. All updating and rendering is performed on the GPU. The programmable graphics pipeline and recent hardware advances affords this flexibility. The previous technique of streaming the vertices of the particles is eliminated freeing the CPU for other computational work.

[ Latta04 ] utilises the above implementation creating up to a million particles on the GPU. The velocities and positions of all particles are stored in textures. [ Kolb04 ] describes a system for collision detection using depth maps that represent the outer shape of an object. These store distance values and normal vectors. Additional textures can store other collision data allowing the pixel shader to bounce or merge objects. [ Kruger05 ] implements system which utilises fluid dynamics techniques. For 2D, it uses pre-computed effect templates to update the particle positions which are stored in textures. By using volumetric extrusion, the techniques are extended to work in 3D, but are only suitable for symmetric effects.

2.2.3 Verlet Integration

In order to update the positions of particles over time, Verlet integration will be used. This states that the velocity of the particle is deduced by comparing the previous position to the one before. [ Latta04 ]

x ‘= 2x – x* + a . Δt^2
x* = x
x* = previous position

Chapter 3 – Method

3.1 Overview

The following section details the intended method to produce the final result. The high level outcomes of this project are

  1. A system which calculates dynamic forces and moves the particles according to these.
  2. Implement a GPU based system on the above factors which controls the entire scene.

3.2 CPU System

3.2.1 Particle Generation

Particles are required for four types of phenomena:

  • Fire
  • Smoke
  • Explosions
  • Lava Flow

The emitter system which creates the particles will store them as texels in two procedural textures. The textures will store their previous and current positions. Another texture stores other static data for each particle as shown in Figure 4. The rate of particles emitted will depend on the type of phenomena. Created particles are assigned a random point and velocity. As the Verlet integration generates a previous and current position, the origin, of the particle, the GPU can manipulate it using this data. Once a particle dies, its texel must be freed for the use of new particles.

FIGURE 4 – PARTICLE DATA STORAGE IN MULTIPLE TEXTURES [ LATTA04 ]

3.2.2 Phenomena Velocities

Each of the phenomena shall have different forces acting on them. Lava will contain higher viscosity and advection values. Viscosity will increase over time changing the texture. Each one will be contained in a 3D grid. Each cell of the grid will have a velocity assigned representing the air flow of that particular cell. The solvers utilised in [ Stam03 ] will form a basis to calculate the velocities which will then be stored in a procedural texture for GPU access. An extra layer of cells will be created to represent the boundary conditions. This will allow realistic interactions between the different phenomena. The procedural texture generated will represent the entire velocity grid.

3.3 GPU System

The shader will determine the new positions of the particles depending on the velocity of the cell in question, current position and the previous position of the particle itself. The outer cells will be affected by the cells from other grids. The new position is stored in a particle position buffer and the vertex shader positions this in the game world. After particles have been translated, they have to be depth sorted into a suitable render order. They will also undergo alpha blending as all phenomena will be inter connected. Point sprites are used to render the particle which is rotated according to the stored value.

3.4 Alternative Method

Perlin Noise

The entire system could be designed with the use of Perlin noise. [ Perlin01 ] utilises a method to create lava flow and fire with smoke. Perlin Noise is extended to incorporate advection. [Tatarinov07 ] creates fire using Perlin Noise to perturb a basic fire shape. This could be extended to the other intended features.

ANIMATED LAVA FLOW

FIGURE 5 – ANIMATED LAVA FLOW [ PERLIN01 ]

BASIC SHAPE (LEFT) USED TO CREATE FIRE EFFECT (RIGHT) [ TATARINOV07 ]

FIGURE 6 – BASIC SHAPE (LEFT) USED TO CREATE FIRE EFFECT (RIGHT) [ TATARINOV07 ]

3.5 Evaluation Criteria

The final application will be accessed according to two factors:

  1. Visual Effect – Qualitative
  2. Application Performance – Quantitative

3.5.1 Qualitative

The main aim of this project is to create an aesthetically pleasing and display the realistic interactions between the different types of phenomena which can then be transferred to a game situation. Comparisons shall be made against projects using similar technology.

3.5.2 Quantitative

The application must run at a minimum of 30 fps. There may be trade-off with the qualitative criteria to achieve this. Code optimisations may have to be performed. This would entail modifications to the algorithms or to the data structures. Tests against running the system only on the CPU will also be performed.

Chapter 4 – Project Plan

Figure 7 gives an overview of the monthly project deliverables. Appendix I contains a detailed Gantt chart.

MONTHLY PROJECT DELIVERABLES OVERVIEWMONTHLY PROJECT DELIVERABLES OVERVIEW

FIGURE 7 – MONTHLY PROJECT DELIVERABLES OVERVIEW

Chapter 5

5.1 References

Latta, L. 2004. Building a Million Particle System, http://www.gdconf.com/conference/archives/2004/latta_lutz_01.pdf (Date Accessed: 08/05/08)
Kolb, A., Latta, L., Rezk-Salama, C. 2004. Hardware-Based Simulation and Collision Detection for Large Particle Systems, http://www.2ld.de/gh2004/HWSimCollisionPS.pdf (Date Accessed: 08/05/08)
Kruger, J., Westermann, R. 2005. GPU Simulation and Rendering of Volumetric Effects for Computer Games and Virtual Environments. http://wwwcg.in.tum.de/Research/data/Publications/eg05.pdf (Date Accessed: 08/05/08)
Stam, J. 2003. Real-Time Fluid Dynamics for Games. http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf (Date Accessed: 08/05/08)
Petersson, S. 2007. Particle System Rendering – The effect on rendering speed when using Geometry Shaders. http://www.bth.se/fou/cuppsats.nsf/all/10165efe6d0ac508c12573470032c16a/$file/070518_DVC001_Stefan_Petersson_850417_Rapport.pdf (Date Accessed: 08/05/08)
Anderson, L. 2005. Real-Time Fluid Dynamics for Virtual Surgery. http://www.ce.chalmers.se/~uffe/xjobb/blood.pdf (Date Accessed: 08/05/08)
Fedkiw, R., Stam, J., Jensen, H.W. 2001. Visual Simulation of Smoke. http://physbam.stanford.edu/~fedkiw/papers/stanford2001-01.pdf (Date Accessed: 08/05/08)
Stora, D., Agliati, P.-O., Cani, M.-P., Neyret, F., Gascuel, J.-D. 1999. Animating Lava Flows. http://www-evasion.imag.fr/Publications/1999/SACNG99/gi99.pdf (Date Accessed: 08/05/08)
Tatarinov, A. 2007. Perlin Fire. http://developer.download.nvidia.com/whitepapers/2007/SDK10/PerlinFire.pdf (Date Accessed: 09/05/08)
Perlin, K. 2001. Flow Noise. ., http://chihara.aist-nara.ac.jp/people/2003/takasi-a/research/short_paper.pdf (Date Accessed: 09/05/08)
Amada, T., Imura, M., Yasumuro, Y., Manabe, Y., Chihara, K., Particle-Based Fluid Simulation on GPU, http://chihara.aist-nara.ac.jp/people/2003/takasi-a/research/short_paper.pdf (Date Accessed: 09/05/08)
Mizuno, R., Dobashi, Y., Nishita, T. 2002. Volcanic Smoke Animation using CML, http://nis-lab.is.s.u-tokyo.ac.jp/nis/cdrom/cgi/ics2002-crf5.pdf (Date Accessed: 09/05/08)
Miyazaki, R., Yosida, S., Dobashi, Y., Nishita, T. 2001. A Method for Modeling Clouds based on Atmospheric Fluid Dynamics, http://nis-ei.eng.hokudai.ac.jp/~doba/papers/ PG2001_CML.pdf (Date Accessed: 09/05/08)
Nvidia Developer, http://developer.nvidia.com/page/home.html, (Date Accessed: 01/05/08)
ATI Developer, http://ati.amd.com/developer/, (Date Accessed: 01/05/08)

5.2 Bibliography

Cartwright, C. 2008, SA1100A Masters Proposal Course Notes, University of Abertay, Dundee, https://portal.abertay.ac.uk/learning/module/co/sa/SA1100A07/
Fernando, R. 2004, GPU Gems, Nvidia Developer Zone, http://developer.nvidia.com/object/gpu_gems_home.html
Dickheiser, M. 2006, Game Programming Gems 6, Charles River Media, Boston.
Nyguyen, H. 2007, GPU Gems 3, Addison Wesley Professional, Boston.
Luna, F. D. 2006, Introduction to 3D Game Programming with DirectX 9.0c : A Shader Approach, Wordware Publishing, Texas.

The .pdf version of this report can be viewed here.

VN:F [1.8.8_1072]
Rating: 0.0/10 (0 votes cast)
VN:F [1.8.8_1072]
Rating: 0 (from 0 votes)