%% Lezione del 11 novembre 2008
%
% Calcolare la trasformata di Fourier discreta (DTFT) X(k)
% di un segnale causale sinusoidale discreto x(nT)
% di frequenza 200 Hz, lungo 2 s, a una frequenza di campionamento
% Fs = 1000 Hz [ricordare: omega = 2*pi*f/Fs = 2*pi*f*T]
%
% Di questa DTFT si osservi nell'intervallo [0,Fs/2-1]
% alla precisione dell'Hz
% 1) il modulo, abs(X)
% 2) il modulo in dB, 20*log10(abs(X))
% 3) la parte reale, real(X)
% 4) la parte immaginaria, imag(X)
% 5) il modulo, stavolta nell'intervallo [0,Fs-1] Hz
%
% Similmente, si osservi il modulo nell'intervallo [0,Fs/2-1] Hz della DTFT di
% un segnale causale sinusoidale x(nT), di frequenza 1300 Hz lungo 2 s
% alla frequenza di campionamento Fs = 1000 Hz. Motivare il risultato.
%
% Da ultimo, si osservi il modulo della DTFT di un
% segnale causale sinusoidale x(nT) esponenzialmente smorzato lungo 2s,
% in cui una sinusoide a frequenza 100 Hz ad ogni secondo si smorza di un
% fattore exp(-10), a una frequenza di campionamento Fs = 1000 Hz
%
%%

clear all

Fs = 1000;            % sampling frequency

n = [0:2*Fs-1];        % discrete-time domain (2s)

freq = 200;           % in Hz

%%%%%%%%%% First signal

x = sin(2*pi*freq/Fs*n);

%%%%%%%%%% Fourier tansform of first signal

K = zeros(0,Fs-1);  % empty frequency vector (memory allocation)

for k = 1:Fs    
    X(k) = x * exp(-j*2*pi*k*n'/Fs);
end

%%%%%%%%%% Figure 1

D = [0:Fs/2-1];        % frequency domain (in Hz)

figure(1)
plot(D,abs(X(1:Fs/2)));
xlabel('frequency (Hz)'), ylabel('magnitude'), title('Frequency magnitude');

pause

plot(D,20*log10(abs(X(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude (dB)');

pause

plot(D,real(X(1:Fs/2)));
xlabel('frequency (Hz)'), ylabel('magnitude'), title('Frequency magnitude (real part)');

pause

plot(D,imag(X(1:Fs/2)));
xlabel('frequency (Hz)'), ylabel('magnitude'), title('Frequency magnitude (imaginary part)');

%%%%%%%%%% Figure 2

D2 = [0:Fs-1];        % frequency domain (in Hz)

figure(2)
plot(D2,abs(X));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude');

%%%%%%%%%% Second signal

freq = 1300;

x = sin(2*pi*freq/Fs*n);

for k = 1:Fs    
    X(k) = x * exp(-j*2*pi*k*n'/Fs);
end

%%%%%%%%%% Figure 3

figure(3)
plot(D,abs(X(1:Fs/2)));
xlabel('frequency (Hz)'), ylabel('magnitude'), title('Frequency magnitude');

%%%%%%%%%% Third signal

freq = 100;
tau = 1/10;

x = sin(2*pi*freq/Fs*n) .* exp(-n/tau/Fs);

X = zeros(1,Fs);    % empty frequency vector (memory allocation)

for k = 1:Fs
    X(k) = x * exp(-j*2*pi*k*n'/Fs);
end

%%%%%%%%%% Figure 4

figure(4)
plot(n/Fs,x);
xlabel('time (s)'), ylabel('magnitude'), title('Discrete-time signal');

figure(5)
plot(D,abs(X(1:Fs/2)));
xlabel('frequency (Hz)'), ylabel('magnitude'), title('Frequency magnitude');
