%相位差校正法
clear all
close all
clc
L=100; %延迟点数
N=1024;
S=L+N;
n=0:S-1;
n1=0:N-1;
fs=1024;
dt=1/fs;
x=cos(2*pi*200.4*n*dt+80*pi/180); %此处在时域取L+N点
wn=1/2-1/2*cos(2*pi*n1/N); %汉宁窗时域表达式
xh=x(1:N).*wn;
%%
Xk1=FFT(xh,N); %对前N点做FFT，即0到N-1点
y=abs(Xk1); %取幅值
subplot(321)
plot((0:N/2-1)*fs/N,4*y(1:N/2)/N),grid on
xlabel('频率'),ylabel('幅值'),title('平移前幅值图')
subplot(322)
plot((0:N/2-1)*fs/N,angle(Xk1(1:N/2))*180/pi),grid on
xlabel('频率'),ylabel('相位'),title('平移前相位图')
%%
wn=1/2-1/2*cos(2*pi*n1/N); %汉宁窗时域表达式
xh1=x(L+1:S).*wn;
Xk2=FFT(xh1,N); %对后N点做FFT运算，即第L点到L+N-1
y=abs(Xk2)/N;
[y1,index]=max(y);
k=index-1;
subplot(323)
plot((L-1:L+N/2-2)*fs/N,4*y(L:L+N/2-1)),grid on
xlabel('频率'),ylabel('幅值'),title('平移后幅值图')
subplot(324)
plot((L-1:L+N/2-2)*fs/N,180*(angle(Xk2(L:L+N/2-1)))/pi),grid on
xlabel('频率'),ylabel('相位'),title('平移后相位图')
%%
%校正
deltaphi=180*(angle(Xk2(201)))/pi-angle(Xk1(201))*180/pi;
phase1=deltaphi/180*pi;
phase2=2*pi*200*100/1024;
delta=mod(phase1-phase2,2*pi); %和2pi取余
%将sigma调整到（-pi,pi）之间
if delta<-pi
delta1=delt+2*pi;
elseif delta>pi
delta1=delt-2*pi;
else delta1=delta;
end
deltaf=delta1*1024/(2*pi*100);
f0=(k+deltaf)*fs/N;
A=2*y(index)/(sinc(deltaf)/(2*(1-deltaf^2)));%汉宁窗的频谱函数化成sinc形函数
theta=(angle(Xk1(201)))*180/pi-deltaf*180;
%%
%校正后的图形
subplot(325)
stem(f0,A),grid on
xlabel('校正后频率'),ylabel('幅值'),title('校正后幅值图')
subplot(326)
stem(f0,theta),grid on
xlabel('校正后频率'),ylabel('相位'),title('校正后相位图')
%%
