%load the file
disp=load('pezzack.dat');
% Sampling frequency is 512
SF=50;

%calculate the FFT of the signal
FFTdisp = fft(disp);
%Calculate the square amplitude of the signal (power spectrum)
PS=FFTdisp.*conj(FFTdisp);
ndata=round(length(PS)/2);
i=1:ndata;
%the first PS number is the average
PSpositive(i)=PS(i+1);
%construct the frequency axis - up to half the sampling frequency
j=(SF/2)/ndata*i'

%calculate the overall sum of power and the % for each harmonic
sumPS=sum(PSpositive);
perc_powerPS=100*PSpositive/sumPS;

%Find the cumulative sum as a percentage for each harmonic
cumsumPS = cumsum(perc_powerPS)'