Digital Signal Processing System Design

Homework 1

Due: March 16, before the class  (선배의 숙제를 보는 것 금지. 어길 경우 F).

1. We need to interpolate an 8KHz sampled input speech signal into 32KHz one.   The conceptual block diagram is shown below, but this does not mean that you have to implement each box separately. You may merge the upsampling and FIR filtering processes.  You can use Matlab or any other tools for filter design.  Use floating-point arithmetic.

a) As for the lowpass FIR filter for interpolation, the transition bandwidth (wstop - wpass) would be less than 10% of the wstop.  The passband ripple allowed is 0.5 dB, and the stopband ripple is -50dB.   Please design the needed FIR filter using Matlab.  You first need to determine the approximate filter order.

b) Develop a C program for interpolation using the above FIR filter.   Employ the computation reduction method for interpolation.

c) Use an IIR filter for the interpolation instead of FIR filter.  Design the IIR filter.

d) Develop a C program for interpolation using the IIR filter and see if there is any advantage in terms of the execution speed.  It would be best if you can measure and compare the CPU time when you execute your programs.   You may use Windows or Linux environment.

2. Prove the following symmetry relation of the discrete time Fourier transform of a complex sequence.  Consult any DSP book.

Let F{x[n]} = X(ejw)

a) x[-n] => X(e-jw) and x*[-n] => X*(ejw)

b) Re{x[n]} => Xcs(ejw) = 1/2{ X(ejw) + X*(e-jw) }

and jIm{x[n]} => Xca(ejw) = 1/2{ X(ejw) - X*(e-jw) }

c) xcs[n] => Xre(ejw) and xca[n] => j Xim(ejw)

Note: 'cs' means complex-symmetric and 'ca' means complex-antisymmetric

3. (a) Using the above property, develop a procedure that computes the FFT of the real 1K-point data (the input data x[n] is real) using 512-point (complex) FFT routine.  Verify this by developing a C-program.  Here, you need to compare the results of the straightforward 1K-point FFT and your program using 512-point FFT and some pre- and post-processing.

(b) You have two real-sequences, x[n] and y[n], both have the length of 1K.  Using 1K-point FFT, computes X[k] and Y[k].  Verify this by developing a C-program.

4. You have a sequence of length 1K whose odd numbered samples are all zeros (x[1] = x[3] = x[5] = .. = x[1023] = 0).   Develop a procedure how you can evaluate X[k] just using 512-point FFT.  Verify this by developing a C-program.

5. Reading materials.

1) [Reading this one is required] RW Schafer, LR Rabiner, "A digital signal processing approach to interpolation, " Proceedings of the IEEE, 1973 - ieeexplore.ieee.org

2) [Recommend, this material is easy] Visit the web site explaining interpolation and decimation.

http://www.dspguru.com/info/faqs/mrfaq.htm

3) [Recommend, good for application to communication system design]  S. Jou, S. WU, and C. Wang, "Low-power multirate Architecture for IF Digital Frequency Down Converter," IEEE Trans. on Circuits and Systems-II, Nov 1998.

Digital Signal Processing System Design

Homework 3

Due: April 6, before the class  (선배의 숙제를 보는 것 금지. 어길 경우 F).

다음의 숙제를 위해서 Intel C++ compiler를 사용하시오.

1. Wikipedia 등을 통해서 현재 본인이 쓰고 있는 컴퓨터에서 다음의 사양을 조사하시오.

# of CPU core, L1 cache size, L2 cache size, integer instruction의 파이프라인 깊이, floating-point instruction의 파이프라인 깊이, 한번에 issue할 수 있는 명령어의 갯수.

2. Matrix vector 연산 C = A*B 를 실행하고자 한다.  C=1,024x1 vector, A=1,024x1,024 matrix, B= 1,024x1 vector이다.

a) 본 matrix 계산을 위해서 필요한 연산 (곱셈과 가산)의 갯수는 총 몇 번인가?

b) 여기서 C[0], C[1], ... C[1024]를 순차적으로 구하는 프로그램을 작성하고, 이 프로그램을 1,000번 돌릴 때 걸리는 시간을 측정하시오.   또 이 프로그램을 실행시에 총 필요한 메모리 load, store의 횟수는 얼마인지 대략 계산을 하시오.

for (j=0; j<1,024; j++){

c[j] = 0.0;

for (i=0; i<1,024; i++)

c[j]+=a[j,i]*b[i];

}

