1 基本概念

相位：在函数y=Acos(ωx+φ)中，ωx+φ称为相位。

相位差：简谐运动中，如果两个简谐运动的频率相等，其初相位分别是φ1，φ2。当φ2>φ1时，他们的相位差是：△φ=(ωt+φ2)-(ωt+φ1)=φ2-φ1。

2 相位差的计算

根据信号处理方式分类，相位差的估计方法主要有时域法和频域法两类，频域法主要用到FFT，而时域法中比较常见的一种方法为互相关函数法。

算例：振动频率f=10Hz；采样频率Fs=1000；相位差phi=π/2；数据长度N=1024；

2.1 频域法

步骤：FFT→求能量最大点的相位（phase）→计算相位差

2.1.1 代码

% 测试相位差——频域法

% 2019-2-18

% By zhipeng feng

% 算例：振动频率f=10Hz；采样频率Fs=1000；相位差phi=π/2；数据长度N=1024；

clc; clear; close all

f=10;

Fs=1000;

phi=pi/2;

N=1024*2;

f1=10;

f2=20;

t=(0:N-1) / Fs; %首先生成时间序列

y1=sin(2*pi*f*t);

y2=sin(2*pi*f*t+phi);

figure(1)

plot(t,y1,'r',t,y2,'b');

legend('信号1','信号2');

grid

%FFT

frequency=Fs*(0:N/2)/N;  %单边功率谱

fft_y1=fft(y1);  

fft_y2=fft (y2);

%求能量最大点

[~,idx1]=max(abs(fft_y1));

[~,idx2]=max(abs(fft_y2));

% 求两信号能量最大点的相位

phase1=phase(fft_y1(idx1))*180/pi

phase2=phase(fft_y2(idx2))*180/pi

fprintf('信号2比信号1相位差=%f度\n',phase2-phase1);

运行结果：

信号2比信号1相位差=88.651114度

其实，对于同频率的两个向量，其相位滞后dt=0.025s，如下图所示，根据三角函数关系：y1=sin(2*pi*f*t); y2=sin(2*pi*f*(t+dt))，那么y1与y2间的相位差即为：2*pi*f*dt=2*pi*10*0.025=0.5pi=90°。

2.1.2 小结

当采样频率与数据长度相接近时，误差较大，为了提高精度，尽量取较多的数据点。其实只要两者不接近，误差都较小。

取2048个点的结果如下：

信号2比信号1相位差=89.745190度

信号4比信号3相位差=-3.690543度

取512个点的结果如下：

信号2比信号1相位差=89.133260度

信号4比信号3相位差=110.452917度

2.2 时域法

2.2.1 理论背景

互相关函数法：设有A(t)和另一个延迟时间τ的波B(t+τ)，在有限时间间隔T内其互相关函数定义为：

式中τ的取值范围为-T~+T。设τ为变量，则相关函数_RAB就是τ的函数，由相关函数的性质可知：A(t)、B(t)同相时_RAB有最大值，即_RAB取得最大值时，对应的τ即为B(t)相对A(t)的时间延迟。

在实际处理中，总是把A(t)、B(t)进行离散化采集以方便计算机的处理，设A(t)=a(n)、B(t) =b(n)，n=1，2……N；

则互相关函数为：

由于T为有限时间间隔，因此式中当n+m>N时，b(n+m)=0。由上述可知，_rAB是m的函数，当_rAB取最大值时，m对应的τ就是B(t)相对 A(t)的延迟。

2.2.2 代码

取前面的算例进行验证：

%测试相位差——时域法

%2019-2-18

%By zhipeng feng

%算例：振动频率f=10Hz；采样频率Fs=1000；相位差phi=π/2；数据长度N=1024；

clc; clear; close all

f=10;

Fs=1000;

phi=pi/2;

N=1024;

t=(0:N-1) / Fs; %首先生成时间序列

y1=sin(2*pi*f*t);

y2=sin(2*pi*f*t+phi);

n=1:N;

% subplot 211; plot(t,y1,'r',t,y2,'b'); grid;

subplot 211; plot(n,y1,'r',n,y2,'b'); grid;

legend('y1','y2'); title('Signals Waveform');

[acor,lag]=xcorr(y1,y2);%lag=-N-1~N-1

nn=-N+1:N-1;

subplot212; plot(nn,acor); axis tight

grid;title('Correlation Function');

[~,I]= max(abs(acor));

tau=lag(I)/Fs;%数据点除以Fs，道理与生成时间序列t一样

phase=2*pi*f*tau/pi*180

Loc=I-N  %lag=-N-1~N-1，而lag的下标范围是1~2*N-1，因此要减去N，即为图中最大值对应的位置

运行结果如下：

phase =

  -90

Loc =

  -25

2.2.3 还有一种代码

将2.2.2的代码合并到了一起。

% 测试相位差——时域法

% 2019-2-18

% By zhipengfeng

% 算例：振动频率f=10Hz；采样频率Fs=1000；相位差phi=π/2；数据长度N=1024；

clc; clear; close all

f=10;

Fs=1000;

phi=pi/2;

N=1024;

t=(0:N-1) / Fs; %首先生成时间序列

y1=sin(2*pi*f*t);

y2=sin(2*pi*f*t+phi);

figure(1)

n=1:N;

% subplot 211; plot(t,y1,'r',t,y2,'b'); grid;

subplot 211; plot(n,y1,'r',n,y2,'b'); grid;

legend('y1','y2'); title('Signals Waveform');

[acor,lag]=xcorr(y1,y2); %lag=-N-1~N-1

nn=-N+1:N-1;

subplot 212; plot(nn,acor); axis tight

grid; title('Correlation Function');

[~,I] = max(abs(acor));

tau= lag(I)/Fs;%数据点除以Fs，道理与生成时间序列t一样

phase=2*pi*f*tau/pi*180

Loc=I-N  %lag=-N-1~N-1，而lag的下标范围是1~2*N-1，因此要减去N，即为图中最大值对应的位置

%%

figure(2)

num=N;

Ix=sum(y1.^2)/num;

Iy=sum(y2.^2)/num;

Ixy=sum(y1.*y2)/num;

c=180*acos(2*Ixy/(4*Ix*Iy)^0.5)/pi;

plot(t,y1,t,y2);

legend('sin(2*pi*f*t)','y2=sin(2*pi*f*t+phi)');

text(0.5,0,strcat('相位差=',strcat(num2str(c),'度')));

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删