信号处理方法 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

Catsayer 发表于 2015-11-28 11:47

全矢谱、全息谱、全频谱、全矢功率谱、全矢倒频谱程序

程序未经验证，如有需要请自行验证

x，y信号时域和频域程序

clear
clc
load('gu.txt','r');%导入数据
x=gu(:,3);
y=gu(:,4);
N=length(x);%采样点数
fs=2048;%采样频率
t=0:1/fs:(N-1)/fs;
f=0:fs/N:fs/2-fs/N;
%双通道信号时域图
subplot(221);
plot(t,x);
title('x时域图');
xlabel('时间/s');
ylabel('幅值/um');
grid on;
axis();

subplot(222);
plot(t,y);
title('y时域图');
xlabel('时间/s');
ylabel('幅值/um');
grid on;
axis();
%双通道信号频域图
x1=fft(x);
X=abs(x1)*2/N;
subplot(223);
plot(f,X(1:N/2));
title('x幅值谱');
xlabel('频率/Hz');
ylabel('幅值/um');
grid on;
axis(); %取前1000个值显示

y1=fft(y);
Y=abs(y1)*2/N;
subplot(224);
plot(f,Y(1:N/2));
title('y幅值谱');
xlabel('频率/Hz');
ylabel('幅值/um');
grid on;
axis(); %取前1000个值显示

Catsayer 发表于 2015-11-28 11:47

全矢谱程序
clear;
clc;
load('gu.txt');
x=gu(:,3);%导入x信号
y=gu(:,4);%导入y信号
N=1024;fs=2048;
df=fs/N;
z=x+i*y;
Z=fft(z)*2;
for k=1:N/2
Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
RL(k)=Xp(k)+Xr(k);%椭圆长轴
RS(k)=Xp(k)-Xr(k);%椭圆短轴
Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
end
f=df*(0:N/2-1);
subplot(2,1,1);
plot(f,RL);
xlabel('频率');ylabel('主振矢幅值/um');
%set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
axis();
grid on ;
title('全矢谱');
subplot(2,1,2),plot(f,Pha);
xlabel('频率');ylabel('相位角度/（°）');
%set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
set(gca,'YTick',(-180:90:180),'YTickLabel',{'-180','-90','0','90','180'});
axis();grid on

Catsayer 发表于 2015-11-28 11:48

全息谱程序
clear;
clc;
load('gu.txt');
x=gu(:,3);%导入x信号
y=gu(:,4);%导入y信号
N=1024;fs=2048;
df=fs/N;
%全息谱
z=x+i*y;
Z=fft(z)*2;
for k=1:N/2
Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
RL(k)=Xp(k)+Xr(k);%椭圆长轴
RS(k)=Xp(k)-Xr(k);%椭圆短轴
Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
end
f=df*(0:N/2-1);
for i=1:8;
t=1:128;
e=20*i;
W=(2*pi)/128;
q=pi/180;
X=Xp(e)*cos(W*t+PhaP(e)*q)+Xr(e)*cos(W*t+PhaR(e)*q)+e;
Y=Xp(e)*sin(W*t+PhaP(e)*q)-Xr(e)*sin(W*t+PhaR(e)*q);
plot(X,Y);% 画出全息谱图
hold on
end
grid on;
xlabel('倍频');
ylabel('幅值/μm');
set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x','6x'});
axis();

Catsayer 发表于 2015-11-28 11:48

全频谱程序
clear;
clc;
load('gu.txt');
x=gu(:,3);%导入x信号
y=gu(:,4);%导入y信号
N=1024;fs=2048;
df=fs/N;
%全频谱
z=x+i*y;
Z=fft(z)*2;
for k=1:N/2
Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
end
f=df*(0:N/2-1);
plot(f,Xp);
hold on;
plot(-f,Xr);
xlabel('频率');
ylabel('幅值/μm');
%set(gca,'XTick',(-80:20:80),'XTickLabel',{'-4x','-3x','-2x','-x','0','x','2x','3x','4x'});
axis([-1000,1000,0,50]);
title('全频谱')
grid on

Catsayer 发表于 2015-11-28 11:49

