%% Lezione del 11 novembre 2009
%
% Si definisca un segnale causale discreto x lungo 1 s costituito dalla
% sovrapposizione di tre componenti sinusoidali, rispettivamente poste a
% 100, 4400 e 8800 Hz, a una frequenza di campionamento Fs = 44100 Hz
%
% 1) si osservi la risposta in frequenza di x su 44100 "bin"
% [sugg.: adoperare la primitiva fft(x,n), vedi anche "help fft"]
%
% Ora si consideri il segnale causale discreto (media mobile)
% h(nT) = [.1, .1, .1, .1, .1, .1, .1, .1, .1, .1] con T=1/Fs
%
% 2) si osservi il modulo della risposta in frequenza H(omega) su 44100 "bin"
%
% 3) si calcoli y(n) = h * x (n) adoperando la definizione di convoluzione
% discreta
%
% 4) si osservi il modulo della risposta in frequenza y(omega) su 44100 "bin"
%
% 5) si ascoltino x e y e si motivino le differenze percepite nei due segnali
% in termini di differenze spettrali
%
%%

clear all

Fs = 44100;           % sampling frequency

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

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

%%%%%%%%%% Tones

freq1 = 100;
x1 = sin(2*pi*freq1/Fs*n);

freq2 = 4400;
x2 = sin(2*pi*freq2/Fs*n);

freq3 = 8800;
x3 = sin(2*pi*freq3/Fs*n);

x = x1+x2+x3;

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

X = fft(x,Fs);    % FFT of x

figure(1)
semilogx(0:Fs/2-1,20*log10(abs(X(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude of x(n)');

%%%%%%%%%% System

h = [.1, .1, .1, .1, .1, .1, .1, .1, .1, .1];

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

H = fft(h,Fs);    % FFT of x

figure(2)
semilogx(0:Fs/2-1,20*log10(abs(H(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude of h(n)');

%%%%%%%%%% Output

y = zeros(1,Fs);      % discrete-time domain (1s)

for k = 1:Fs
    l = min(k,9);
    y(k) = h(1:l) * x(k:-1:k-l+1)';
end

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

Y = fft(y,Fs);
Z = H.*X;    % FFT of x

figure(3)
semilogx(0:Fs/2-1, 20*log10(abs(Y(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude of y(n)');
figure(4)
semilogx(0:Fs/2-1, 20*log10(abs(Z(1:Fs/2))));
xlabel('frequency (Hz)'), ylabel('magnitude (dB)'), title('Frequency magnitude of y(n)');

wavwrite(x/max(abs(x)*1.01),Fs,'x.wav');
wavwrite(y/max(abs(x)*1.01),Fs,'y.wav');
