Exploring DTFT concept using MATLAB

I am trying to understand the concept of DTFT using MATLAB. So I am doing the following experiments:

1. plot X1[n]=COS(2*pi*201/1024*n) where n=0 to1023
2. plot X2[n]=rect(N) where N=1024 i.e. X2[n]=1 for n=0 to 1023; else X2[n]=0
3. plot multiplication of both i.e. X1(n).*X2(n)
4. plot magnitude plot of Y1(w)=DTFT of { X1(n)*X2(n)} for -pi<w<pi

In figure 7, I expect the peak of sinc function (due to convolution of cosine and rectangular pulse) exactly at the location of the tone freq, but it does not. What may be the issue ???

My code is as below:

clear all; clc; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nfft = 2^12;          %  FFT size

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Generate single tone sine wave %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tone_num=Nfft/5;
tone_bin =  (2*floor(tone_num/2)+1)+0.5;  % prime wrt Nfft 
n = [0:Nfft-1];
xsin = cos(2*pi*tone_bin/Nfft*n);

%%%%%%%%%%% Plot single tone sine wave %%%%%%%%%%%%%%%%%
figure(1); plot(n,xsin);

%%%%%%%%%%% Plot DTFT of single tone sine wave %%%%%%%%%
[X1 , W1] = dtft ( xsin , Nfft) ;
figure(2);
plot( W1/2/pi , abs (X1) ) ; grid , title( ' MAGNITUDE RESPONSE ')
xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Generate rectangular pulse %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Nsampl = 200;     % number of samples
 width=100;        % width of the rectangular pulse 
 sampl_offest=40;  % offset from first sampple from where pulse will start
 xrect = [zeros(1,sampl_offest) ones(1,width) zeros(1,Nsampl-width-sampl_offest)];

%%%%%%%%%%%%%% PLOT rectangular pulse %%%%%%%%%%%%%%%%%
figure(3);stem(xrect);

%%%%%%%%%%% Plot DTFT of rectangular pulse %%%%%%%%%%%%
[X2 , W2] = dtft ( xrect , Nfft) ;
figure(4);
plot( W2/2/pi , abs (X2) ) ; grid , title( ' MAGNITUDE RESPONSE ')
xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )

●

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Generate product of sine & rectangular pulse %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xpdt=xsin(1:Nsampl).*xrect;

%%%%%% Plot product of sine & rectangular pulse %%%%%%%
figure(5); plot(xpdt);

%%% Plot DTFT of product of sine & rectangular pulse %%
[X3 , W3] = dtft ( xpdt , Nfft) ;
figure(6);
plot( W3/2/pi , abs (X3) ) ; grid , title( ' MAGNITUDE RESPONSE ')
xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )
ylim([-5 max(ylim)+5])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(7);
plot( W3/2/pi , abs (X3) ) ; grid , title( ' MAGNITUDE RESPONSE ')
xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )

k_tone=tone_bin/Nfft; % normalized freq. corresponding to sine tone
xline(k_tone,':r',{'sine tone'}); % plot normalized freq. corresponding to sine tone

wk=1/(Nfft);          % freq. bin width

ylim([-5 max(ylim)+5])
xlim([k_tone-30*wk k_tone+30*wk]) % zoom around sine tone freq
%xlim([0 0.5])

%%% generate xticks at every freq. bin w(k) = 2*pi/Nfft %%%%%%%%%%
 xwtk=[-0.5  0  0.5];
 i=1;
 while i*wk < 0.5
     xwtk = [xwtk i*wk];
    i=i+1;
 end
xwtk=sort(xwtk); 
  xticks(xwtk);

  %xline([-k_tone k_tone],':','c','sin_tone');
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%      END           %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% DEFINE FUNCTION dtft subroutine %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [H, W] = dtft(h, N)
%DTFT   calculate DTFT at N equally spaced frequencies
%----
%   Usage:   [H, W] = dtft(h, N)
%
%      h : finite-length input vector, whose length is L
%      N : number of frequencies for evaluation over [-pi,pi)
%              ==> constraint: N >= L 
%      H : DTFT values (complex)
%      W : (2nd output) vector of freqs where DTFT is computed

%---------------------------------------------------------------
% copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim,
% T.W. Parks, R.W. Schafer, & H.W. Schussler.  For use with the book
% "Computer-Based Exercises for Signal Processing Using MATLAB"
% (Prentice-Hall, 1994).
%---------------------------------------------------------------

N = fix(N);
L = length(h);  h = h(:);  %<-- for vectors ONLY !!!
if( N < L )
   error('DTFT: # data samples cannot exceed # freq samples')
end
W = (2*pi/N) * [ 0:(N-1) ]';
mid = ceil(N/2) + 1;
W(mid:N) = W(mid:N) - 2*pi;   % <--- move [pi,2pi) to [-pi,0)
W = fftshift(W);
H = fftshift( fft( h, N ) );  %<--- move negative freq components
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Paul]

Hi Pankaj,

Not an answer to the question, but ... why doesn't the code follow the instructions in the problem? Neither xsin nor xrect are as defined in steps 1 and 2.

