相位差校正法程序.txt

该文件是用相位差校正技术来检测简谐波

%  本程序理论依据：《通用的离散频谱相位差校正方法》，丁康 %



N=256;
m=256;
f=49.4988;
t1=4/f;
global k;
L=256;
%t=0:t1/N:(t1-t1/N);
fs=3167.9232;
n=0:1:N-1;
x1=310*cos(2*pi*49.5*n/fs+(pi/36))+6.2*cos(2*pi*2*49.5*n/fs+(pi/6))+93*cos(2*pi*3*49.5*n/fs+(pi/9))+62*cos(2*pi*5*49.5*n/fs+(pi/2))+37.2*cos(2*pi*7*49.5*n/fs-(5*pi/6))+15.5*cos(2*pi*13*49.5*n/fs+(2*pi/3));
                              %第一段采样后的信号
x2=310*cos(2*pi*(n+L)*49.5/fs+(pi/36))+6.2*cos(2*pi*(n+L)*2*49.5/fs+(pi/6))+93*cos(2*pi*(n+L)*3*49.5/fs+(pi/9))+62*cos(2*pi*(n+L)*5*49.5/fs+(pi/2))+37.2*cos(2*pi*(n+L)*7*49.5/fs-(5*pi/6))+15.5*cos(2*pi*(n+L)*13*49.5/fs+(2*pi/3));
                              %滞后第一段信号LTs的采样信号
for i=0:N-1
     %w(i+1)=1;           %矩形窗
     w(i+1)=0.5*(1-cos(2*pi*i/(N-1)));    %hanning窗
     w1(i+1)=0.54-0.46*cos(2*pi*i/(N-1));   %hamming窗
     %w(i+1)=0.42-0.5*cos(2*pi*i)/(N-1))+0.08*cos(4*pi*i/(N-1));  %blackman窗
end
xx1=x1.*w;
xx2=x2.*w1;
y1=fft(xx1,m);
y2=fft(xx2,m); 
I1=imag(y1(1:m/2));
R1=real(y1(1:m/2));
z1=atan(I1./R1);           % 第一段信号相位角
I2=imag(y2(1:m/2));
R2=real(y2(1:m/2));
z2=atan(I2./R2);         %第二段信号相位角
y3=abs(y1(1:m/2))*4/N;
y3(1)=y3(1)/2;
y4=abs(y2(1:m/2))*4/N;
y4(1)=y4(1)/2;
q1=angle(y1(1:m/2))*180/pi; % 加窗FFT后的相位角
%q2=angle(y2(1:m/2))*180/pi;
d=(z2-z1)/(2*pi);    %频率校正量 
da=-d*pi;      %相位角校正量
Z=da+q1;      %校正后的相位角
K=0:(m/2-1);
F=(K+d)*fs/N;   %校正后的频率
k=(0:(m/2-1)).*(N/t1)./(N*f); %谐波次数
%subplot(2,2,1);
%plot(x1);
%subplot(2,2,2);
stem(k,y3);
%stem(k,y1,'.');
%subplot(2,2,3);
%plot(x2);
%subplot(2,2,4);
%stem(k,y4);
%stem(k,y2,'.');



N=256;
m=256;
f=50.4994;
t1=4/f;
global k;
L=256;
%t=0:t1/N:(t1-t1/N);
fs=3231.9616;
n=0:1:N-1;
x1=310*cos(2*pi*50.5*n/fs+(pi/36))+6.2*cos(2*pi*2*50.5*n/fs+(pi/6))+93*cos(2*pi*3*50.5*n/fs+(pi/9))+62*cos(2*pi*5*50.5*n/fs+(pi/2))+37.2*cos(2*pi*7*50.5*n/fs-(5*pi/6))+15.5*cos(2*pi*13*50.5*n/fs+(2*pi/3));
                              %第一段采样后的信号
x2=310*cos(2*pi*(n+L)*50.5/fs+(pi/36))+6.2*cos(2*pi*(n+L)*2*50.5/fs+(pi/6))+93*cos(2*pi*(n+L)*3*50.5/fs+(pi/9))+62*cos(2*pi*(n+L)*5*50.5/fs+(pi/2))+37.2*cos(2*pi*(n+L)*7*50.5/fs-(5*pi/6))+15.5*cos(2*pi*(n+L)*13*50.5/fs+(2*pi/3));
                              %滞后第一段信号LTs的采样信号
for i=0:N-1
     %w(i+1)=1;           %矩形窗
     w(i+1)=0.5*(1-cos(2*pi*i/(N-1)));    %hanning窗
     w1(i+1)=0.54-0.46*cos(2*pi*i/(N-1));   %hamming窗
     %w(i+1)=0.42-0.5*cos(2*pi*i)/(N-1))+0.08*cos(4*pi*i/(N-1));  %blackman窗
end
xx1=x1.*w;
xx2=x2.*w1;
y1=fft(xx1,m);
y2=fft(xx2,m); 
I1=imag(y1(1:m/2));
R1=real(y1(1:m/2));
z1=atan(I1./R1);           % 第一段信号相位角
I2=imag(y2(1:m/2));
R2=real(y2(1:m/2));
z2=atan(I2./R2);         %第二段信号相位角
y3=abs(y1(1:m/2))*4/N;
y3(1)=y3(1)/2;
y4=abs(y2(1:m/2))*4/N;
y4(1)=y4(1)/2;
q1=angle(y1(1:m/2))*180/pi; % 加窗FFT后的相位角
%q2=angle(y2(1:m/2))*180/pi;
d=(z2-z1)/(2*pi);    %频率校正量 
da=-d*pi;      %相位角校正量
Z=da+q1;      %校正后的相位角
K=0:(m/2-1);
F=(K+d)*fs/N;   %校正后的频率
k=(0:(m/2-1)).*(N/t1)./(N*f); %谐波次数
%subplot(2,2,1);
%plot(x1);
%subplot(2,2,2);
stem(k,y3);
%stem(k,y1,'.');
%subplot(2,2,3);
%plot(x2);
%subplot(2,2,4);
%stem(k,y4);
%stem(k,y2,'.');