最近使用王兆华教授的全相位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)也就是生成未知信号的频率，但在实际应用中，并不知道未知信号的频率是多少，所以怎么弄呀？

最好解决后能实现：对于任意的未知频率（范围30：300）的信号解调后误差小于0.001

如您已解决请联系我修改酬金至100~500(看最后的误差多少，误差小于0.1时100,误差小于0.01时200，误差小于0.001时500)

这个不太了解，但感觉你的数据长度1024是不是有点小了，频率分辨率有点低。不知道你说的精度指什么，这个全相位FFT很出名吗……可能未必效果很好呢


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=parzenwin(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             %误差