% E03_1.m                10mar17   19mar20  05feb21
% DTFT Computation & properties
clear all;

%% 03.1-22   Principal value of the phase function and 2*pi discontinuities
%    due to the inverse tangent of (imag/real) 
%    (see Matlab doc for atan2)
[X,Y] = meshgrid(-10:0.2:10,-10:0.2:10);
P = atan2(Y,X);
surf(X,Y,P);  % a 2*pi wide discontinuity exists at Y=0 for all X<0.
xlabel('real'); ylabel('imag'); zlabel('tan^{-1}(imag/real)')
axis('tight'); view(30,45);
return

%% 03.1-26    DTFT of [alpha^n * mu(n)]
alpha = .9;  % test alpha = [.5, .9]
ome = linspace(0,pi,1000); % dense sampling of continuous ang.freq.

X = 1 ./ (1 - alpha * exp(-j*ome)); % exact transform

N = 20; % finite sequence [truncated impulse response]
X_hat = 0;
for n = 0:N-1,
    X_hat = X_hat + (alpha*exp(-j*ome)).^n;
end;

plot(ome, abs(X), 'b', ome, abs(X_hat), 'r');

%%  03.1-44   Gibbs
clear all;
omeC = 0.3*pi;
M = 1000;  % dense sampling of continuous ang.freq. in [0 pi]
ome = linspace(0,pi,M);
for N = [5 10 20 40 100] % progressively longer summations
H = zeros(1,M);
for n = -N:N
    HH = sin(omeC*n)/pi./n .* exp(-j*ome*n);
    if n==0, HH = ones(1,M)/pi; end   % sin(0)/0 --> 1
    H = H + HH;
end
plot(real(H)); axis([0 M-1 -0.2 1.2]); pause()
hold on;
end
hold off

%%  03.1-56 and 03.2.14    Frequency-Shifting Property of DTFT
num = [1  3  5  7  9  11  13  15  17];  
h1 = freqz(num, 1, ome); 
n = 0:length(num)-1;
h2 = freqz(exp(j*.3*pi*n).*num, 1, ome);   % num modulates a complex exp
subplot(2,2,1);  plot(ome/pi,abs(h1));grid
title('Magnitude Spectrum of Original Sequence')
subplot(2,2,2);  plot(ome/pi,abs(h2));grid
title('Magnitude Spectrum of Frequency-Shifted Sequence') % shifted by .3 pi
subplot(2,2,3);  plot(ome/pi,angle(h1)/pi);grid
title('Phase Spectrum of Original Sequence')
subplot(2,2,4);  plot(ome/pi,angle(h2)/pi);grid
title('Phase Spectrum of Frequency-Shifted Sequence')

%%  03.2-15   Example
L = 256;  % Number of frequency points (length of the DTFT)
num = [.008  -.033  .05  -.033  .008];    % Numerator coefficients
den = [1  2.37  2.7  1.6  .41];   % Denominator coefficients
% num = 1;  den = [1 -.99];  % exp. weighted accumulator
ome = linspace(0,pi,L);
h = freqz(num, den, ome);
subplot(2,2,1); plot(ome/pi,real(h)); grid; xlabel('\omega/\pi'); ylabel('Real')
subplot(2,2,2);  plot(ome/pi,imag(h)); grid; xlabel('\omega/\pi'); ylabel('Imag.')
subplot(2,2,3);  plot(ome/pi,abs(h)); grid; xlabel('\omega/\pi'); ylabel('Magnitude')
subplot(2,2,4);  plot(ome/pi,angle(h)/pi);grid; xlabel('\omega/\pi'); ylabel('Phase, radians /pi')

%%  03.2-20   Unwrapping (& Time-Shifting Property of DTFT)
num = [1  3  5  7  9  11  13  15  17];  %(not the same as in the slides)
ome = linspace(0,pi,100);  %  note errors if few freq. samples are used
h1 = freqz(num, 1, ome);  % note den = 1; 
h2 = freqz([zeros(1,7) num], 1, ome);   % num is time-shifted by 7 samples
subplot(3,2,1);  plot(ome/pi,abs(h1));grid
title('Magn. Spectrum of Original Seq.')
subplot(3,2,2);  plot(ome/pi,abs(h2));grid
title('Magn. Spectrum of Time-Shifted Seq.')
subplot(3,2,3);  plot(ome/pi,angle(h1)/pi);grid
title('Phase Spectrum of Original Seq.')
subplot(3,2,4);  plot(ome/pi,angle(h2)/pi);grid
title('Phase Spectrum of Time-Shifted Seq.')  % original phase + omega*7
subplot(3,2,5);  plot(ome/pi,unwrap(angle(h1))/pi);grid
title('Unwrapped Phase of Original Seq.')
subplot(3,2,6);  plot(ome/pi,unwrap(angle(h2))/pi);grid
title('Unwrapped Phase of Time-shifted Seq.')

