ADSP Lab Manual

a)Linear convolution of two sequences

% Description : Script to find Linear convolution of two given sequences

% Prepare by clearing the screen

clc;clear all;close all;

% inputs: x[n], h[n], N the number of samples

xn = input('Enter the first sequence: '); %[1 1 0 1]
hn = input('Enter the second sequence: '); %[1 2 3 0]

% Specifying the position of samples in each sequence

n1 = -2 : 1;
n2 = -1 : 2;

% length and positions of the output sequence helps in plotting

ybegin = n1(1) + n2(1);
yend = n1(length(xn)) + n2(length(hn));
ny = ybegin : yend;

% Perform convolution
y = conv(xn,hn);

% display output sequence

disp('Linearly convoluted output sequence: ')
disp(y);

% Plot the input sequence 1

subplot(3,1,1);
stem(n1,xn);
xlabel('Time in samples');
VTU CPGSB Muddenhalli Page 1
 ADSP Lab Manual

ylabel('Amplitude');
title('The first input sequence');

% Plot the input sequence 2

subplot(3,1,2);
stem(n2,hn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The second input sequence');
% Plot the output sequence
subplot(3,1,3);
stem(ny,y);
xlabel('Time in samples');
ylabel('Amplitude');
title('The linear convoluted output sequence');
(OR)

Program

clc;

x1=input('enter the first sequence');

subplot(3,1,1);
stem(x1);
ylabel('amplitude');
title('plot of the first sequence');
x2=input('enter 2nd sequence');
subplot(3,1,2);
stem(x2);
ylabel('amplitude');
title('plot of 2nd sequence');
f=conv(x1,x2);

VTU CPGSB Muddenhalli Page 2

 ADSP Lab Manual

disp('output of linear conv is');

disp(f);
xlabel('time index n');
ylabel('amplitude f');
subplot(3,1,3);
stem(f);
title('linear conv of sequence');

Results:

Enter the first sequence: [1 1 0 1]

Enter the second sequence: [1 2 3 0]

Linearly convoluted output sequence:

1 3 5 4 2 3 0

VTU CPGSB Muddenhalli Page 3

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 4

 ADSP Lab Manual

(b) Circular convolution

% Description : Script to find Circular convolution of two given sequences

% Prepare by clearing the screen

clc;clear all;close all;

% inputs: x[n],h[n]
xn = input('Enter the first sequence: '); % [1 1 0 1]
hn = input('Enter the second sequence: '); % [1 2 3 0]

% To adjust for different length sequences

N1 = length(xn);
N2 = length(hn);
N = max(N1, N2);

xn = [xn, zeros(1,N - N1)];

hn = [hn, zeros(1,N - N2)];
y = zeros(1,N);

% Generate circular convolution of two sequences

for i = 0:N-1
for j = 0:N-1
k = mod((i-j),N);
y(i+1) = y(i+1) +(xn(j+1) * hn(k+1));
end
end
n = 1:N;

% display output sequence

disp('Circularly convoluted output sequence: ')
disp(y);

VTU CPGSB Muddenhalli Page 5

 ADSP Lab Manual

% Plot the input sequence 1

subplot(3,1,1);
stem(n,xn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The first input sequence');

% Plot the input sequence 2

subplot(3,1,2);
stem(n,hn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The second input sequence');

% Plot the output sequence

subplot(3,1,3);
stem(n,y);
xlabel('Time in samples');
ylabel('Amplitude');
title('The circularly convoluted output sequence');

VTU CPGSB Muddenhalli Page 6

 ADSP Lab Manual

(OR)
clc;
clear all;
x1=input('enter the first sequence');
x2=input('enter the second sequence');
n1 = length(x1);
n2 = length(x2);
subplot(3,1,1);
stem(x1,'filled');
title('plot of first sequence');
subplot(3,1,2);
stem(x2,'filled');
title('plot the second sequnce');
y1=fft(x1,n);
y2=fft(x2,n);
y3=y1.*y2;
y=ifft(y3,n);
disp('the circular convolution result is ......');
disp(y);
subplot(3,1,3);
stem(y,'filled');

title('plot of circularly convoluted sequence'

Results:

VTU CPGSB Muddenhalli Page 7

 ADSP Lab Manual

Enter the first sequence: [1 1 0 1]

Enter the second sequence: [1 2 3 0]

Circularly convoluted output sequence:

3 6 5 4

VTU CPGSB Muddenhalli Page 8

 ADSP Lab Manual

(c) Linear convolution using circular convolution

clc;
clear
all;
x1=input(‘enter the first sequence’);
x2=input(‘enter the second sequence’);
n=input(‘enter the no of points of the
dft’); subplot(3,1,1);
stem(x1,’filled’);
title(‘plot of first sequence’);

subplot(3,1,2);
stem(x2,’filled’);
title(‘plot the second sequnce’);
n1 =length(x1);
n2 = length(x2);
m = n1+n2-1; % Length of linear convolution
x = [x1 zeros(1,n2-1)]; % Padding of zeros to make it of
% length m
y = [x2 zeros(1,n1-1)];
x_fft = fft(x,m);
y_fft = fft(y,m);
dft_xy =
x_fft.*y_fft;

y=ifft(dft_xy,m);
disp(‘the circular convolution result is ......’);
disp(y);
subplot(3,1,3);
stem(y,’filled’);
title(‘plot of circularly convoluted sequence’);
VTU CPGSB Muddenhalli Page 9
 ADSP Lab Manual

Output

enter the first sequence[1 2 1 2 1 2]

enter the second sequence[1 2 3 4]

enter the no of points of the dft 6

the circular convolution result is ......

1.0000 4.0000 8.0000 14.0000 16.0000 14.0000 15.0000 10.0000 8.0000

VTU CPGSB Muddenhalli Page 10

 ADSP Lab Manual

2)(a) DFT of Given Sequences

% Description : Script to find DFT of a sequence

% Prepare by clearing the screen

clc;clear all;close all;

% Take the input sequence

xn = input('Enter the input sequence: '); % [1 3 -1 -2]
N = length(xn);

% Initialize output sequence with zeros

Xk = zeros(1,N);

% Compute DFT{x[n]} using formula

for k = 0:N-1
for n = 0:N-1
Xk(k+1) = Xk(k+1) + xn(n+1) * exp((-j*2*pi*n*k)/N);
end
end

% Display output sequence

disp('The DFT X(k) for given sequence is: ');
disp(Xk); % [1 2-5i -1 2+5i]

n = 0:N-1;

% Plot the input sequence

subplot(3,1,1);
stem(n,xn);
xlabel('Samples');

VTU CPGSB Muddenhalli Page 11

 ADSP Lab Manual

ylabel('Amplitude');
title('INPUT SEQUENCE');

% Plot the magnitude of output

subplot(3,1,2);
stem(n,abs(Xk));
xlabel('Samples');
ylabel('Amplitude');
title('MAGNITUDE');
% Plot the degree of output sequence
subplot(3,1,3);
stem(n,angle(Xk));
xlabel('Samples');
ylabel('Amplitude');
title('DEGREE');
Results:

Enter the input sequence: [1 3 -1 -2]

The DFT X(k) for given sequence is:

1.0000 2.0000 - 5.0000i -1.0000 + 0.0000i 2.0000 + 5.0000i

VTU CPGSB Muddenhalli Page 12

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 13

 ADSP Lab Manual

(b) IDFT of Given sequences

% Description : Script to find IDFT of a sequence

% Prepare by clearing the screen

clc;clear all;close all;

% Take the input sequence

Xk = input('Enter the input sequence: '); % [1 2-5i -1 2+5i]
N = length(Xk);

% Initialize output sequence with zeros

xn = zeros(1,N);

% Compute IDFT{x[n]} using formula

for n = 0:N-1
for k = 0:N-1
xn(n+1) = xn(n+1) + Xk(k+1) * exp((j*2*pi*n*k)/N);
end
end

xn = xn/N;

% Display output sequence

disp('The IDFT x(n) for given sequence is: ');
disp(round(real(xn))); % [1 3 -1 -2]

n = 0:N-1;

% Plot the input sequence

subplot(3,1,1);
stem(n,abs(Xk));

VTU CPGSB Muddenhalli Page 14

 ADSP Lab Manual

xlabel('Samples');
ylabel('Amplitude');
title('INPUT SEQUENCE');
% Plot the magnitude of output
subplot(3,1,2);
stem(n,round(real(xn)));
xlabel('Samples');
ylabel('Amplitude');
title('MAGNITUDE');
% Plot the phase of output sequence
subplot(3,1,3);
stem(n,angle(xn));
xlabel('Samples');
ylabel('Amplitude');
title('DEGREE');
Results:

Enter the input sequence: [1 2-5i -1 2+5i]

The IDFT x(n) for given sequence is:

1 3 -1 -2

VTU CPGSB Muddenhalli Page 15

 ADSP Lab Manual

(C) Circular convolution using DFT and IDFT

% Description : Script to find circular convolution of two given sequences

% using DFT and IDFT

% Prepare by clearing the screen

clc;clear all;close all;

% inputs: x[n], h[n]

xn = input('Enter the first sequence: '); % [1 1 0 1]
hn = input('Enter the second sequence: '); % [1 2 3 0]

% Specifying the position of samples in each sequence

% To adjust for different length sequences
N1 = length(xn);
N2 = length(hn);
N = max(N1,N2);
xn = [xn, zeros(1,(N-N1))];
hn = [hn, zeros(1,(N-N2))];

% Compute DFT of both sequences

Xk = fft(xn);
Hk = fft(hn);

% Compute output and IDFT

Yk = Xk .* Hk;
yn = ifft(Yk);

% Note: The output of convolution is obtained as a double since fft is used

% to match the output with normal convolution program abs, uint8
% conversions are used
disp('The Circularly convoluted output(using DFT & IDFT) is: ');

VTU CPGSB Muddenhalli Page 16

 ADSP Lab Manual

disp(yn);

n = 1:N;

% Plot the input sequence 1

subplot(3,1,1);
stem(n,xn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The first input sequence');

% Plot the input sequence 2

subplot(3,1,2);
stem(n,hn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The second input sequence');

% Plot the output sequence

subplot(3,1,3);
stem(n,yn);
xlabel('Time in samples');
ylabel('Amplitude');
title('The Circular convoluted output(using DFT & IDFT) sequence');

Results

Enter the first sequence: [1 1 0 1]

Enter the second sequence: [1 2 3 0]

The Circularly convoluted output(using DFT & IDFT) is:

3 6 5 4

VTU CPGSB Muddenhalli Page 17

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 18

 ADSP Lab Manual

3) comparison of the DFT and DCT (in terms of energy compactness)

(a)Let X[k]=DFT{x[n]}. For various values of L , set to zero “high frequency coefficients ”
X[64-1]=………….X[64]=……X[64+L]=0 and take the inverse DFT .plot the results

(b)Let XDCT[k]=DCT(x[n]).For the same values of L ,set to zero “high frequency coefficient ”
XDCT[127-L]=….. XDCT[127]. Take the inverse DCT for each case and compare the
reconstruction with the previous case.

% Description : Script to compare DCT and DFT in terms of energy

% compactness

% Prepare by clearing the screen

clc;clear all;close all;

% Generate the sequence x[n] = n-64 for n = 0 ..... 127

n = 0:127;
x = n-64;
z = x;

% X[k] = DFT(x[n])
Xk = fft(x);

% For various values of L set to zero "High Frequency coefficients"

% Xk[64-L] =...Xk[64]=...Xk[64+L] = 0 and find inverse DFT
L = 32;
Xk((64 - L):(L + 64)) = 0;
Yf = ifft(Xk);

% Plot input Sequence

subplot(2,2,1);

VTU CPGSB Muddenhalli Page 19

 ADSP Lab Manual

plot(x);
title('INPUT SEQUENCE');

% Plot IFFT Output Sequence

subplot(2,2,2);
plot(real(Yf(1:128)));
title('IFFT OUTPUT SEQUENCE');

% XDCT[k] = DCT(x[n])
XDCTk = dct(z);

% For same valuse of L, set to zero "High frequency Coefficients"

% XDCT[127-L]=...XDCT[127] and take IDCT
XDCTk((128-L):(128+L)) = 0;
Xl = idct(XDCTk);

% Plot input squence

subplot(2,2,3);
plot(z);
title('INPUT SEQUENCE');

% Plot the IDCT Output Sequence

subplot(2,2,4);
plot(Xl);
title('IDCT OUTPUT SEQUENCE');

VTU CPGSB Muddenhalli Page 20

 ADSP Lab Manual

(OR)
Program:

clc;
clear all;
close all;
n=0:127;
x=n-64;
z=x;
xk=fft(x);
L=32;
xk((64-L):(L+64))=0;
yf=ifft(xk);
figure;
subplot(2,2,1);
plot(x);
title('input sequence');
subplot(2,2,2);
plot(real(yf(1:128)));
title('IFFT output sequence');
xk=dct(z);
xk((128-L):(128+L))=0;
x1=idct(xk);
subplot(2,2,3);
plot(x);
title('input sequence');
subplot(2,2,4);
plot(x1);
title('IDCT output sequence');

VTU CPGSB Muddenhalli Page 21

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 22

 ADSP Lab Manual
4) Power Spectral Density
AIM: To verify Power Spectral Density
EQUIPMENTS:
Operating System - Windows family
Constructor - Simulator
Software - Matlab
THEORY: The power spectral density(P.S.D) is a measurement of the energy at various
frequencies
In the previous section we saw how to unwrap the FFT and get back the sine and cosine
coefficients. Usually we only care how much information is contained at a particular frequency
and we don’t really care whether it is part of the sine or cosine series. Therefore, we are
interested in the absolute value of the FFT coefficients. The absolute value will provide you with
the total amount of information contained at a given frequency, the square of the absolute value
is considered the power of the signal. Remember that the absolute value of the Fourier
coefficients are the distance of the complex number from the origin. To get the power in the
signal at
each frequency (commonly called the power spectrum) you can try the following commands.
>> N = 8; %% number of points
>> t = [0:N-1]’/N; %% define time
>> f = sin(2*pi*t); %%define function
>> p = abs(fft(f))/(N/2); %% absolute value of the fft
>> p = p(1:N/2).^2 %% take the positve frequency half, only
This set of commands will return something much easier to understand, you should get 1 at a
frequency of 1 and zeros everywhere else. Try substituting cos for sin in the above commands,
you should get the same result. Now try making >>f = sin(2*pi*t) + cos(2*pi*t). This change
should result in twice the power contained at a frequency of 1.
Thus far we have looked at data that is defined on the time interval of 1 second, therefore we
could interpret the location in the number list as the frequency.

VTU CPGSB Muddenhalli Page 23

 ADSP Lab Manual

PROGRAM:
%Power spectral density
t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
subplot(2,1,1);
% figure(1);
plot(1000*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)');
Y = fft(y,512);
%The power spectral density, a measurement of the energy at various frequencies, is:
Pyy = Y.* conj(Y) / 512;
f = 1000*(0:256)/512;
subplot(2,1,2);
% figure(2);
plot(f,Pyy(1:257))
title('Frequency content of y');
xlabel('frequency (Hz)');

OUTPUT

VTU CPGSB Muddenhalli Page 24

 ADSP Lab Manual

5) IMPLEMENTATION OF DECIMATION PROCESS

AIM: To implementation of decimation of given sequence by factor M.
EQUIPMENTS:
Operating System – Windows Family
Constructor - Simulator
Software - MATLAB
THEORY :
Decimation is the process of reducing the sampling frequency of a signal to a lower sampling
frequency that differs from the original frequency by an integer value. Decimation also is known
as down-sampling. The lowpass filtering associated with decimation removes high-frequency
content from the signal to accommodate the new sampling frequency.
Decimation is useful in applications in which the Nyquist frequency of a signal is much higher
than the highest frequency of the signal. Decimation filters help you remove the excess
bandwidth and reduce the sampling frequency of the signal. Decimation filters also help you
reduce the computational resources required for processing and storing the signal. During the
analog-to-digital (A/D) conversion process, decimation filters also can reduce the variance of
quantization noise in a signal and maintain the signal power, thus improving the signal-to-noise
ratio (SNR).
The following figure shows a typical M-fold decimation filter, where M is the integer value by
which you want to decrease the sampling frequency. This filter contains a lowpass FIR filter
H(z). This lowpass FIR filter is an anti-aliasing filter followed by anM-fold decimator. The
decimator passes every Mth sample and discards the other samples. After this operation, the
decimation filter changes the sampling frequency fs of the input signal x(n) to a new sampling
frequency fs/M. The decimation filter then returns an output signal y(n) with the new sampling
frequency.

VTU CPGSB Muddenhalli Page 25

 ADSP Lab Manual

equation for decimation

To prevent aliasing, this system uses the lowpass filter H(z) before the M-fold decimator to
suppress the frequency contents above the frequency fs/(2M), which is the Nyquist frequency of
the output signal. This system produces the same results as an analog anti-aliasing filter with a
cutoff frequency of fs/(2M) followed by an analog-to-digital (A/D) converter with a sampling
frequency of fs/M. Because the system shown in the figure above is in the digital domain, H(z) is
a digital anti-aliasing filter.

IMPLEMENTATION OF INTERPOLATION PROCESS

THEORY:
The process of increasing the sampling rate is called interpolation. Interpolation is up sampling

followed by appropriate filtering. obtained by interpolating , is generally represented

as:
The simplest method to interpolate by a factor of L is to add L-1 zeros in between the samples,
multiply the amplitude by L and filter the generated signal, with a so-called anti-imaging low
pass filter at the high sampling frequency.

Program
clc;
clear all;
close all;
t=0:0.01:1;
d=input('enter number D/I');
x=sin(2*pi*5*t);
y=decimate(x,d);
subplot(3,1,1);
stem(x(1:25));
title('original input signal');
xlabel('samples');
ylabel('amplitude');

VTU CPGSB Muddenhalli Page 26

 ADSP Lab Manual

subplot(3,1,2);
stem(y);
title('decimated output');
xlabel('samples');
ylabel('amplitude');
y1=interp(y,d);
subplot(3,1,3);
stem(y1(1:25));
title('interpolation output');
xlabel('samples');
ylabel('amplitude');

Results:

enter number D/I 5

VTU CPGSB Muddenhalli Page 27

 ADSP Lab Manual

6) Image Denoising
The denoising method described for the one-dimensional case applies also to images and applies
well to geometrical images. The two-dimensional denoising procedure has the same three steps
and uses two-dimensional wavelet tools instead of one-dimensional ones. For the threshold
selection, prod(size(y)) is used instead of length(y) if the fixed form threshold is used.

Generate a noisy image.

Program:

load woman
init = 2055615866; rng(init);
x = X + 15*randn(size(X));

[thr,sorh,keepapp] = ddencmp('den','wv',x);
thr =557.9838

xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);
figure('Color','white')
colormap(pink(255)), sm = size(map,1);
image(wcodemat(X,sm)), title('Original Image')
figure('Color','white')
colormap(pink(255))
image(wcodemat(x,sm)), title('Noisy Image')
image(wcodemat(xd,sm)), title('De-Noised Image')