Also, were you required to write your own DTFT function?

[Paul]

Running the code as is and only showing the final plot ....

Is the peak not exactly where it's expected?

clear all; clc; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nfft = 2^12;          %  FFT size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Generate single tone sine wave %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tone_num=Nfft/5;
tone_bin =  (2*floor(tone_num/2)+1)+0.5;  % prime wrt Nfft 
n = [0:Nfft-1];
xsin = cos(2*pi*tone_bin/Nfft*n);
%%%%%%%%%%% Plot single tone sine wave %%%%%%%%%%%%%%%%%
%figure(1); plot(n,xsin);

%%%%%%%%%%% Plot DTFT of single tone sine wave %%%%%%%%%
[X1 , W1] = dtft ( xsin , Nfft) ;
%figure(2);
%plot( W1/2/pi , abs (X1) ) ; grid , title( ' MAGNITUDE RESPONSE ')
%xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Generate rectangular pulse %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Nsampl = 200;     % number of samples
 width=100;        % width of the rectangular pulse 
 sampl_offest=40;  % offset from first sampple from where pulse will start
 xrect = [zeros(1,sampl_offest) ones(1,width) zeros(1,Nsampl-width-sampl_offest)];
%%%%%%%%%%%%%% PLOT rectangular pulse %%%%%%%%%%%%%%%%%
%figure(3);stem(xrect);

%%%%%%%%%%% Plot DTFT of rectangular pulse %%%%%%%%%%%%
[X2 , W2] = dtft ( xrect , Nfft) ;
%figure(4);
%plot( W2/2/pi , abs (X2) ) ; grid , title( ' MAGNITUDE RESPONSE ')
%xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Generate product of sine & rectangular pulse %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xpdt=xsin(1:Nsampl).*xrect;
%%%%%% Plot product of sine & rectangular pulse %%%%%%%
%figure(5); plot(xpdt);

%%% Plot DTFT of product of sine & rectangular pulse %%
[X3 , W3] = dtft ( xpdt , Nfft) ;
%figure(6);
%plot( W3/2/pi , abs (X3) ) ; grid , title( ' MAGNITUDE RESPONSE ')
%xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )
%ylim([-5 max(ylim)+5])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7);
plot( W3/2/pi , abs (X3) ) ; grid , title( ' MAGNITUDE RESPONSE ')
xlabel( ' NORMALIZED FREQUENCY \omega/(2\pi)' ) , ylabel( ' I H(\omega) I ' )
k_tone=tone_bin/Nfft; % normalized freq. corresponding to sine tone
xline(k_tone,':r',{'sine tone'}); % plot normalized freq. corresponding to sine tone
wk=1/(Nfft);          % freq. bin width
ylim([-5 max(ylim)+5])
xlim([k_tone-30*wk k_tone+30*wk]) % zoom around sine tone freq

%xlim([0 0.5])
%%% generate xticks at every freq. bin w(k) = 2*pi/Nfft %%%%%%%%%%
%{
 xwtk=[-0.5  0  0.5];
 i=1;
 while i*wk < 0.5
     xwtk = [xwtk i*wk];
    i=i+1;
 end
xwtk=sort(xwtk); 
xticks(xwtk);
%}
  %xline([-k_tone k_tone],':','c','sin_tone');
  

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%      END           %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% DEFINE FUNCTION dtft subroutine %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [H, W] = dtft(h, N)
%DTFT   calculate DTFT at N equally spaced frequencies
%----
%   Usage:   [H, W] = dtft(h, N)
%
%      h : finite-length input vector, whose length is L
%      N : number of frequencies for evaluation over [-pi,pi)
%              ==> constraint: N >= L 
%      H : DTFT values (complex)
%      W : (2nd output) vector of freqs where DTFT is computed
%---------------------------------------------------------------
% copyright 1994, by C.S. Burrus, J.H. McClellan, A.V. Oppenheim,
% T.W. Parks, R.W. Schafer, & H.W. Schussler.  For use with the book
% "Computer-Based Exercises for Signal Processing Using MATLAB"
% (Prentice-Hall, 1994).
%---------------------------------------------------------------
N = fix(N);
L = length(h);  h = h(:);  %<-- for vectors ONLY !!!
if( N < L )
   error('DTFT: # data samples cannot exceed # freq samples')
end
W = (2*pi/N) * [ 0:(N-1) ]';
mid = ceil(N/2) + 1;
W(mid:N) = W(mid:N) - 2*pi;   % <--- move [pi,2pi) to [-pi,0)
W = fftshift(W);
H = fftshift( fft( h, N ) );  %<--- move negative freq components
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Pankaj]

Actually I am trying to understand and emulate the theory of spectral leakage when the signal tone does not fall in a freq. bin.

The flow that I want to emulate in MATLAB is as given below

[Paul]

Again, the code doesn't follow from the question. For starters, it's apparent from the question and the graphic that N = 1024. Hence the code should start with something like:

N = 1024;
n = 0:N-1;
w0 = 201/N*2*pi; % case 1
xsin = cos(w0*n);

But the code in the question uses

Nfft = 2^12

Nfft = 4096

to define the length of xsin.