% Here we solve Chapra 4th Edition Problem 16.12

%******************************************
function w = f16_12(t)
w=2.5.*cos(10.*t)-3.1.*sin(3.*t);
% main periods 2*pi/10 = 0.6283 and 2*pi/3 = 2.0944
%******************************************
function w = randomizer(y,mn,s)
%takes a matrix y and adds a normally distributed random value to each entry. The added vales have mean mn and st. dev. s.
w=y+mn+s*randn(size(y));
%******************************************




t=linspace(0,2*pi,64);
y=f16_12(t);
yrand=randomizer(y,0,1);
plot(t,y,t,yrand,'*');
[z power] = DFTtest(yrand);



% The main period is identified in Figure 39 as 21.3333 sample periods.
% ach sample period is (2*pi)/64, and 21.3333*(2*pi)/64=2.0944.

% The second most-important period is at 6.4 sample intervals, 
% and 6.4*(2*pi)/64 = 0.6283, as it should be.

Y=fft(y);

n=length(Y);

figure(12)
plot(t,y,'b.-')
title('original function values')

%Y(1) is the sum of the data. Y(2) is the coefficient of the longest possible period (in this case 2*pi) visible in the data
%Y(1+n/2) is the coeff. of the Nyquist frequency term, (frequency 16/pi, wavelength pi/16)
%Y(64-j) is the conjugate of Y(j+2).

figure(13)
plot(Y(2:n),'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');

figure(14)
power = abs(Y(1:(1+n/2))).^2;
freq = (1:(n/2))*(1/(2*pi));
figure(14)
plot(freq(1:(n/2)),power(2:(1+n/2)))
xlabel('frequency')
ylabel('Power');

figure(15)
period=1./freq;
plot(period(1:(n/2)),power(2:(1+n/2)));

ylabel('Power');
xlabel('Period (Time/Cycle)');

figure(16)
plot(period(1:(n/2)),power(2:(1+n/2)));

ylabel('Power');
xlabel('Period (Time/Cycle)');
grid

hold on;
index=find( power(2:(1+n/2))==max(power(2:(1+n/2))) );
mainPeriodStr=num2str(period(index));
plot(period(index),power(index+1),'r.', 'MarkerSize',25);
text(period(index)+1,power(index+1),[' Period = ',mainPeriodStr]);
hold off;

h=(10e-2)*max(abs(Y(2:(n/2))));
w=1;
while w<length(freq)+1
if    norm(Y(w))<h   phase(w)=0;
elseif  (imag(Y(w)) >0)  phase(w)=acos( real(Y(w))/norm(Y(w)));
else  phase(w)=-acos( real(Y(w))/norm(Y(w)));
end
w=w+1;
end
figure(17)
stem(freq,phase);
xlabel('frequency')
ylabel('phase angle')
title('phase versus frequency (when power > 1% max)','FontSize',20,'FontWeight','bold')
grid
%******************************************