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

Catsayer

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

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

1. 
   
2. clear
   
3. clc
   
4. load('gu.txt','r');%导入数据
   
5. x=gu(:,3);
   
6. y=gu(:,4);
   
7. N=length(x);%采样点数
   
8. fs=2048;%采样频率
   
9. t=0:1/fs:(N-1)/fs;
   
10. f=0:fs/N:fs/2-fs/N;
    
11. %双通道信号时域图
    
12. subplot(221);
    
13. plot(t,x);
    
14. title('x时域图');
    
15. xlabel('时间/s');
    
16. ylabel('幅值/um');
    
17. grid on;
    
18. axis([0,0.2,-100,100]);
    
19. 
    
20. subplot(222);
    
21. plot(t,y);
    
22. title('y时域图');
    
23. xlabel('时间/s');
    
24. ylabel('幅值/um');
    
25. grid on;
    
26. axis([0,0.2,-100,100]);
    
27. %双通道信号频域图
    
28. x1=fft(x);
    
29. X=abs(x1)*2/N;
    
30. subplot(223);
    
31. plot(f,X(1:N/2));
    
32. title('x幅值谱');
    
33. xlabel('频率/Hz');
    
34. ylabel('幅值/um');
    
35. grid on;
    
36. axis([0,1000,0,100]); %取前1000个值显示
    
37. 
    
38. y1=fft(y);
    
39. Y=abs(y1)*2/N;
    
40. subplot(224);
    
41. plot(f,Y(1:N/2));
    
42. title('y幅值谱');
    
43. xlabel('频率/Hz');
    
44. ylabel('幅值/um');
    
45. grid on;
    
46. axis([0,1000,0,100]); %取前1000个值显示

复制代码

Catsayer

全矢谱程序

1. clear;
   
2. clc;
   
3. load('gu.txt');
   
4. x=gu(:,3);%导入x信号
   
5. y=gu(:,4);%导入y信号
   
6. N=1024;fs=2048;
   
7. df=fs/N;
   
8. z=x+i*y;
   
9. Z=fft(z)*2;
   
10. for k=1:N/2
    
11.     Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
    
12.     Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
    
13.     PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
    
14.     PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
    
15.     RL(k)=Xp(k)+Xr(k);%椭圆长轴
    
16.     RS(k)=Xp(k)-Xr(k);%椭圆短轴
    
17.     Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
    
18. end
    
19. f=df*(0:N/2-1);
    
20. subplot(2,1,1);
    
21. plot(f,RL);
    
22. xlabel('频率');ylabel('主振矢幅值/um');
    
23. %set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
    
24. axis([0,1000,0,60]);
    
25. grid on ;
    
26. title('全矢谱');
    
27. subplot(2,1,2),plot(f,Pha);
    
28. xlabel('频率');ylabel('相位角度/（°）');
    
29. %set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
    
30. set(gca,'YTick',(-180:90:180),'YTickLabel',{'-180','-90','0','90','180'});
    
31. axis([0,1000,-180,180]);grid on
    

复制代码

Catsayer

全息谱程序

1. clear;
   
2. clc;
   
3. load('gu.txt');
   
4. x=gu(:,3);%导入x信号
   
5. y=gu(:,4);%导入y信号
   
6. N=1024;fs=2048;
   
7. df=fs/N;
   
8. %全息谱
   
9. z=x+i*y;
   
10. Z=fft(z)*2;
    
11. for k=1:N/2
    
12.     Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
    
13.     Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
    
14.     PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
    
15.     PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
    
16.     RL(k)=Xp(k)+Xr(k);%椭圆长轴
    
17.     RS(k)=Xp(k)-Xr(k);%椭圆短轴
    
18.     Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
    
19. end
    
20. f=df*(0:N/2-1);
    
21. for i=1:8;
    
22. t=1:128;
    
23. e=20*i;
    
24. W=(2*pi)/128;
    
25. q=pi/180;
    
26. X=Xp(e)*cos(W*t+PhaP(e)*q)+Xr(e)*cos(W*t+PhaR(e)*q)+e;
    
27. Y=Xp(e)*sin(W*t+PhaP(e)*q)-Xr(e)*sin(W*t+PhaR(e)*q);
    
28. plot(X,Y);% 画出全息谱图
    
29. hold on
    
30. end
    
31. grid on;
    
32. xlabel('倍频');
    
33. ylabel('幅值/μm');
    
34. set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x','6x'});
    
35. axis([0,200,-5,5]);
    

复制代码

Catsayer

全频谱程序

1. clear;
   
2. clc;
   
3. load('gu.txt');
   
4. x=gu(:,3);%导入x信号
   
5. y=gu(:,4);%导入y信号
   
6. N=1024;fs=2048;
   
7. df=fs/N;
   
8. %全频谱
   
9. z=x+i*y;
   
10. Z=fft(z)*2;
    
11. for k=1:N/2
    
12.     Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
    
13.     Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
    
14. end
    
15. f=df*(0:N/2-1);
    
16. plot(f,Xp);
    
17. hold on;
    
18. plot(-f,Xr);
    
19. xlabel('频率');
    
