Dsp Simulation Assignment Using Octave

DSP Simulation Assignment using Octave V. A. Amarnath November 17, 2009

1

Contents I Generate following sequences

4

1 Unit Impulse Function

4

2 Unit Step Function

5

3 Ramp Function

6

4 Exponential Sequence

7

5 Sinusoidal Sequence

8

6 given x[n]

9

II Linear Convolution

10

III Circular Convolution

13

IV FFT Algorithms

16

V Symmetry Property of DFT

20

VI Scaling Property of DFT

22

VII Overlap Save Method

24

VIII Overlap Add Method

27

IX Linear Convolution using DFT

30

X Analyse an audio le

32

2

List of Figures 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

Unit Impulse Response . . . . . . . . . . . . Unit Step Function . . . . . . . . . . . . . . Ramp Function . . . . . . . . . . . . . . . . Exponential Sequence . . . . . . . . . . . . Sinusoidal Sequence . . . . . . . . . . . . . Given x[n] . . . . . . . . . . . . . . . . . . . x[n] . . . . . . . . . . . . . . . . . . . . . . . y[n] . . . . . . . . . . . . . . . . . . . . . . . z[n] = x[n]*y[n] . . . . . . . . . . . . . . . . x[n] . . . . . . . . . . . . . . . . . . . . . . . y[n] . . . . . . . . . . . . . . . . . . . . . . . z[n] = x[n]~y[n] . . . . . . . . . . . . . . . . x[n] . . . . . . . . . . . . . . . . . . . . . . . real{DF T (x[n])} . . . . . . . . . . . . . . . imag{DF T (x[n])} . . . . . . . . . . . . . . Symmetry Property of DFT . . . . . . . . . Scaling Property of DFT . . . . . . . . . . . x[n] . . . . . . . . . . . . . . . . . . . . . . . h[n] . . . . . . . . . . . . . . . . . . . . . . Result of Overlap Save Method using l = 4 x[n] . . . . . . . . . . . . . . . . . . . . . . . h[n] . . . . . . . . . . . . . . . . . . . . . . Result of Overlap Add Method using l = 3 x[n] . . . . . . . . . . . . . . . . . . . . . . . y[n] . . . . . . . . . . . . . . . . . . . . . . . x[n] * y[n] . . . . . . . . . . . . . . . . . . . Audio Sequence . . . . . . . . . . . . . . . . Frequency Spectrum . . . . . . . . . . . . .

3

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 5 6 7 8 9 11 11 12 14 14 15 18 18 19 21 23 25 25 26 28 28 29 30 31 31 32 33

Part I

Generate following sequences 1

Unit Impulse Function