c) 위의 프로그램을 아래와 같이 변형한다.  위의 경우와 마찬가지로 프로그램을 1,000번 돌릴 때 걸리는 시간과 메모리 load, store의 회수를 대략 계산하시오.

for (j=0; j<1,024; j++)

c[j] = 0.0;

for (j=0; j<1,024; j=j+4){

for (i=0; i<1,024; i++){

c[j]+=a[j,i]*b[i];

c[j+1]+=a[j+1,i]*b[i];

c[j+2]+=a[j+2,i]*b[i];

c[j+3]+=a[j+3,i]*b[i];

}

}

d) 위의 프로그램을 변형해서 소위 outer product방식으로 바꾸고, 앞서와 마찬가지로 실행시간과 메모리 load, store의 회수를 대략 계산하시오.  참고할 C program은 아래와 같다.

for (j=0; j<1,024; j++)

c[j] = 0.0;

for (i=0; i<1,024; i++) {

for (j=0; j<1,024; j++)

c[j]+=a[j,i]*b[i];

}

e) 위의 측정결과를 분석하여서 Pentium의 어떤 특성 때문에 프로그램 실행 시간의 차이가 생겼는지 기술하시오.

3. 위의 matrix, vector multiplication에서 각각의 크기가 c=1,024x1 vector, A=1,024x100,000 B=100,000x1 으로 바뀌었다.  가장 시간이 적게 걸리는 프로그램을 짜시오. 단 이 때 이 프로그램 실행 시에 L1 cache miss를 어떻게 작게 할지를 고려하여야 한다.

4. 아래는 아주 간단한 1st order recursive equation이다. 이식을 C 프로그램으로 작성해서 총 1M input sample을 1000번 돌렸을 때 걸리는 시간을 재시오.  a는 -1과 +1 사이의 상수이다.   아래에 이해를 돕기 위해서 C 프로그램을 준다 (사소한 오류가 있을 수 있음.)  데이타 타입은 부동소수점이다.

y[n] = a*y[n-1] + x[n]

// 아래를 천번 돌려서 시간 측정

y[0]=0.0;

for (i=1; i=<1,000,000; i++)

y[i] = acoef*y[i-1];

5. SSE3 명령어 세트를 이용하여 다음의 식을 프로그램하시오.  이 때 intrinsic function을 사용하면 된다.  하나의 SSE3 명령어로 최대 4개의 부동소수점 연산을 할 수 있다.   단 SSE3를 위한 intrinsic function을 공부할 여유가 없는 학생은 보통의 C program으로 작성해서 시간을 측정하시오.  이 경우 절반의 감점이 있다 (학점 맞는데 큰 문제없는 감점.)

위의 식을 이용하여 역시 총 1M개의 output sample (y[1]~y[1,000,000])을 구하는 프로그램을 1000번 돌리는 시간을 구해서 4번 문제의 경우와 비교를 하기 바란다.  모든 데이타 연산은 부동소수점이며, a, a2, a3, a4은 미리 계산한 상수 값으로서 -1과 +1 사이의 값이다.  데이타 타입은 부동소수점이다.

참고로 위식의 계산 결과는 4의 결과와 같아야 한다.

이 식은 다음과 같이 유도가 된 것이다.

y[2] = a*y[1]+x[2]

-> y[2]=a*{a*y[0]+x[1]}+x[2]

-> y[2]=a*{a*(a*y[-1]+x[0])+x[1]}+x[2] = a3 y[-1] + a2 x[0] + a x[1] + x[2]

6. 4와 5의 결과를 의논하시오.  5의 경우 연산회수가 훨씬 더 많아졌다.  실행 시간도 길어졌는지?

참고:  본 실험은 Intel C++ compiler를 이용하여 실행시키고, 자세한 성능 측정은 아래의 방법으로 할 수 있다.  참고로 Intel C++ compiler는 한달짜리를 무료로 다운로드 받아서 본인의 PC에서 실행을 시키거나, 연구실의 Linux server를 이용할 수 있다.

또 프로파일관련은 아래의 조교의 메일을 참고하기 바란다. 추가사항은 유기선 박사과정에게 문의 ksyou@dsp.snu.ac.kr

단순히 수행시간만을 측정한다면 rdtsc로 충분할것으로 생각하고, 다른 사항들까지 측정해야한다면 VTune이나 OProfile이 적합할것으로 생각됩니다.

VTune은 GUI라서 사용이 편리한대신, 학생들이 사용하려면 nomachine이나 vnc같은 클라이언트 툴을 이용 GUI환경으로 접속하여 사용하여야합니다.

