|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
由自功率和互功率求相干系数,但是最终结果有问题,相干系数太小,希望大家能帮指点一下。谢谢
其中1Z,5Z为电动水泵上的两个振动测点的时域信号。
clear;
x1=load('5z.txt');
x2=load('1z.txt');
n=5120;
fs=2560;%sampling frequency
N=5120;%数据点长度
%%%%%%%%%%%%%%%%%%%% 自功率Gxx=A(f).*A'(F) %%%%%%%%%%%%%%%
t=x1([1:5120],1);
y1=x1([1:5120],2);
y2=x2([1:5120],2);
%n=0:length(y)-1;
%f=fs*n/length(y);
Y1=fft(y1(1:N),N)/N;
Y2=fft(y2(1:N),N)/N;
f=fs/2*linspace(0,1,N/2);
%%%%%%%%%求自相关函数,互相关函数
[R1,tao1]=xcorr(y1,5120,'coeff');
YR1=fft(R1(N/2:N),N);
[R2,tao2]=xcorr(y2,5120,'coeff');
YR2=fft(R2(1:N),N)/N;
[R,tao3]=xcorr(y1,y2,5120,'coeff');
YR=fft(R(N/2:N),N)/N; %互相关函数
%%%%%%%%%y1相关变换%%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(311);
plot(t(1:5120),y1(1:5120));%绘制y1原始时域信号
xlabel('时间(t)');
ylabel('幅值');
title('y1时域信号');
subplot(312);
plot(f,abs(Y1(1:N/2)));%绘制y1频谱图
xlabel('频率(Hz)');
ylabel('幅值');
title('y1频谱图');
subplot(313);
plot(f,abs(YR1(1:N/2)));%绘制y1自相关函数
xlabel('频率(Hz)');
ylabel('幅值');
title('y1自相关函数');
%%%%%%%%%y2相关变换%%%%%%%%%%%%%%%%%%%%%%%
figure;
subplot(311);
plot(t(1:5120),y2(1:5120));%绘制y2原始时域信号
xlabel('时间(t)');
ylabel('幅值');
title('y2时域信号');
subplot(312);
plot(f,abs(Y2(1:N/2)));%绘制y2频谱图
xlabel('频率(Hz)');
ylabel('幅值');
title('y2频谱图');
subplot(313);
plot(f,abs(YR2(1:N/2)));%绘制y2自相关函数
xlabel('频率(Hz)');
ylabel('幅值');
title('y2自相关函数');
%%%%%%%%%%%%%%y1y2互相关%%%%%%%
figure;
plot(f,abs(YR(1:N/2)));%自相关函数
title('y1y2互相关');
%%%%%%%%%%%求均方根谱
sq1mm=1000.*Y1;%以mm为单位
sq1=abs(sq1mm);%以m为单位
%求功率谱
Y1power=sq1.^2;
figure;
plot(f,Y1power(1:N/2));
xlabel('频率(Hz)');
ylabel('功率谱');
title('y1信号功率谱');
grid;
%求均方根谱
sq2mm=1000.*Y2;%以mm为单位
sq2=abs(sq2mm);%以m为单位
%求功率谱
Y2power=sq2.^2;
figure;
plot(f,Y2power(1:N/2));
xlabel('频率(Hz)');
ylabel('功率谱');
title('y2信号功率谱');
grid;
%logY2power=20.*log(sq2./10^-6);
%%%%%%%%% Y1 Y2互功率谱 %%%%%%%%%
%求互谱均方根
sq3mm=1000.*YR;%以mm为单位
sq3=abs(sq3mm);%以m为单位
%sq3=abs(Y2);%以m为单位
Y1Y2power=sq3.^2;%以mm为单位
figure;
plot(f,Y1Y2power(1:N/2));
xlabel('频率(Hz)');
ylabel('功率谱');
title('信号互功率谱');
grid;
%%%%%%%%%%%%%%%%% Y1 Y2 相干系数 %%%%%%%%
%r=(Y1Y2power.^2)./(Y1power.*Y2power);
m=Y1Y2power.^2;
n=Y1power.*Y2power;
r=m./n;
figure;
plot(f,r(1:N/2));
xlabel('频率(Hz)');
ylabel('相干');
title('相干系数');
grid; |
|