最近使用王兆华教授的全相位fft去计算未知频率信号的频率时，发现得出的结果频率值与期望值相差有点大，是不是程序哪里出错了？还是apfft使用的方式错了。
如您已解决请联系我修改酬金至100~500(看最后的误差多少，<0.1s时100,<0.01时200，<0.001时500)
附上预期频率和计算结果的误差图：
未知信号的频率 求解出的频率 误差
100 100.2298641 -0.22986
101 101.2079725 -0.20797
102 102.1855332 -0.18553
103 103.1514857 -0.15149
104 104.1285342 -0.12853
105 105.1248622 -0.12486
106 106.1111395 -0.11114
107 107.0767472 -0.07675
108 108.0531438 -0.05314
109 109.0329303 -0.03293
matlab：
clc;
clear;
times = 10; %验证不同频率w的次数
term = 100; %测试集频率w的起始点
dterm = 1; %测试集频率w的间距
real_value_w = (term : dterm : term + dterm*(times-1))';%测试集频率w
get_w = (1:times)';%结果
d_w = (1:times)';%误差
for column= 1 : times
% column
%产生原始信号//表达式是固定的
x_prim=2*pi()*1000/1590:0.0002:2*pi()*1000/1510; %原始信号时间序列
simdata = cos(100*x_prim+0*pi()/180); %原始信号
%进行插值到2*N-1个点
N = 2048-1; %目标插值采样间隔 204
%插值
data_x_end = x_prim(end);
data_x_start = x_prim(1);
d_u = (data_x_end - data_x_start)/(N - 1);
x_interp = data_x_start : d_u : data_x_end;
fs = 1/d_u;
Y_interp = interp1(x_prim,simdata,x_interp,'spline');
%插值后的信号以x_interp、Y_interp表示
%%%%
%使用apfft解算原始信号的real_value_w的值
%准备数据
fs=1000; %采样频率
nfft=1024; %FFT计算点数
n=2*nfft-1; %数据采样点数
f1 = real_value_w(column);
% x = cos(2 * pi *x_interp * f1 / fs);
%cos(real_value_w(column)*x_interp+0*pi()/180);
x = Y_interp;
% =============加窗预处理================
win=hanning(nfft)';
win1=conv(win,win);
win1=win1/sum(win1);
win2=win/sum(win);
%加窗
y1=x.*win1;
y2=x(nfft:2*nfft-1).*win2;
%移位相加
y1a=y1(nfft:end) + [0 y1(1:nfft-1)];
%计算FFT
Y1=fft(y1a, nfft);
Y2=fft(y2, nfft);
%幅度
A1=abs(Y1);
A2=abs(Y2);
%相位
P1=mod(angle(Y1)*180/pi,360);
P2=mod(angle(Y2)*180/pi,360);
%频偏
fe=mod((P2-P1)/180/(1-1/nfft),1);
ae=(A2.^2)./A1*2;
%校正后的值
r=round(f1*nfft/fs);
fre=(floor(f1*nfft/fs)+fe(r+1))*fs/nfft;
a=ae(r+1);
p=P1(r+1);
get_w(column)=fre;
d_w(column)=real_value_w(column)-fre;
end
real_value_w %未知信号的频率
get_w %求解出的频率
d_w %误差
我想知道在求解未知信号的频率（范围30：300）时，这个f1、fs怎么设置呀，上面的代码里面f1被设置成了real_value_w(column)也就是生成未知信号的频率，但在实际应用中，并不知道未知信号的频率是多少，所以怎么弄呀？
