1 基本概念

相位：在函數(shù)y=Acos(ωx+φ)中，ωx+φ稱為相位。

相位差：簡諧運動中，如果兩個簡諧運動的頻率相等，其初相位分別是φ1，φ2。當φ2>φ1時，他們的相位差是：△φ=(ωt+φ2)-(ωt+φ1)=φ2-φ1。

2 相位差的計算

根據(jù)信號處理方式分類，相位差的估計方法主要有時域法和頻域法兩類，頻域法主要用到FFT，而時域法中比較常見的一種方法為互相關(guān)函數(shù)法。

算例：振動頻率f=10Hz；采樣頻率Fs=1000；相位差phi=π/2；數(shù)據(jù)長度N=1024；

 2.1 頻域法

步驟：FFT→求能量最大點的相位（phase）→計算相位差

2.1.1 代碼

% 測試相位差——頻域法

% 2019-2-18

% By zhipeng feng

% 算例：振動頻率f=10Hz；采樣頻率Fs=1000；相位差phi=π/2；數(shù)據(jù)長度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);

 

運行結(jié)果：

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

其實，對于同頻率的兩個向量，其相位滯后dt=0.025s，如下圖所示，根據(jù)三角函數(shù)關(guān)系：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 小結(jié)

當采樣頻率與數(shù)據(jù)長度相接近時，誤差較大，為了提高精度，盡量取較多的數(shù)據(jù)點。其實只要兩者不接近，誤差都較小。

取2048個點的結(jié)果如下：

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

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

取512個點的結(jié)果如下：

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

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

2.2 時域法

2.2.1 理論背景

互相關(guān)函數(shù)法：設(shè)有A(t)和另一個延遲時間τ的波B(t+τ)，在有限時間間隔T內(nèi)其互相關(guān)函數(shù)定義為：

式中τ的取值范圍為-T~+T。設(shè)τ為變量，則相關(guān)函數(shù)R_AB就是τ的函數(shù)，由相關(guān)函數(shù)的性質(zhì)可知：A(t)、B(t)同相時R_AB有最大值，即R_AB取得最大值時，對應(yīng)的τ即為B(t)相對A(t)的時間延遲。

在實際處理中，總是把A(t)、B(t)進行離散化采集以方便計算機的處理，設(shè)A(t)=a(n)、B(t) =b(n)，n=1，2……N；

則互相關(guān)函數(shù)為：

由于T為有限時間間隔，因此式中當n+m>N時，b(n+m)=0。由上述可知，r_AB是m的函數(shù)，當r_AB取最大值時，m對應(yīng)的τ就是B(t)相對 A(t)的延遲。

2.2.2 代碼

取前面的算例進行驗證：

%測試相位差——時域法

%2019-2-18

%By zhipeng feng

%算例：振動頻率f=10Hz；采樣頻率Fs=1000；相位差phi=π/2；數(shù)據(jù)長度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;%數(shù)據(jù)點除以Fs，道理與生成時間序列t一樣

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

Loc=I-N  %lag=-N-1~N-1，而lag的下標范圍是1~2*N-1，因此要減去N，即為圖中最大值對應(yīng)的位置

運行結(jié)果如下：

phase =

   -90

Loc =

   -25

2.2.3 還有一種代碼

將2.2.2的代碼合并到了一起。

% 測試相位差——時域法

% 2019-2-18

% By zhipengfeng

% 算例：振動頻率f=10Hz；采樣頻率Fs=1000；相位差phi=π/2；數(shù)據(jù)長度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;%數(shù)據(jù)點除以Fs，道理與生成時間序列t一樣

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

Loc=I-N  %lag=-N-1~N-1，而lag的下標范圍是1~2*N-1，因此要減去N，即為圖中最大值對應(yīng)的位置

%%

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),'度')));

3 參考資料

1. 羅映霞，馬君，朱青松.基于MATLAB的信號相位差的互相函數(shù)求法[J].科技情報開發(fā)與經(jīng)濟，2013,13(7):149.

2. https://baike.baidu.com/item/%E7%9B%B8%E4%BD%8D/2391710?fr=aladdin

3. https://baike.baidu.com/item/%E7%9B%B8%E4%BD%8D%E5%B7%AE/4671188

4. https://zhidao.baidu.com/question/1950485492629291508.html

5. http://www.ilovematlab.cn/thread-481755-2-1.html

來源：FSI Center、

作者：zhipeng feng