phase.txt, 共357个采样

由於不知第一列时间单位, 设数据取样频率fs=N

用下面fft/apfft法校正,  得图一 N=128时fft/apfft谱分析图

图一a 为 输入信号波形

图一b为 apfft的振幅谱, 在f=-4处有一峰值

图一c为 apfft的相位谱, 在f=4处为p2(4)= 322.482度, apfft测的是序列中点(第128点)的相位值.

图一d为校正频偏谱,在f=4处为ee(4)=0.1984
图一e 为校正振幅谱, 在f=4处为  aa(4)=0.55102

原信号一开始有直流分量, 信号是一个含直流分量0.076975的余弦波, 去掉直流分量,

给的数据的是余弦信号,振幅是0.55102,频率3.1984,中间点相位322.482度, 可用下式产生
t=-127:127;
A=0.55102;f=3.1984;phi=322.482;
Y=A*cos(2*pi*t*f/128+phi*pi/180);

图二为(a)原数据波形和(b)重构波形,除直流和漂移外, 二者波形一致. 表明谱分析正确.
而楼主要求初相位是笫一个采样的相位, 即 t= -127 点的相位(t= -127代入上Y的相位值)

即用语句ppp=mod(p2(4)-fff*127/128*360,360)    得初相位为260.04度

信号第1点相位为260.04 度(校正值为余弦波相位, 相对正弦波是350.04) ,符合图二(b)图示正弦信号的相位

这个实例说明apfft测的是序列中点的相位, 如何求序列第1点的相位

close all;clc;clear all;
N=128;
yy=load('phase.txt');
yy=yy(:,2);
y0=mean(yy)%直流分量
y=yy'-mean(yy);
y=y(1:2*N-1);%2N-1个输入数据
y1 = y(N:2*N-1);%后N个输入数据
win =  hanning(N)';;
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
p1 = mod(angle(y11_fft)*180/pi,360);;%FFT相位谱
y2 = y(1:2*N-1);%2N-1个输入数据
win =  hanning(N)';;
winn =  conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
p2=mod( angle(y2_fft)*180/pi,360);%apFFT相位谱
ee=(p1-p2)/180/(1-1/N);%频率偏离校正值
aa=(a1.^2)./a2*2;%振幅校正值
subplot(5,1,1);stem(y,'.');title('input signal(1;255)');ylim([-1,1]);xlim([0 N*2]);grid
subplot(5,1,2),stem(a2,'.');title('apFFTamplitude spectrum');ylim([0,0.5]);xlim([0 N/2]);grid
subplot(5,1,3),stem(p2,'.');title('128点apFFT phase spectrum');ylim([0,400]);xlim([0 N/2]);grid
subplot(5,1,4);stem(ee,'.');title('frequency correction spectrum');ylim([-1,1]);xlim([0 N/2]);grid
subplot(5,1,5);stem(aa,'.');title('amplitude correction spectrum');ylim([0,1]);xlim([0 N/2]);grid
disp('频率校正值')
fff=3+ee(4)
disp('振幅校正值')
aaa=aa(4)
disp('初相位校正值')
pp=p2(4)%第128点相位
ppp=mod(p2(4)-fff*127/128*360,360) %由128点相位计算第1点相位

                  图一  fft/apfft谱分析图

                           图二  原数据波形和重构波形