Chapter 5

function y = circonv(x1, x2)

% Develops a sequence y obtained by the circular

% convolution of two equal-length sequences x1 and x2

L1 = length(x1);

L2 = length(x2);

if L1 ~= L2, error('Sequences of unequal length'), end

y = zeros(1,L1);

x2tr = [x2(1) x2(L2:-1:2)];

for k = 1:L1,

sh = circshift(x2tr', k-1)';

h = x1.*sh;

disp(sh);

y(k) = sum(h);

end

function [Y] = haar_1D(X)

% Computes the Haar transform of the input vector X

% The length of X must be a power-of-2

% Recursively builds the Haar matrix H

v = log2(length(X)) - 1;

H = [1 1;1 -1];

for m = 1:v,

A = [kron(H,[1 1]);

2^(m/2).*kron(eye(2^m),[1 -1])];

H = A;

end

Y = H*double(X(:));

% Program 5_1

% Illustration of DFT Computation

%

% Read in the length N of sequence and the desired

% length M of the DFT

N = input('Type in the length of the sequence = ');

M = input('Type in the length of the DFT = ');

% Generate the length-N time-domain sequence

u = [ones(1,N)];

% Compute its M-point DFT

U = fft(u,M);

% Plot the time-domain sequence and its DFT

t = 0:1:N-1;

stem(t,u)

title('Original time-domain sequence')

xlabel('Time index n'); ylabel('Amplitude')

pause

subplot(2,1,1)

k = 0:1:M-1;

stem(k,abs(U))

title('Magnitude of the DFT samples')

xlabel('Frequency index k'); ylabel('Magnitude')

subplot(2,1,2)

stem(k,angle(U))

title('Phase of the DFT samples')

xlabel('Frequency index k'); ylabel('Phase')

% Program 5_2

% Illustration of IDFT Computation

%

% Read in the length K of the DFT and the desired

% length N of the IDFT

K = input('Type in the length of the DFT = ');

N = input('Type in the length of the IDFT = ');

% Generate the length-K DFT sequence

k = 0:K-1;

V = k/K;

% Compute its N-point IDFT

v = ifft(V,N);

% Plot the DFT and its IDFT

stem(k,V)

xlabel('Frequency index k'); ylabel('Amplitude')

title('Original DFT samples')

pause

subplot(2,1,1)

n = 0:N-1;

stem(n,real(v))

title('Real part of the time-domain samples')

xlabel('Time index n'); ylabel('Amplitude')

subplot(2,1,2)

stem(n,imag(v))

title('Imaginary part of the time-domain samples')

xlabel('Time index n'); ylabel('Amplitude')

% Program 5_3

% Numerical Computation of Fourier transform Using DFT

k = 0:15; w = 0:511;

x = cos(2*pi*k*3/16);% Generate the length-16 sinusoidal sequence

X = fft(x); % Compute its 16-point DFT

XE = fft(x,512); % Compute its 512-point DFT

% Plot the frequency response and the 16-point DFT samples

plot(k/16,abs(X),'o', w/512,abs(XE))

xlabel('\omega/\pi'); ylabel('Magnitude')

% Program 5_4

% Linear Convolution Via the DFT

%

% Read in the two sequences

x = input('Type in the first sequence = ');

h = input('Type in the second sequence = ');

% Determine the length of the result of convolution

L = length(x)+length(h)-1;

% Compute the DFTs by zero-padding

XE = fft(x,L); HE = fft(h,L);

% Determine the IDFT of the product

y1 = ifft(XE.*HE);

% Plot the sequence generated by DFT-based convolution and

% the error from direct linear convolution

n = 0:L-1;

subplot(2,1,1)

stem(n,y1)

xlabel('Time index n');ylabel('Amplitude');

title('Result of DFT-based linear convolution')

y2 = conv(x,h);

error = y1-y2;

subplot(2,1,2)

stem(n,error)

xlabel('Time index n');ylabel('Amplitude')

title('Error sequence')

% Program 5_5

% Illustration of Overlap-Add Method

%

% Generate the noise sequence

R = 64;

d = rand(R,1)-0.5;

% Generate the uncorrupted sequence and add noise

k = 0:R-1;

s = 2*k.*(0.9.^k);

x = s + d';

% Read in the length of the moving average filter

M = input('Length of moving average filter = ');

% Generate the moving average filter coefficients

h = ones(1,M)/M;

% Perform the overlap-add filtering operation

y = fftfilt(h,x,4);

% Plot the results

plot(k,s,'r-',k,y,'b--')

xlabel('Time index n');ylabel('Amplitude')

legend('s[n]','y[n]')