1. Write a MATLAB program to find the N point DFT of any finite discrete time sequence using direct computation (without using inbuilt commands).
a. Evaluate the DFT for x[n]= [1 2 3 4 4 3 2 1] and plot its magnitude and phase response
%MRINALINI
%AM.EN.U4ELC20022
x_n= [1 2 3 4 4 3 2 1] ;
N = 8;
L = 8 ;
zeros(1,N);
xk = [];
for k = 0:N-1
x_k=0;
for n=0:N-1
x_k =x_k+x_n(n+1)*exp(-i*k*2*pi*n/N);
end
xk= [xk x_k];
end
n = 0:7;
XK = [xk]
magn = abs(XK)
phase = angle(XK)
figure(1)
subplot(2,1,1)
stem(n,magn)
xlabel('n')
ylabel('magnitude X(K)')
title('Magnitude plot')
subplot(2,1,2)
stem(n,phase)
xlabel('n')
ylabel('Phase')
title('Phase Plot')
b. Verify the output by taking the inverse DFT.
zeros(1,N-1);
xn = [];
for n = 0:N-1
xx=0;
for k=0:N-1
xx = xx+(1/N)*(XK(k+1)*exp(i*n*k*2*pi/N));
end
xn = [xn xx];
end
X = [xn]
n = 0:7;
subplot(2,1,1)
stem(n,xn)
xlabel('n')
ylabel('x(n)')
title('Given sequence')
subplot(2,1,2)
stem(n,real(X))
xlabel('n')
ylabel('x(n)')
title('Sequence after DFT and IDFT')
Write a MATLAB program to find 4-point discrete Fourier transform of the signal x[n]={1,-1,2,-2}(Use the command dftmtx). Compute the IDFT of the result obtained.
%MRINALINI 20022
xn = [1 -1 2 -2];
N = 4;
b = dftmtx(N);
xk = b.*xn;
XK = [xk]
d = conj(dftmtx(N));
x = (1/N)*(d)*(XK)
3. Compute a 32-point DFT for the above problem. Plot the discrete magnitude and phase spectrum and comment on the frequency resolution compared to 4-point DFT.
%MRINALINI - 22
x= [1 -1 2 -2 ] ;
n1 = length(x);
L = 32;
N = 32;
x_n = [x zeros(1,(N-n1))];
zeros(0,N-1);
xk = [];
for k = 0:N-1
x_k = 0;
for n=0:N-1
x_k =x_k+x_n(n+1)*exp(-1i*k*2*pi*n/N);
end
xk = [xk x_k];
end
n = 0:31;
x_k_final = [xk]
y = abs(x_k_final)
figure(1)
subplot(2,1,1)
stem(n,y)
xlabel('n')
ylabel('magnitude')
title('magnitude plot - Mrinalini 22')
z = angle(x_k_final)
subplot(2,1,2)
stem(n,z)
xlabel('n')
ylabel('phase')
title('phase plot - Mrinalini 22')
1. Let x[n] = 10(0.8) n defined in the interval 0 ≤ n ≤ 10.
(a) Sketch x[((n + 4))11], that is, a circular shift by 4 samples toward the left.
%Mrinalini
% AM.EN.U4ELC20022
a = 0:10;
y_n = (10)*(0.8).^a
m=-4;
N=length(a);
k_shift=mod(a-m,N);
subplot(2,1,1)
stem(a,y_n)
xlabel('time')
ylabel('magnitude')
title("x_n")
subplot(2,1,2)
y2_n=y_n(k_shift+1)
stem(a,y2_n)
title("nshift")
xlabel('time')
ylabel('magnitude')
(b) Sketch x[((n − 3))15], that is, a circular shift by 3 samples toward the right, where x[n] is assumed to be a 15-point sequence. (Hint: n1 = mod(n-m,N); %y(n) = x((n-m) mod N))
%Mrinalini
% AM.EN.U4ELC20022
a = 0:14;
y_n = (10)*(0.8).^a
m=3;
N=length(a);
k_shift=mod(a-m,N);
subplot(2,1,1)
stem(a,y_n)
xlabel('time')
ylabel('magnitude')
title("x_n")
subplot(2,1,2)
y2_n=y_n(k_shift+1)
stem(a,y2_n)
title("nshift")
xlabel('time')
ylabel('magnitude')
(c) Write a MATLAB function called cshift(n,m,N) to implement circular shifting operations done in part (a) and (b)
function y = cshift(n,m,N)
n_shift=mod(n-m,N)
stem(n,n_shift)
2. Let x[n] = {1, 2, 2} and h[n] = {1, 2, 3, 4}. (Use ‘input’ command to give the data).
(a) Compute the 4 - point circular convolution by time domain approach (by using shifting and folding operations)
%MrinaliniL
% AM.EN.U4ELC20022
x = [1 2 2];
h = [1 2 3 4];
n_1 = length(x);
n_2 = length(h);
N = max(n_1,n_2);
y = [x zeros(1,(N-n_1))];
g = [h zeros(1,(N-n_2))];
for i = 1:N
k = i;
for j = 1:n_2
H(i,j)=y(k)*g(j);
k = k-1;
if (k == 0)
k = N;
end
end
end
Y = zeros(1,N);
M=H' ;
for j = 1:N
for i = 1:n_2
Y(j)=M(i,j)+Y(j)
end
end
n = 0:3;
stem(n,Y)
xlabel("n")
ylabel("Y(n)")
title("Circular Convolution Mrinalini [AM.EN.U4ELC20022]")
(b) Compute the 4 - point circular convolution by frequency domain approach (The frequency-domain implementation uses the fact that the circular convolution of x[n] and h[n] is equivalent to the multiplication of their DFTs, i.e., Y [k] = X[k].H[k] and take the inverse DFT of Y [k]). (use dftmtx command)
%Mrinalini
% AM.EN.U4ELC20022
x = [1 2 2];
h = [1 2 3 4];
n_1 = length(x);
n_2 = length(h);
N = max(n_1,n_2);
L = 4;
y = [x zeros(1,(N-n_1))];
g = [h zeros(1,(N-n_2))];
b = dftmtx(4);
X_K = b*y'
H_K = b*g'
f = (X_K) .* (H_K);
C = [f]
d = conj(dftmtx(4));
Y_n = (1/L)*(d)*(f)
p = 0:3;
stem(p,Y_n)
xlabel("n")
ylabel("Y_n")
title("Circular Convolution using DFT and IDFT")
(c ) verify the results using cconv command.
%MRINALINI
%AM.EN.U4ELC20022
x = [1 2 2];
h = [1 2 3 4];
n_1 = length(x);
n_2 = length(h);
N = max(n_1,n_2);
a = [x zeros(1,(N-n_1))];
h = [h zeros(1,(N-n_2))];
y = cconv(a,h,4)
Find the discrete Fourier transform of the signal x[n]={2,2,2,2,1,1,1,1} using FFT algorithm and verify it by taking the inverse of the result. Plot the magnitude and frequency spectrum.
%MRINALINI
%AM.EN.U4ELC20022
x=[2,2,2,2,1,1,1,1];
y=fft(x)
magnitude=abs(y);
t=0:7;
subplot(2,1,1)
stem(t,magnitude)
title('FFT Magnitude')
xlabel('time')
ylabel('magnitude')
p=angle(y);
phase=180*p/pi;
subplot(2,1,2)
stem(t,phase)
title('Phase Plot')
xlabel('time')
ylabel('phase')
X=ifft(y)
2. Spectral Analysis of Sinusoidal Signals
This application discusses the use of the DFT in the determination of the spectral (frequency) contents of a continuous time signal. The effect of DFT length is examined in detail.
•
0 ≤ n ≤ N −1. Let the normalized frequencies of the two 16 length sequences be f1 = 0.22 and f2 = 0.34. Write a script le to compute the DFT of their sum x[n] by varying the values of the DFT length N from 16 to 128. Use the built-in command fft for computing DFT. Note: When N > 16, the M file fft automatically zero pads the sequence x.
%MRINALINI
%AM.EN.U4ELC20022
N = 16;
n = 0:15;
r = 0:31;
t = 0:63;
z = 0:127;
k = 0:N-1;
f_s = 1;
f_1 = 0.22;
f_2 = 0.34;
s_1 = (1/2)*sin(2*pi*f_1*n);
s_2 = sin(2*pi*f_2*n);
s = s_1 + s_2
a = fft(s,16);
b = fft(s,32);
e = fft(s,64);
f = fft(s,128);
c = abs(a);
d = abs(b);
x = abs(e);
y = abs(f);
subplot(4,1,1)
stem(n,c)
subplot(4,1,2)
stem(r,d)
subplot(4,1,3)
stem(t,x)
subplot(4,1,4)
stem(z,y)

%MRINALINI
%AM.EN.U4ELC20022
close all
n=0:15;
f_1=0.22;
f_2 = 0.63;
s_1 = (1/2)*sin(2*pi*f_1*n);
s_2 = sin(2*pi*f_2*n);
s = s_1 + s_2;
N = 16;
k=0:N-1;
w=(2*pi*k)/N;
w1=0:(2*pi)/1000:((2*pi)-((2*pi)/1000));
w2=0:(2*pi)/16:((2*pi)-((2*pi)/16));
X_dtft=mydtft(s,n,w);
figure
plot(w,abs(X_dtft));
hold on
X_dft=mydtft(s,n,w1);
stem(w1,abs(X_dft),'r')
xlabel('t')
ylabel('mag')
title('DTFT vs 16-point DFT')