参考:点击参考
本人使用方法二,利用公式:
Matable代码如下:
fid1=fopen(strcat('C:\Users\yxz\Desktop\CNG-WPF123\works\20190619_dn80_复相关\DN80\流量-10\','流量10-0通道-增益3-延时150.dat'),'r');
%打开二
fid2=fopen(strcat('C:\Users\yxz\Desktop\CNG-WPF123\works\20190619_dn80_复相关\DN80\流量-10\','流量10-3通道-增益3-延时150.dat'),'r');
data=[];
data2=[];
%转换文件1的AD值
while ~feof(fid1) % while ~feof 表示 若 未读到文件末尾 则 继续 循环
m=fscanf(fid1,'%2x%3x',[1 inf]);%读取文件
data=[data;m]; %将读出的数据存入data矩阵中
end
i = 1;
k = data;
xtemp = data;
while(i<440)
k(i) = data(2*i+1)*256 + data(2*i+2)-8;
xtemp(i) = i;
if(k(i)>32768)
k(i) = -(65536 - k(i));
end
i = i+1;
end
%转换文件2的AD值
while ~feof(fid2) % while ~feof 表示 若 未读到文件末尾 则 继续 循环
m2=fscanf(fid2,'%2x%3x',[1 inf]);%读取文件
data2=[data2;m2]; %将读出的数据存入data矩阵中
end
j = 1;
k2 = data2;
xtemp2 = data2;
while(j<440)
k2(j) = data2(2*j+1)*256 + data2(2*j+2)-8;
xtemp2(j) = j;
if(k2(j)>32768)
k2(j) = -(65536 - k2(j));
end
j = j+1;
end
y1=k;%信号1
y2=k2;%信号2
%%%%
hs=360/(2*pi);% 将弧度转化为度数
Cc=xcorr(y1);%求互相关函数
%用互相关函数数学表达式求相位差:rm(n)=Rxy(n)/sqrt(Rxx(0)*Ryy(0))
C0=sum(y1.*y2);%y1和y2互相关
A0=sum(y1.*y1);%y1自相关
B0=sum(y2.*y2);%y2自相关
a=sqrt(A0);
b=sqrt(B0);
jiaocha=acos(C0/(a*b))*hs %相位差计算公式*hs*fen
参考是利用方法一,直接调用xcorr
Matable代码
T=8;% 采样周期
fen=60;%将度转化为分:1度=60分
hs=360/(2*pi);% 将弧度转化为度数
t=0:2*pi/(n-1):2*T*pi;% 采样数,每隔2*pi/(n-1)采数据从0 到 2*T*pi结束
N=length(t);%采样长度
y1=4*sin(t);%信号1
y2=4*sin(t*pi/(180*6));%信号2
Cc=xcorr(y1);%求互相关函数
%用互相关函数数学表达式求相位差:rm(n)=Rxy(n)/sqrt(Rxx(0)*Ryy(0))
C0=sum(y1.*y2);%y1和y2互相关
A0=sum(y1.*y1);%y1自相关
B0=sum(y2.*y2);%y2自相关
a=sqrt(A0);
b=sqrt(B0);
Zhunquezhi=(pi/(180*6))*hs*fen %准确值
jiaocha=acos(C0/(a*b))*hs*fen %相位差计算公式
jdwucha=jiaocha-Zhunquezhi %绝对误差
xdwucha=((jiaocha-Zhunquezhi)/Zhunquezhi)*100 %相对误差
%用有效值计算比值差
youxiaoy1=sqrt((1/(N-1))*sum(y1.^2)) %计算y1有效值
youxiaoy2=sqrt((1/(N-1))*sum(y2.^2)) %计算y2有效值
bizhicha=abs((youxiaoy1-youxiaoy2)/youxiaoy1)*100%计算比值差
%作图
m=(-N+1):(N-1);
subplot(311);%画到同一个图
plot(t,y1);
xlabel('t');ylabel('y1');
grid;
subplot(312);
plot(t,y2);
xlabel('t');ylabel('y2');
grid;
subplot(313);
plot(m,Cc);
xlabel('m');ylabel('Cc');
grid;