请教频域信号转化到时域

lxw

下面的程序是关于陈果《铁路轨道不平顺随机过程的数值模拟》论文中将信号从频域转化到时域的代码，但是出来的结果和他在论文中给出的差距挺大的，主要是开始和结尾的值太大了，请高手帮忙看一下问题出在哪里，另外这个程序也是这个论坛别人发过的，当时讨论也没得出什么结论，我感觉也没什么问题，谢谢。
% 频域离散点采样
%nextpow2() 取最接近的较大2次幂
clear;
clc;
Av=[1.2107,1.0181,0.6816,0.5376,0.2095,0.0339];
Aa=[3.3634,1.2107,0.4128,0.3027,0.0762,0.0339];
Ws=[0.6046,0.9308,0.8520,1.1312,0.8209,0.4380];
Wc=[0.8245,0.8245,0.8245,0.8245,0.8245,0.8245];
k=0.25;
rand('state',0);%%把正态随机数发生器置0
%yl=1.524;yu=304.8;
yl=0.5;yu=50;
pr1={'V='};
in1=inputdlg(pr1,'线路参数模拟信息',1,{'100'});
V=str2num(in1{1});
pr2={'选择路面等级i=','输入不平顺种类(1-高低不平顺，2－方向不平顺，3－水平及轨距不平顺j='};
in2=inputdlg(pr2,'线路参数模拟信息',1,{'6','1'});
i=str2num(in2{1});
j=str2num(in2{2});
v=V/3.6;
switch j
case 1
s=@(f)(k*Av(i)*Wc(i)^2./((f.*2*pi/v).^2.*((f*2*pi/v).^2+Wc(i)^2)))/v;%%s=k*Av(i)*Wc(i)^2/(w^2*(w^2+Wc(i)^2));
case 2
s=@(f)(k*Aa(i)*Wc(i)^2./((f*2*pi/v).^2.*((f*2*pi/v).^2+Wc(i)^2)))/v;
case 3
s=@(f)(4*k*Av(i)*Wc(i)^2./((f*2*pi/v).^2+Ws(i)^2)./((f*2*pi/v).^2+Wc(i)^2))/v;
end
fu=v/yl;fl=v/yu;
co=1;%判断点的个数和Nr/2的关系
while(co>0)
Ts=inputdlg('模拟时间Ts=','模拟的时间长度',1,{'10'});
Ts=str2num(Ts{1});
dt=0.0001;
N=Ts/dt;
Nr=2^(nextpow2(N));
Xm=zeros(1,Nr/2+1);
df=1/(Nr*dt);
n0=ceil(fl/df);
nf=ceil(fu/df);
if nf<=Nr/2
S=feval(s,(n0+1:nf)*df);
Xm(n0+1:nf)=Nr*sqrt(S*df);
break;
% N0+1～N0+Nf数值非零，其余为零。
end
end
n=nf-n0;
q=2*pi*rand(1,Nr/2+1);
X=Xm.*q;
temp=fliplr(X);
X=[X,temp(2:end)];
x=ifft(X,Nr);
t=(1:Nr)*dt;
%%-->频域法模拟结果
figure(1);
subplot(2,1,1);
plot(t,x);
grid on;
Fs=1/dt;
fn=(1:Nr/2)/Nr*Fs;
Xnn=fft(x,Nr);
Xn=Xnn(1:Nr/2);
Sn=(abs(Xn)*2/Nr).^2;
figure(1);
subplot(2,1,2);
loglog(fn,Sn,'* g',fn,s(fn),'r');
legend('模拟值','解析值');
xlim([0,100]);
ylim([1e-8,1]);

linlin820

程序估计是没有什么问题，不过取值是不是有所不同？
现在国内文章都这样，很多取值都不给全，让人无法再现