VTU CPGSB Muddenhalli Page 28

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 29

 ADSP Lab Manual

VTU CPGSB Muddenhalli Page 30

 ADSP Lab Manual

7)Haar Wavelet Transform for Image compression

clc; clear; close all;

disp('2D-DWT Image file *Check workspace')
% disp('2D-DWT using dwt function')
% process sample image file
y=imread('Tulips.jpg);
y=rgb2gray(y);
x=single(y); %very important if not x remain in 8bit

h0=[0.7071 0.7071];
h1=[0.7071 -0.7071];
% disp('2D-DWT Decomposition')
[A,H,V,D]=dwt2(x,h0,h1);

[m,n]=size(x); %[m=row n=column]

for i=1:m/2
for j=1:n/2
% disp('2D-DWT Decomposition of A using formula')
Ad(i,j)=(x(2*i-1,2*j-1)+x(2*i-1,2*j)+x(2*i,2*j-1)+x(2*i,2*j))*h0(1)*h0(1);

% disp('2D-DWT Decomposition of H using formula')

Hd(i,j)=(x(2*i,2*j-1)+x(2*i,2*j))*h0(1)*h0(1)+(x(2*i-1,2*j-1)+x(2*i-1,2*j))*h0(1)*h1(2);

% disp('2D-DWT Decomposition of V using formula')

Vd(i,j)=(x(2*i-1,2*j)+x(2*i,2*j))*h0(1)*h0(1)+(x(2*i-1,2*j-1)+x(2*i,2*j-1))*h0(1)*h1(2);

% disp('2D-DWT Decomposition of D using formula')

Dd(i,j)=(x(2*i,2*j))*h0(1)*h0(1)+(x(2*i-1,2*j-1))*h1(2)*h1(2)+(x(2*i-1,2*j)+x(2*i,2*j-
1))*h0(1)*h1(2);

VTU CPGSB Muddenhalli Page 31

 ADSP Lab Manual

end
end

figure
imshow(y)
figure
Atemp=uint8(Ad);
imshow(Atemp)

OutPut

Original Image

VTU CPGSB Muddenhalli Page 32

 ADSP Lab Manual

Compressed Image

VTU CPGSB Muddenhalli Page 33