function [z power] = DFTtest(y)

Y=fft(y);
n=length(Y);
t=0:(n-1);
figure(36)
plot(t,y,'b.-',t,y,'ro')
title('original function values')



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

%The following branches depending on n even or n odd. This is the even case.
if floor(n/2)==n/2

power = abs(Y(1:(1+n/2))).^2;

% Possibly this script exaggerates the relative power at 
% the Nyquist frequency Y(1+n/2) by a factor of 2. Check this.

freq = (1:(n/2))./n;
figure(38)
plot(freq(1:(n/2)),power(2:(1+n/2)))
ax = axis;
set (gca, 'xscale', 'log')
axis (ax)
xlabel('frequency, multiples of sample frequency')
ylabel('Power');



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

ylabel('Power');
xlabel('Period (Time/Cycle, Multiples of Sample Interval)');
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(40)
stem(freq,phase);
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'})
ax = axis;
set (gca, 'xscale', 'log')
axis (ax)
xlabel('frequency')
ylabel('phase angle')
title('phase versus frequency (when power > 1% max)','FontSize',20,'FontWeight','bold')
grid


%The following is the second branch, which depends on n even or n odd. 
%This is the odd case.

else

power = abs(Y(1:(1+(n-1)/2))).^2;
freq = (1:((n-1)/2))./n;
figure(38)
plot(freq(1:((n-1)/2)),power(2:(1+(n-1)/2)))
ax = axis;
set (gca, 'xscale', 'log')
axis (ax)
xlabel('frequency, multiples of sample frequency')
ylabel('Power');



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

ylabel('Power');
xlabel('Period (Time/Cycle, Multiples of Sample Interval)');
grid

hold on;
index=find( power(2:(1+(n-1)/2))==max(power(2:(1+(n-1)/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-1)/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(40)

stem(freq,phase);
yticks([-pi -pi/2 0 pi/2 pi])
yticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'})
ax = axis;
set (gca, 'xscale', 'log')
axis (ax)
xlabel('frequency')
ylabel('phase angle')
title('phase versus frequency (when power > 1% max)','FontSize',20,'FontWeight','bold')
grid



end

z=Y;