% powerspec.m
% Periodogram Example

actual_f_bin = 1 + 13/2
   % This will be the exact (non-integer) bin for the frequency.
   % The 1+ is because Octave/Matlab arrays are unit offset.
   % The second term is a frequency in in cycles per domain.
N = 200; % Number of bins
x = (0:N-1)*2*pi*(actual_f_bin - 1)/(N-1);
y = cos(x);
y2 =  hanning(N)' .* y;

figure(1)
hold off
plot(1:N, y, 'r', 1:N, y2, 'k');
xlabel('x (spatial bin)')
ylabel('signal y(x)')
legend('raw data','windowed')
%axis([1,N])
print('powerspec1.pdf','-dpdfwrite')

yf = fft(y);
y2f = fft(y2);
ps = abs(yf).^2;
ps2 = abs(y2f).^2;

figure(2)
hold off
normalization = sum(y.^2)/sum(y2.^2)
plot(ps,'r')
hold on
plot(normalization*ps2,'k');
xlabel('k (spectral bin)')
ylabel('power spectrum')
legend('Raw power spectrum','PS of windowed data')
axis([1,N,0,max(ps)])
print('powerspec2.pdf','-dpdfwrite')

figure(3)
hold off
plot(ps(1:20),'r')
hold on
plot(normalization*ps2(1:20),'k');
hold on
xlabel('k (spectral bin)')
ylabel('power spectrum')
plot(actual_f_bin*[1,1], [0,1.1*max(ps2*normalization)],'b-');
centroid  = sum(  ps(1:20).*(1:20) )./sum( ps(1:20))
centroid2 = sum( ps2(1:20).*(1:20) )./sum(ps2(1:20))
plot(centroid*[1,1], [0,1.1*max(ps2*normalization)],'r--');
plot(centroid2*[1,1], [0,1.1*max(ps2*normalization)],'k--');
legend('Raw power spectrum','PS of windowed data','Exact signal frequency','Raw PS centroid','Windowed PS centroid')
axis([1,20,0,max(ps)])
print('powerspec3.pdf','-dpdfwrite')
