Power Spectrum via fft computation

Good Morning,

I am having problems with a Matlab script.

I want to create a Power Spectrum starting from a discrete signal, using the theory's formulas.

The problem should be in the computational frame of fft, after the calculation of the autocorrelation.

clear

load C:\MATLAB6p1\work\mareggiate\file3\file3.dat


clc
close all

dt=input('Inserisci il dt della registrazione: ');
tmax=2*dt;    %teorema di shannon
a='La frequenza di campionamento è: %0.2f Hz';
fmax=1/tmax;        
sprintf(a,fmax)
disp('---------------------------------')

y=file3;      % modificare anche il save!
y=y';
media=mean(y);
npts=length(y);

%cacolo del tempo di registrazione
t_record=npts*dt
FS=1/t_record;
periodo_Max=(1/FS)/60;
sprintf('Il tempo di regitrazione è= %0.2f minuti',periodo_Max)
disp('---------------------------------')
for i=1:npts
    tempo(i)=(i*dt)/60;
end
input k

%plot della mareggiata
figure
plot(tempo,y),grid on, title('Mareggiata'),xlabel('lunghezza registrazione [minuti]'),ylabel('ampiezza [m]')
input k

b=npts; %+24;
taumax=fix(npts/2);

%autocorrelation, Bendat pg 381
som=0;
for k=1:taumax
    som=0;
    for i=1:taumax
        som=som+y(i)*(y(i+(k-1)));
    end
    autoc(k)=som/(taumax);     
end 
autoc=autoc';
figure
plot(autoc), grid on, title ('funzione di autocorrelazione'), xlabel('npts/2')
input k

%calcolo trasformata fourier
L=length(autoc)
n=2*nextpow2(L);
Xxx = fft(autoc,npts);        %two sided spectral density b
Xabs= (abs(Xxx).^2)/n;
Pxx = [Xabs(1);4*Xabs(2:(npts/2))];  %one sided spectral density, moltiplico per 2 pg 120 Bendat
df=fmax/(length(Pxx)-1)         %definizione delta f
f = 0 : df : fmax;              %1/2 1/b-1, parametri precedenti
f=f';

%calcolo del punto di picco
np=length(Pxx);
v=Pxx(2:np);
[a,c]=max(v);                    %posizione del picco
sprintf('La posizione del picco è= %0.1f punto',c)
disp('---------------------------------')
fp=c*df;                         %calcolo frequenza di picco
sprintf('La frequenza di picco è= %0.2f Hz',fp)
t_max=1/fp;                      %calcolo periodo massimo
sprintf('Il periodo massimo della mareggiata è= %0.2g secondi',t_max)
input k
%
%plot della PSD
figure
subplot(211)
plot(v), grid on, title('spettro di potenza mareggiata')
subplot(212)
plot(f,Pxx), grid on, title ('PSD Mareggiata'), xlabel('frequenze [Hz]'), ylabel ('potenza [Hz/m^2]')
input k

%funzione smooth
z=Pxx;
for k=1:100
    z=smoonew(z);
end
%azzeramento primi valori
for i=1:5
    z(i)=0;
end

%plot PSD e funzione smooth
input k
figure
plot(Pxx), grid on, title ('PSD file e funzione smooth')
hold on
plot(z,'r'), legend('PSD file','smooth'), ylabel('Ampiezza'), xlabel('npts/2')
input k
%
%calcolo dell'energia
somma=0;
for i=2:np                      %(b/2)
    somma=somma+z(i)*df;        %FS;
end 
energia_PSD=somma;
sprintf('Energia della PSD= %0.2f Hz/m^2',energia_PSD)

area=sqrt(energia_PSD);
Hs=4*sqrt(energia_PSD);
sprintf('La altezza significativa è= %0.2f m',Hs)
Hmax=1.86*Hs;
sprintf('La altezza del picco è= %0.1f m',Hmax)
input k




%save C:\MATLAB6p1\work\mareggiate\file3\PSD_File3X4.dat z  -ascii