% E03_2.m              20mar18   22mar22
%  DFT and IDFT Computation
clear all;
return

%% note transformed sequence is NOT exactly symmetrical
%
L = 17;  % try L even & L odd
 x = (rand(1,L)-.5)*10;  % generate random data
% x = 10 * [1:L] .* exp(-[1:L]/3); % or determin. data
L2 = floor(L/2);
L_even = 2*(L2-L/2)+1; % =1 if L even, =0 if L odd

x, sum(x)
X = fft(x);  % calculate DFT
Xm = abs(X); Xp = angle(X);  % magnitude and phase
[Xm; Xp]'                    % (could use also real and imag)

     % keep non-redundant transform terms only (L real numbers)
Xm1 = Xm(1:L2+1);  Xp1 = Xp(1:L2+1);
[Xm1; Xp1]'     % note if L even, nonredund. terms are L/2+1 !

     % lossy compression of the original sequence
     %%%%%%%%%%% try also with larger L, different trunc.length,
% Xm1(end-1:end)=0;  % discard hf components

     % recreate "symmetric" terms
 Xm2 = [Xm1, fliplr(Xm1(2:end-L_even))];
 Xp2 = [Xp1, -fliplr(Xp1(2:end-L_even))];
 [Xm2; Xp2]'

x_hat = ifft(Xm2 .* exp(j*Xp2))
mse = mean((x-x_hat).^2)


%%  DFT of important sequences
N = 16;            % length of the sequence
M = N;  % = N, or e.g. 128;  % length of the DFT  (>= N)
 u = [1 zeros(1,N-1)];          %  impulse
% u = [zeros(1,4) 1 zeros(1,N-5)];  %  delayed impulse
 r = 3; n = 0:N-1;  u = cos(2*pi*r*n/N);  %  cosine --> 2 pulses in the DFT iff r=integer (and M=N)
U = fft(u,M);   % Compute its M-point DFT
figure(1);
subplot(3,1,1);   n = 0:1:N-1;   stem(n,u);
title('Original time-domain sequence'); xlabel('Time index n'); ylabel('Amplitude')
subplot(3,1,2);   k = 0:1:M-1;   stem(k,abs(U));
title('Magnitude of the DFT samples'); xlabel('Frequency index k'); ylabel('Magnitude')
subplot(3,1,3);   stem(k,angle(U));
title('Phase of the DFT samples'); xlabel('Frequency index k'); ylabel('Phase')

uu = ifft(U,M);
figure(2);
subplot(2,1,1);   n = 0:1:M-1;   stem(n,real(uu))
title('Real part of the IDFT'); xlabel('Time index n'); ylabel('Amplitude')
subplot(2,1,2);   stem(n,imag(uu))
title('Imaginary part of the IDFT'); xlabel('Time index n'); ylabel('Amplitude')


%% time-domain aliasing if sampling rate is too low (slide 03.2-52+)
clear all;
x = [1 2 3 4 5 6];
M = length(x);
N = M;  % = M, or e.g. M-2, or M+2;
x = [x  zeros(1,N-M)];  % in case N>M
X = zeros(1,N);
for k = 0:N-1; % freq. sampling is insufficient if N<M
for n = 0:M-1;
    X(k+1) = X(k+1) + x(n+1).*exp(-j*2*pi*k*n/N);
end, end
x
x_hat = real(ifft(X,N))

XX = fft(x,N); % note Matlab picks up only the first N data
xx_hat = ifft(XX,N)  % ...and the result is 'correct'

%% ...and DTFT approximation if N >> M

x = [0 1 2 3 4 5 4 3 2 1];
M = length(x);
% r = 1.6; n = 0:M-1;  x = cos(2*pi*r*n/M);
for N = M:10:M*10;  
  x = [x  zeros(1,N-M)];  % in case N>M
  X = zeros(1,N);
  for k = 0:N-1;
  for n = 0:M-1;
    X(k+1) = X(k+1) + x(n+1).*exp(-j*2*pi*k*n/N);
    end, end  
    % NOTE simple X = fft(x,N) command exists instead
  stem(abs(X)); pause();
  end;
hold off


