从网上下载了frft的matlab程序，没有进行归一化处理，在仿真的过程中角度是对的，的，所以斜率是对的，但是初始频率一直是错的，请问是我的方法出现了什么问题吗？
附程序段：
ph=0; %正弦信号的初始相位
B=60000;
N=12800;
fs=1e5;
alpha=468750;
t = (0:N-1)/fs; %n倍的采样周期
dt=1/fs;
y = exp(j*2*pi*(-B/2*t+alpha/2*t.*t));
p=-1:0.01:1;
u=1:1:length(t);
p3=zeros(length(u),length(p));
for k=1:1:length(p)
p3=frft(y,p(k));
Pp3(:,k)=p3.*conj(p3)/length(y);
end
surf(p,u,Pp3);
[CX,CI]=max(max(Pp3));
p(CI)
计算出来p=-0.66，用公式 fm=-cotθf0*f0/N 斜率对的上
但是出来的u=6362 用公式f0=ucscθf0/N 初始频率对应不上
请问我公式错了还是怎么回事
FrFt
function Faf = frft(f, a)
% The fast Fractional Fourier Transform
% input: f = samples of the signal
% a = fractional power
% output: Faf = fast Fractional Fourier transform
error(nargchk(2, 2, nargin));
f = f(:);
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;%rem（）取余数；fix（）取整数部分;总体是右边的一半数移到左边；
sN = sqrt(N);
a = mod(a,4);
% do special cases
if (a==0), Faf = f; return; end;
if (a==2), Faf = flipud(f); return; end;%flipud turn oppsite
if (a==1), Faf(shft,1) = fft(f(shft))/sN; return; end
if (a==3), Faf(shft,1) = ifft(f(shft))*sN; return; end
% reduce to interval 0.5 < a < 1.5
if (a>2.0), a = a-2; f = flipud(f); end
if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end
if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end
% the general case for 0.5 < a < 1.5
alpha = a*pi/2;
tana2 = tan(alpha/2);
sina = sin(alpha);
f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];%increase sampling rate
% chirp premultiplication
chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);
f = chrp.*f;
% chirp convolution
c = pi/N/sina/4;
Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.^2),f);
Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);
% chirp post multiplication
Faf = chrp.*Faf;
% normalizing constant
Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);
function xint=interp(x)
% sinc interpolation
N = length(x);
y = zeros(2*N-1,1);
y(1:2:2*N-1) = x;
xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2));
xint = xint(2*N-2:end-2*N+3);
function z = fconv(x,y)
% convolution by fft
N = length([x(:);y(:)])-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);