20. ylabel('幅值/μm');
    
21. %set(gca,'XTick',(-80:20:80),'XTickLabel',{'-4x','-3x','-2x','-x','0','x','2x','3x','4x'});
    
22. axis([-1000,1000,0,50]);
    
23. title('全频谱')
    
24. grid on
    

复制代码

Catsayer

全矢功率谱程序

1. clear;
   
2. clc;
   
3. load('gu.txt');
   
4. x=gu(:,3);%导入x信号
   
5. y=gu(:,4);%导入y信号
   
6. N=1024;fs=2048;
   
7. F1=fft(x);
   
8. P1=F1.*conj(F1)/N;
   
9. F2=fft(y);
   
10. P2=F2.*conj(F2)/N;
    
11. df=fs/N;
    
12. z=x+i*y;
    
13. Z=fft(z)*2;
    
14. for k=1:N/2
    
15.     Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
    
16.     Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
    
17.     PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
    
18.     PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
    
19.     RL(k)=Xp(k)+Xr(k);%椭圆长轴
    
20.     RS(k)=Xp(k)-Xr(k);%椭圆短轴
    
21.     Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
    
22.     P(k)=RL(k)*RL(k)+RS(k)*RS(k);%功率谱
    
23. end
    
24. f=df*(0:N/2-1);
    
25. subplot(3,1,1);%x通道信号功率谱
    
26. plot(f,log10(P1(1:N/2)));
    
27. title('x通道信号功率谱');
    
28. xlabel('频率');ylabel('E/dB');
    
29. axis([0,1000,-5,10]);
    
30. grid on;
    
31. subplot(3,1,2);%y通道信号功率谱
    
32. plot(f,log10(P2(1:N/2)));
    
33. title('y通道信号功率谱');
    
34. xlabel('频率');ylabel('E/dB');
    
35. axis([0,1000,-5,10]);
    
36. grid on;
    
37. subplot(3,1,3);%全矢功率谱
    
38. plot(f,log10(P));
    
39. title('全矢功率谱');
    
40. xlabel('频率');ylabel('E/dB');
    
41. %set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
    
42. axis([0,1000,-5,10]);
    
43. grid on
    

复制代码

Catsayer

全矢倒频谱程序

1. clear;
   
2. clc;
   
3. load('gu.txt');
   
4. x=gu(:,3);%导入x信号
   
5. y=gu(:,4);%导入y信号
   
6. N=1024;fs=2048;
   
7. F1=fft(x);
   
8. P1=F1.*conj(F1)/N;
   
9. F2=fft(y);
   
10. P2=F2.*conj(F2)/N;
    
11. df=fs/N;
    
12. z=x+i*y;
    
13. Z=fft(z)*2;
    
14. for k=1:N/2
    
15.     Xp(k)=abs(Z(k))/(2*N);%正进动圆半径
    
16.     Xr(k)=abs(Z(N+1-k))/(2*N);%反进动圆半径
    
17.     PhaP(k)=angle(Z(k))*180/pi;%正进动圆起始相位角
    
18.     PhaR(k)=angle(Z(N+1-k))*180/pi;%反进动圆起始相位角
    
19.     RL(k)=Xp(k)+Xr(k);%椭圆长轴
    
20.     RS(k)=Xp(k)-Xr(k);%椭圆短轴
    
21.     Pha(k)=(PhaP(k)+PhaR(k))/2;%相位角
    
22.     P(k)=RL(k)*RL(k)+RS(k)*RS(k);%功率谱
    
23. end
    
24. f=df*(0:N/2-1);
    
25. n=0:N/4-1;
    
26. subplot(3,1,1);%x通道信号倒频谱
    
27. Gx1=abs(fft(abs(log10(P1))));
    
28. plot(n,log10(Gx1(1:N/4)));
    
29. title('x通道信号倒频谱谱');
    
30. xlabel('倒频率');ylabel('FVCP');
    
31. axis([0,250,-3,3]);
    
32. grid on;
    
33. subplot(3,1,2);%y通道信号倒频谱
    
34. Gx2=abs(fft(abs(log10(P2))));
    
35. plot(n,log10(Gx2(1:N/4)));
    
36. title('y通道信号倒频谱');
    
37. xlabel('倒频率');ylabel('FVCP');
    
38. axis([0,250,-3,3]);
    
39. grid on;
    
40. subplot(3,1,3);%全矢倒频谱
    
41. Gx=abs(fft(abs(log10(P))));
    
42. plot(n,log10(Gx(1:N/4)));
    
43. title('全矢倒频谱');
    
44. xlabel('倒频率');ylabel('FVCP');
    
45. %set(gca,'XTick',(0:20:100),'XTickLabel',{'0','x','2x','3x','4x','5x'});
    
46. axis([0,250,-3,3]);
    
47. grid on ;
    

复制代码

13162515140

挺深奥的！

dsp2008

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

dsp2008

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

又一个“中国创造”？

猫头鹰先生

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

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

猫头鹰先生

分享是美德！大赞楼主

猫头鹰先生

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

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

猫头鹰先生

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

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

woshiqiao

好贴！

osbertbovey

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