全矢功率谱程序
clear;
clc;
load('gu.txt');
x=gu(:,3);%导入x信号
y=gu(:,4);%导入y信号
N=1024;fs=2048;
F1=fft(x);
P1=F1.*conj(F1)/N;
F2=fft(y);
P2=F2.*conj(F2)/N;
df=fs/N;
z=x+i*y;
Z=fft(z)*2;
for k=1:N/2
Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
RL(k)=Xp(k)+Xr(k);%椭圆长轴
RS(k)=Xp(k)-Xr(k);%椭圆短轴
Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
P(k)=RL(k)*RL(k)+RS(k)*RS(k);%功率谱
end
f=df*(0:N/2-1);
subplot(3,1,1);%x通道信号功率谱
plot(f,log10(P1(1:N/2)));
title('x通道信号功率谱');
xlabel('频率');ylabel('E/dB');
axis();
grid on;
subplot(3,1,2);%y通道信号功率谱
plot(f,log10(P2(1:N/2)));
title('y通道信号功率谱');
xlabel('频率');ylabel('E/dB');
axis();
grid on;
subplot(3,1,3);%全矢功率谱
plot(f,log10(P));
title('全矢功率谱');
xlabel('频率');ylabel('E/dB');
%set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
axis();
grid on

Catsayer 发表于 2015-11-28 11:49

全矢倒频谱程序
clear;
clc;
load('gu.txt');
x=gu(:,3);%导入x信号
y=gu(:,4);%导入y信号
N=1024;fs=2048;
F1=fft(x);
P1=F1.*conj(F1)/N;
F2=fft(y);
P2=F2.*conj(F2)/N;
df=fs/N;
z=x+i*y;
Z=fft(z)*2;
for k=1:N/2
Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
RL(k)=Xp(k)+Xr(k);%椭圆长轴
RS(k)=Xp(k)-Xr(k);%椭圆短轴
Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
P(k)=RL(k)*RL(k)+RS(k)*RS(k);%功率谱
end
f=df*(0:N/2-1);
n=0:N/4-1;
subplot(3,1,1);%x通道信号倒频谱
Gx1=abs(fft(abs(log10(P1))));
plot(n,log10(Gx1(1:N/4)));
title('x通道信号倒频谱谱');
xlabel('倒频率');ylabel('FVCP');
axis();
grid on;
subplot(3,1,2);%y通道信号倒频谱
Gx2=abs(fft(abs(log10(P2))));
plot(n,log10(Gx2(1:N/4)));
title('y通道信号倒频谱');
xlabel('倒频率');ylabel('FVCP');
axis();
grid on;
subplot(3,1,3);%全矢倒频谱
Gx=abs(fft(abs(log10(P))));
plot(n,log10(Gx(1:N/4)));
title('全矢倒频谱');
xlabel('倒频率');ylabel('FVCP');
%set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
axis();
grid on ;

13162515140 发表于 2015-11-28 12:37

挺深奥的！

dsp2008 发表于 2015-11-29 09:09

看上去很玄，造文章用的？

dsp2008 发表于 2015-11-29 09:12

这个玩意，感觉与天津大学的那个“全相位FFT” 神似。

又一个“中国创造”？

猫头鹰先生 发表于 2015-11-29 12:33

dsp2008 发表于 2015-11-29 09:09
看上去很玄，造文章用的？

你这个人，可否积点口德！感觉不好就不要看，有多远滚多远

猫头鹰先生 发表于 2015-11-29 12:38

分享是美德！大赞楼主

猫头鹰先生 发表于 2015-11-29 15:56

dsp2008 发表于 2015-11-29 09:09
看上去很玄，造文章用的？

呵呵，我看是戳了你的笼子了吧，这就叫不打自招。

猫头鹰先生 发表于 2015-11-30 13:36

dsp2008 发表于 2015-11-29 09:09
看上去很玄，造文章用的？

跟一个住笼子的没什么还说的。

woshiqiao 发表于 2015-11-30 15:20

好贴！{:{39}:}

osbertbovey 发表于 2015-11-30 17:24

感谢楼主分享！！！！！！！！！！！

页: [1] 2

声振论坛's Archiver

全矢谱、全息谱、全频谱、全矢功率谱、全矢倒频谱程序