大家帮我看看这个程序做了个啥，谢谢

xiuhuwang

%This program is for observing power spectral density
format long
Nfft=2^10;
omega=0.1;
A=0.1;
gamma=0.1;
PT=2.0*pi/omega;
M=Nfft; %the mumber of discrete points in esch period
dt=PT/M;
N=20;
NM=50;
D=0.2;
x0=1.0;
for i=1:1:N*M
    R=rand(2,1);
    gw=(-4*D*dt*dt*log(R(1,1)))^0.5*cos(2.0*pi*R(2,1));
    t=(i-1)*dt;
    x=x0+(x0*(1+A*cos(omega*t))-x0^3.0)*dt+gw;
    x0=x;
end

Fs=1.0/(8*dt);%define the sampling frequency
P=zeros(Nfft/2+1,1);%initialize
for kk=1:1:NM
    disp('1=,kk=')
    kk
    XX=zeros(Nfft,1);
    k=0;
    for i=1:1:8*M;
        R=rand(2,1);
        gw=(-4*D*dt*log(R(1,1)))^0.5*cos(2.0*pi*R(2,1));
        t=(i-1)*dt;
        x=x0+(x0*(1+A*cos(omega*t))-x0^3.0)*dt+gw;
        if (mod(i,8)==0)
            k=k+1;
            TT(k,1)=t;
            XX(k,1)=x;
        end
        x0=x;
    end
    [Pxx,F]=psd(XX,Nfft,Fs);
    P=P+Pxx;
end
P=P./NM;
plot(F(1:40),P(1:40))