OProfile은 GUI가 필요 없지만, 사용법이 복잡하고 수퍼유저 권한이 필요하기때문에 연구실 서버를 이용하게 하기에는 부적합할 것 같습니다.

1. Intel VTune Performance Analyzer.

CPU안의 Performance Counter를 사용하여 상대적으로 정확한 성능 결과를 얻을 수 있다.

(각 event마다 일정 수준의 Count Unit을 정해놓고 그 이상이 되는 값에 대해서 계측을 하는 방식이므로 Cycle-accurate Simulator처럼 매 사이클의 정보까지는 얻어낼수는 없다)

Clock cycles, Cache misses, Bus utilization 등등 많은 종류의 event에 대하여 성능 측정이 가능하다.  GUI환경에서 동작하여 사용이 편리하다.

단점은 Event의 종류가 매우 많아서 실험에 적합한 event를 알아내기가 쉽지 않다.

(단순히 수행 시간, cache miss등을 볼 때에는 straightforward하므로 문제가 되진 않는다)

2. OProfile

VTune과 마찬가지로 Performance Counter를 사용한다.

VTune과의 차이점은 GUI환경이 아닌, 커맨드라인 명령어를 사용하여야 하므로 사용법 숙지가 어렵다. 그리고  측정을 위해서는 수퍼유저 권한이 필요하다.

단점은 VTune과 마찬가지

3. rdtsc 명령어

x86 instruction set중에 현재 clock tick값을 리턴해주는 instruction인 rdtsc를 사용한다.

clock tick값의 차이를 이용하여 측정하고자 하는 부분의 execution cycle을 계측 가능하다.

C에서 사용되는 표준 time.h 해당 함수들보다는 정확하나, 파이프라인, 수퍼스칼라등등 요인때문에 이 역시 완벽한 숫자는 아니다.

다음과 같은 inline assembly function을 사용한다.

static __inline__ long long getticks( void )

{

unsigned a, d;

asm volatile("rdtsc" : "=a" (a), "=d" (d));

return ((long long)a) | (((long long)d) << 32);

}

Digital Signal Processing System Design

Homework 4

Due: April 15, before the class

1. 1K by 1K matrix 형태의 real data가 x[n,m]이 있다.  이것을 2D FFT해서 X[j,k]를 얻으려 한다.

---------------->  row

|

|

|

v  column

입력 데이타, x[n,m], 로는 row 방향으로 주기가 64, column 방향으로는 주기가 32인 square wave를 준다 (크기 -0.5와 +0.5 사이). 즉 이것을 visualize한다면 벽돌 형태의 것이 주기적으로 튀어나온 형태가 될 것이다.  이 데이타 생성은 별로 중요한 것이 아니다.

(a) 2D FFT는 row 방향의 1D FFT를 모든 column 방향에 대해서 (병렬)수행하고, 다시 얻어진 데이타를 이용하여 column 방향의 1D FFT를 모든 row 방향에 대해서 (병렬)수행하면 된다.    이를 OpenMP 프로그램을 작성하고, 위의 첫번째 과정 (row 방향 1D FFT)과 두번째 과정에서 걸리는 시간을 각각 재시오.  Scheduling 방법은 가장 간단한 static으로 하시오.   Number of threads는 자기가 사용하는 CPU에서 지원하는 갯수만큼 이용하시오.

(b) 만약 첫번째 과정과 두번째 과정의 시간이 많이 차이가 난다면 왜 그런지 설명하시오.

(c) 시간이 더 많이 걸리는 과정을 개선하시오.  한 가지는 2D data를 transpose하는 방법이 될 것이다 (이 때 transpose방법을 잘 생각해야). 다른 더 좋은 방법도 생가할 수 있다.

(d) 여러가지 스케쥴링 방법을 이용하여 수행하여, 어떤 차이가 발생했는지 관찰하시오.  Static, dynamic, guided 를 적용해보고, chunk size도 1부터 크기를 키워보시오 (chunk size option은 세가지 정도).

(d) 본인이 사용하고 있는 CPU의 L1 cache size, L1 cache의 line size, L2 cache size등을 조사하시오.

(e) 입력이 real data임을 이용하여 위의 계산을 더 간단히 하도록 프로그램을 바꾸고, 시간이 얼마나 줄었는지 보시오.

2. 다음을 읽으시오.

W. Sung and S. K. Mitra, "Implementation of digital filtering algorithms using pipelined vector processors," Proc. IEEE, vol. 75, no. 9, pp. 1293-1302, Sep. 1987.

