**********************
**********************
*LABORATORY SESSION 7*
**********************
**********************


%task-a
Fs = 1000; F = 250; A1 = 3; theta = 0; n = 0:511;
x = 3*cos(2*pi*n*F/Fs);
y = 3*sin(2*pi*n*F/Fs);

ftx = fft(x); %Fourier transform of signal x
N = length(x);
f = (0:(N-1))/N*Fs; %vector of frequencies

fty = fft(y);
N = length(y);
f = (0:(N-1))/N*Fs; %vector of frequencies

%plot real and imaginary part of sin
figure; 
subplot(2,1,1); plot(f,real(ftx)); %f tells Matlab to put frequencies on x axis
subplot(2,1,2); plot(f,imag(fty));
%plot real and imaginary part of cos
figure;
subplot(2,1,1); plot(f,real(ftx));
subplot(2,1,2); plot(f,imag(fty));



%task-b
%amplitude spectrum:
Fs = 8000; N = 1024;
A1 = 3; F1 = 220; theta1 = 0;
A2 = 0.5; F2 = 440; theta2 = 3*pi/4;

x1 = A1*sin(2*pi*F1*(0:N-1)/Fs + theta1);
x2 = A2*sin(2*pi*F2*(0:N-1)/Fs + theta2);
x3 = x1 + x2;
f = (0:(N-1))/N*fs;

ftx1 = fft(x1); %Fourier transform
ftx2 = fft(x2);
ftx3 = fft(x3);
%plot amplitude spectrum
figure; plot(abs(ftx1))
figure; plot(abs(ftx2))
figure; plot(abs(ftx3))


[x,fs] = wavread('aaa.wav');
N = length(x);
f = (0:(N-1))/N*fs;
ftx = fft(x);
figure; plot(f,abs(ftx));

[x,fs] = wavread('iii.wav');
N = length(x);
f = (0:(N-1))/N*fs;
ftx = fft(x);
figure; plot(f,abs(ftx));

[x,fs] = wavread('clarinet.wav');
N = length(x);
f = (0:(N-1))/N*fs;
ftx = fft(x);
figure; plot(f,abs(ftx));



%task-c
%DFT of rectangular signal
x1 = zeros(1,20); x1(1:5) = 1;
x2 = zeros(1,20); x2(16:20) = 1;
x3 = zeros(1,20); x3(2:6) = 1;
x4 = zeros(1,20); x4(6:10) = 1;

ftx1 = fft(x1,256); %tell Matlab to calculate Fourier transform using 256 points
ftx2 = fft(x2,256);
ftx3 = fft(x3,256);
ftx4 = fft(x4,256);

%amplitude spectrum
subplot(2,1,1)
hold on; 
plot(abs(ftx1),'b'); 
plot(abs(ftx2),'r');
plot(abs(ftx3),'m');
plot(abs(ftx4),'k');
hold off;
%phase spectrum
subplot(2,1,2)
hold on; 
plot(angle(ftx1),'b'); 
plot(angle(ftx2),'r');
plot(angle(ftx3),'m');
plot(angle(ftx4),'k');
hold off;


%task-d
%roots of the unit circle in the case N = 4.
%k-th root: W^k_N = e^{(-j 2 pi k)/N}; k = 0,1,...,N-1
N = 4; k = 0:N-1;
roots = exp(-j*2*pi*k/N);

%Calculate the four-point DFT of the following signal x(n) = {0, 1, 2, 3}, x(0) = 0, calculate also the amplitude and phase spectrum.
x = [0 1 2 3];
Xk = fft(x);  %6, -2 + 2i, -2, -2 - 2i

%calculation by "hand"
%use equation in lab7.pdf under Discrete Fourier Transform
k = 0;
X(k) = x(0)(cos(2 pi k 0/4) - jsin(2 pi k 0/4)) + 
     + x(1)(cos(2 pi k 1/4) - jsin(2 pi k 1/4)) + 
     + x(2)(cos(2 pi k 2/4) - jsin(2 pi k 2/4)) + 
     + x(3)(cos(2 pi k 3/4) - jsin(2 pi k 3/4)) = 
     = 0(1-j0) + 1(1-j0) + 2(1-j0) + 3(1-j0) = 
     = 6
     