Code clf; y = [linspace(0,0,10),1,linspace(0,0,10)]; l = length(y); x = [-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([-9, 9, -2, 2]); grid on; hold on; print impulse.jpg -djpg;

Figure

Figure 1: Unit Impulse Response

4

2

Unit Step Function

Code clf; y = [linspace(0,0,9),linspace(1,1,10)]; l = length(y); x = [-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([-9, 9, -2, 2]); grid on; hold on; print step.jpg -djpg;

Figure

Figure 2: Unit Step Function

5

3

Ramp Function

Code clf; y = [linspace(0,0,10),linspace(0,10,11)]; l = length(y); x = [-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([-10, 10, -10, 10]); grid on; hold on; print ramp.jpg -djpg;

Figure

Figure 3: Ramp Function

6

4

Exponential Sequence

Code clf; y = [linspace(0,5,101),exp(linspace(0,5,101))]; l = length(y); x = 0.1*[-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([0, 11, -1, 150]); grid on; hold on; print expo.jpg -djpg;

Figure

Figure 4: Exponential Sequence

7

5

Sinusoidal Sequence

Code clf; y = [linspace(0,0,101),sin(linspace(0,10,101))]; l = length(y); x = [-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([-100, 100, -1, 1]); grid on; hold on; print sine.jpg -djpg;

Figure

Figure 5: Sinusoidal Sequence

8

6

given x[n]

x[n] = 5 cos{ 23 π n + 43 π } + 4 cos{ 53 π n} - sin{ 12 π n}

Code clf; function p = z(n) p = 5*cos(1.5*pi*n+.75*pi)+4*cos(.6*pi*n)-sin(.5*pi*n); endfunction y = [linspace(0,0,11),z(linspace(0,10,11))]; l = length(y); x = [-(l-1)/2:1:(l-1)/2]; stem (x,y); axis ([-10, 10, -10, 10]); grid on; hold on; print xn.jpg -djpg;

Figure

Figure 6: Given x[n]

9

Part II

Linear Convolution Code function [res, num] = linconv(a,b,n1,n2) r=zeros(length(a),length(b)); for m = 1 : length(a) for n = 1 : length (b) r(n,m) = a(m)*b(n); end end for k = 2 : length(a)+length(b) res(k-1) = 0; for m = 1 : length(a) for n = 1 : length(b) if (m + n) == k res(k-1) = r(n,m) + res(k-1); end end end end num = n1+n2:1:length(a)+length(b)+n1+n2-2; end

10

Figures

Figure 7: x[n]

Figure 8: y[n]

11

Figure 9: z[n] = x[n]*y[n]

12

Part III

Circular Convolution Code function all = circonv(a,b) n = max(length(a),length(b)); if length(a) ~= n for k = length(a)+1 : n a(k) = 0; end; end; if length(b) ~= n for k = n - length(b) : n b(k) = 0; end; end; t = zeros(n); temp = zeros(n,n); for k = 1 : n t = shift(b,k-1); for l = 1 : n temp(k,l)=t(l); end; end; all = zeros(n); all = (a*temp); end;

13

Figures

Figure 10: x[n]

Figure 11: y[n]

14

Figure 12: z[n] = x[n]~y[n]

15

Part IV

FFT Algorithms DIT FFT Algorithm Code function z = ditfft(x) N = length(x); n = log2(N); clear t1; z = zeros(1, N); for p = 1 : N t = dec2bin(p-1,n); for q = 1 : n t1(q) = t(n-q+1); end; t = t1; z (p) = x(bin2dec(t)+1); end; for p = 1:n m = 2^(p); w= exp(-i*2*pi/m); for l = 1: N/m temp = z(1+(l-1)*m:l*m); for k= 1:m/2 z(k+(l-1)*m) = temp(k)+temp(m/2+k)*w^(k-1); z(m/2+k+(l-1)*m) = temp(k)-temp(m/2+k)*w^(k-1); end; end; end; end;

16

DIF FFT Algorithm Code function y = diffft(x) len = length(x); n = log2(len); g = x; for p = n:-1:1 lt = 2^(p); Wn= exp(-i*2*pi/lt); for l = 1: len/lt temp = g(1+(l-1)*lt:l*lt); for k= 1:lt/2 g(k+(l-1)*lt) = temp(k)+temp(lt/2+k); g(lt/2+k+(l-1)*lt)=(temp(k)-temp(lt/2+k))*Wn^(k-1); end; end; end; tmp=g; for p = n:-1:1 lt = 2^(p); m=1; for l = 1: len/lt for k = 1:2:lt y(m) = tmp(k+lt*(l-1)); y(lt/2+m)=tmp(k+lt*(l-1)+1); m=m+1; end; m=l*lt+1; end; tmp=y; end; end;

17

Figures

Figure 13: x[n]

Figure 14: real{DF T (x[n])}

18

Figure 15: imag{DF T (x[n])}

19

Part V

Symmetry Property of DFT Code x = [1 2 3 4 4 3 2 1]; N = length(x); n= 0:N-1; y = fft(x); z = [y y]; for k = 1:N sym(k) = conj(z(N+1-(k-1))); end; subplot(2,2,1) stem(n,abs(y)); grid on; xlabel('time'); ylabel('amplitude'); title('magnitude of X(k)'); subplot(2,2,2) stem(n,angle(y)); grid on; xlabel('time'); ylabel('amplitude'); title('phase of X(k)'); subplot(2,2,3) stem(n,abs(sym)); grid on; xlabel('time'); ylabel('amplitude'); title('magnitude of X*(N-k)'); subplot(2,2,4) stem(n,angle(sym)); grid on; xlabel('time'); ylabel('amplitude'); title('phase of X*(N-k)');

20

Figure

Figure 16: Symmetry Property of DFT

21

Part VI

Scaling Property of DFT Code x=[8 7 6 5 4 3 2 1]; len = length(x); y=0:1:len-1; subplot(3,1,1) stem(y,x,'fillled'); grid on; axis([0 8 0 10]); xlabel('Time '); ylabel('Amplitude'); title('Input sequence'); z=fft(x,len); subplot(3,1,2) stem(y,abs(z),'fillled'); axis([0 8 0 40]); grid on xlabel('Time '); ylabel('Amplitude'); title('8 point'); z=fft(x,4*len); y=0:1:4*len-1; subplot(3,1,3); stem(y,abs(z),'fillled'); axis([0 32 0 40]); grid on xlabel('Time '); ylabel('Amplitude'); title('32 point');

22

Figure

Figure 17: Scaling Property of DFT

23

Part VII

Overlap Save Method Code function res = oversave(x, h,l) m = length (h); L = length (x); y = [zeros(1, m-1), x]; r = rem (length (y), l); y = [y, zeros(1, r)] if r == 0 y = [y, zeros(1, l/2)]; end; xn = length (y); h1 = [h, zeros(1, l - m)]; res = zeros (1, L + m -1); k = 1; pos = 0; p = 0; count = 0; while k <= length(res) temp = zeros (1,l); for p = 0 : l-1 temp(p+1) = y(k+p); end; t = circonv (temp, h1);1 p = 1; while p <= l - m + 1 res (p + pos) = t (p + m - 1); p = p + 1; end; pos = pos + l - m + 1; k = k + l - m + 1; end; end; 1 circonv function used from Q. 3

24

Figures

Figure 18: x[n]

Figure 19: h[n]

25

Figure 20: Result of Overlap Save Method using l = 4

26

Part VIII

Overlap Add Method Code function res = overadd(x,h,l) n = length(x); m = length(h); r = rem(n, l); if r ~= 0 x = [x, zeros(1, r)]; end; xn = length(x); k = 1; h = [h, zeros(1, l - 1)]; res = zeros(1, n + m -1); rp = 1; while k <= xn t = zeros(1, l + m - 1); for p = 1 : l t(p) = x(k + p - 1); end; t = circonv(t, h);2 k = k + l; for p = 1 : l + m - 1 res(rp) = res(rp) + t(p); rp = rp +1; if rp > n + m - 1 break; end; end; rp = rp - m + 1; end; end; 2 circonv

function used from Q. 3

27

Figures

Figure 21: x[n]

Figure 22: h[n]

28

Figure 23: Result of Overlap Add Method using l = 3

29

Part IX

Linear Convolution using DFT Code function res = lindft(x, y) n1 = length(x); n2 = length(y); x1 = [x, zeros(1,n2-1)]; y1 = [y, zeros(1,n1-1)]; res = ifft(fft(x1).*fft(y1)); end;

Figures

Figure 24: x[n]

30

Figure 25: y[n]

Figure 26: x[n] * y[n]

31

Part X

Analyse an audio le Code x = wavread(`mohanlal.wav'); stem (x); y = fft(x); N = length(x); z = zeros(1,N/2); for k = 1:N/2 z(k) = y(k); end c = linspace(0,pi,N/2); c = c/(2*pi); c = c*8000; stem(c,abs(z)); grid on;

Figures

Figure 27: Audio Sequence

32

Figure 28: Frequency Spectrum

33