1. 다음의 프로그램을 OpenMP로 작성한다.  우선 자기가 가지고 있는 CPU의 최대 thread 개수를 파악한다.  이 thread 개수의 절반은      c[i] = a[i] + b[i]; 또 다른 절반은  d[i] = a[i] * b[i];의 계산에 할당하도록 프로그램을 작성하시오.

initalization 과정은 serial로 해도 된다.  굳이 이렇게 할 필요는 없겠지만, program 연습목적.

#define N     1000

main ()

{

int i;

float a[N], b[N], c[N], d[N];

/* Some initializations */

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

a[i] = i * 1.5;

b[i] = i + 22.35;

{

{

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

c[i] = a[i] + b[i];

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

d[i] = a[i] * b[i];

}

}

}

2. 다음의 프로그램을 run 하여 결과를 얻고, 그 까닭을 설명하시오.

#include <omp.h>

int  a, b, i, tid;

float x;

#pragma omp threadprivate(a, x)

main ()  {

/* Explicitly turn off dynamic threads */

omp_set_dynamic(0);

printf("1st Parallel Region:\n");

#pragma omp parallel private(b,tid)

{

tid = omp_get_thread_num();

a = tid;

b = tid;

x = 1.1 * tid +1.0;

printf("Thread %d:   a,b,x= %d %d %f\n",tid,a,b,x);

}  /* end of parallel section */

printf("************************************\n");

printf("Master thread doing serial work here\n");

printf("************************************\n");

printf("2nd Parallel Region:\n");

#pragma omp parallel private(tid)

{

tid = omp_get_thread_num();

printf("Thread %d:   a,b,x= %d %d %f\n",tid,a,b,x);

}  /* end of parallel section */

}

3. a[i]*b[i], i=0~100,000의 dot product (즉 sum for all i, a[i]*b[i] )를 구하는 프로그램을 다음의 방법으로 작성한 후 시간을 표시하시오.  위의 프로그램을 정확한 시간 측정을 위해서 일부러 여러번 돌릴 수도 있음.

(a) Critical section을 이용한다.

(b) 로컬 sum을 각각 구한 후에 나중에 이들 sum의 합을 구한다.  이 경우 프로그램의 변경(단지 pragma 넣기가 아니라)이 있을 수 있다.

(c) Reduction clause를 이용한다.

Digital Signal Processing System Design

Homework 5

Due: May 18, before the class

1. 다음의 FIR filtering 프로그램을 CUDA를 이용해서 작성하고 performance를 구하시오.

프로그램의 내용

입력 데이타: 길이 16M samples의 floating-point input data

출력 데이타: 길이 16M samples + 알파 (필터 길이 -1)의 floating-point data

filter:

Filter 길이 (M=) 8, 16, 32, 64, 128, 256에 대해서 실험한다.

필터의 계수는 Matlab을 이용하여 low-pass filter로 설계한다 (fp = 0.2 pi정도 되게).

입력 데이타는 random number generator를 이용하여 +1/-1 사이의 부동소수점 데이타를 만든다.

실험 목표: 최소한 두가지 이상의 FIR filtering program을 만들어서 실험을 한다. 그리고 성능을 구하는데, 최고 성능에 도달하지 못하게 하는 원인을 분석한다.

이 때 input/output time (PC<-> global memory), global memory<-> shared memory, 순수 computing time 등을 파악한다.

처음에 PC에서 프로그램을 짜서 성능을 파악하고 (시간), 다음에 CUDA 모의 실험기에서 돌린 후에 최종적으로 CUDA hardware에서 돌린다.

가능한 프로그램의 종류의 예:

- inner product method: 여러개의 thread들이 하나의 출력샘플을 계산한다.

- outer product method: 여러개의 출력샘플을 동시에 계산한다.

- 기타 mixed방법

- indirect computation method: filter tap의 길이가 클 때 좋은 방법으로 FFT를 이용한다. (보너스 점수)

- Linear phase filter이니 만치 symmetry를 이용하는 것이 도움이 되는지 살필 수 있다.

중요한 것: 평균 성능을 보여주어야 한다 (Gflops등). 또 성능 저하가 발생한 원인 등을 파악해야 한다.  프로그램을 돌리는 것 보다 성능 파악이 더 중요.  Thread, block의 갯수를 여러개로 실험해 본다.

*Inner product, outer product, indirect method는 전에 읽으라고 준 내 논문 (1987, Proceedings of the IEEE)에 있다.