k = 1;
X(k) = x(0)(cos(2 pi k 0/4) - jsin(2 pi k 0/4)) + 
     + x(1)(cos(2 pi k 1/4) - jsin(2 pi k 1/4)) + 
     + x(2)(cos(2 pi k 2/4) - jsin(2 pi k 2/4)) + 
     + x(3)(cos(2 pi K 3/4) - jsin(2 pi k 3/4)) = 
     = 0(1-j0) + 1(0-j1) + 2(-1-j0) + 3(0-j(-1)) = 
     = 0 -j -2 + 3j = -2 +2j
     
k = 2;
X(k) = x(0)(cos(2 pi k 0/4) - jsin(2 pi k 0/4)) + 
     + x(1)(cos(2 pi k 1/4) - jsin(2 pi k 1/4)) + 
     + x(2)(cos(2 pi k 2/4) - jsin(2 pi k 2/4)) + 
     + x(3)(cos(2 pi k 3/4) - jsin(2 pi k 3/4)) = 
     = 0(1-j0) + 1(-1-j0) + 2(1-j0) + 3(-1-j0) = 
     = 0 -1 + 2 -3 = -2     
     
k = 3;
X(k) = x(0)(cos(2 pi k 0/4) - jsin(2 pi k 0/4)) + 
     + x(1)(cos(2 pi k 1/4) - jsin(2 pi k 1/4)) + 
     + x(2)(cos(2 pi k 2/4) - jsin(2 pi k 2/4)) + 
     + x(3)(cos(2 pi k 3/4) - jsin(2 pi k 3/4)) = 
     = 0(1-j0) + 1(0-j(-1)) + 2(-1-j0) + 3(0-j1) = 
     = 0 +j -2 - 3j = -2 -2j     

X(k) = [6, -2+2j, -2, -2-2j]     


%amplitude spectrum
abs(fft(x,4)); %6, 2.8284, 2, 2.8284

%calculation by "hand" 
%|X(k)| = sqrt(XR(k)^2 + XI(k)^2)
|X(0)| = sqrt(6^2 + 0^2) = 6
|X(1)| = sqrt((-2)^2 + (2)^2) = sqrt(8) = 2.8284
|X(2)| = sqrt((-2)^2 + 0^2) = 2
|X(3)| = sqrt((-2)^2 + (-2)^2) = sqrt(8) = 2.8284

|X(k)| = [6, 2.8284, 2, 2.8284]     


%phase spectrum
angle(fft(x,4)); %0, 2.3562, 3.1416, -2.3562

%calculation by "hand" 
%theta(k) = arctan(XI(k)/XR(k))
theta(1) = arctan(0/6) = 0
theta(2) = arctan(2/(-2)) = 2.3562
theta(3) = arctan(0/(-2)) = 3.1416
theta(4) = arctan((-2)/(-2)) = -2.3562

theta(k) = [0, 2.3562, 3.1416, -2.3562]     



%Calculate the eight-point DFT of the following signal x(n) = {1, 1, 1, 1, 0, 0, 0, 0}, x(0) = 1, calculate also the amplitude and phase characteristic
x = [1, 1, 1, 1, 0, 0, 0, 0];
Xk = fft(x,8); %4, 1 - 2.4142i, 0, 1 - 0.4142i, 0, 1 + 0.4142i, 0, 1+ 2.4142i

%amplitude spectrum
abs(fft(x,8)); %4, 2.6131, 0, 1.0824, 0, 1.0824, 0, 2.6131

%phase spectrum
angle(fft(x,8)); %0, -1.1781, 0, -0.3927, 0,  0.3927, 0, 1.1781









