帮忙调试一个振动信号频域积分的程序！

unknowno

源程序如下：
数据以及数据输入格式见附件中的txt文件
振动的幅值加速度为5m/ss，频率为35Hz
运行出的结图见附件

问题是：
理论上加速度信号进行一次积分后，幅值加速度和幅值速度的关系应该是：v=a/(2*pi*f)
但一次积分出来结果完全不同！不过这个程序进行二次积分结果还是蛮吻合的。
请不惜赐教！

1. %频域积分
   
2. %%%%%%%%%%%%%%%%%%%%%%%%%%
   
3. clear; clc; close all hidden
   
4. %%%%%%%%%%%%%%%%%%
   
5. fni=input('频域积分-输入数据文件名：','s');
   
6. fid=fopen(fni,'r');
   
7. sf=fscanf(fid,'%f',1);%采样频率
   
8. fmin=fscanf(fid,'%f',1);%最小截止频率
   
9. fmax=fscanf(fid,'%f',1);%最大截止频率
   
10. c=fscanf(fid,'%f',1);%单位变换系数
    
11. it=fscanf(fid,'%f',1);%积分次数
    
12. sx=fscanf(fid,'%s',1);%横向坐标轴的标注
    
13. sy1=fscanf(fid,'%s',1);%纵向坐标轴输入单位的标注
    
14. sy2=fscanf(fid,'%s',1);%纵向坐标轴输出单位的标注
    
15. fno=fscanf(fid,'%s',1);%输出数据文件名
    
16. x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
    
17. status=fclose(fid);
    
18. n=length(x);
    
19. %建立时间向量
    
20. t=0:1/sf:(n-1)/sf;
    
21. %大于并最接近n的2的幂次方为FFT长度
    
22. nfft=2^nextpow2(n);
    
23. %FFT变换
    
24. y=fft(x,nfft);
    
25. %计算频率间隔（Hz/s）
    
26. df=sf/nfft;
    
27. %计算指定频带对应频率数组的下标
    
28. ni=round(fmin/df+1);
    
29. na=round(fmax/df+1);
    
30. %计算圆频率间隔（rad/s）
    
31. dw=2*pi*df;
    
32. %建立正的离散圆频率向量
    
33. w1=0:dw:2*pi*(0.5*sf-df);
    
34. %建立负的离散圆频率向量
    
35. w2=2*pi*(0.5*sf-df):-dw:0;
    
36. %将正负圆频率向量组合成一个向量
    
37. w=[w1,w2];
    
38. %以积分次数为指数，建立圆频率变量向量
    
39. w=w.^it;
    
40. %进行积分的频域变换
    
41. a=zeros(1,nfft); a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);
    
42. if it == 2
    
43.    y=-a; %进行二次积分的相位变换
    
44. else
    
45.    a1=imag(a); a2=real(a); y=a1-a2*i; %进行一次积分的相位变换
    
46. end
    
47. a=zeros(1,nfft);
    
48. %消除指定正频带外的频率成分
    
49. a(ni:na)=y(ni:na);
    
50. %消除指定负频带外的频率成分
    
51. a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
    
52. y=ifft(a,nfft); %IFFT变换
    
53. %取逆变换的实部n个元素并乘以单位变换系数为积分结果
    
54. y=real(y(1:n))*c;
    
55. subplot(2,1,1); plot(t,x); xlabel(sx); ylabel(sy1); grid on; %绘制几分钱的时程曲线图形
    
56. subplot(2,1,2); plot(t,y); xlabel(sx); ylabel(sy2); grid on; %绘制积分后的时程曲线图形
    
57. %打开文件输出积分后的数据
    
58. fid=fopen(fno,'w');
    
59. for k=1:n, fprintf(fid,'%f \n',y(k)); end
    
60. status=fclose(fid);

复制代码

VibrationMaster

原帖由 unknowno 于 2007-7-1 21:30 发表
理论上加速度信号进行一次积分后，幅值加速度和幅值速度的关系应该是：v=a/(2*pi*f)
但积分出来结果完全不同！不过这个程序进行二次积分结果还是蛮吻合的。
请不惜赐教！

理论上这个关系仅对单一频率有效，而你的信号大概不是这样

unknowno

我的信号就是正弦信号，是实测模拟正弦振动的振动台的信号！
再者，用这个程序对这个信号进行二次积分是，结果与理论就蛮吻合的！如下图：
振动频率为35，振动幅值大约为5m/ss，位移x=a/w^2=5/(2*pi*35)^2=1.0339e-004。

[ 本帖最后由 unknowno 于 2007-7-2 14:47 编辑 ]

songzy41

楼主的程序是取自王济和胡晓编的 “MATLAB在振动信号处理中的应用”一书的程序5.5。细看一下程序，发现在圆频率设置上有错误。原程序为

1. %建立正的离散圆频率向量
   
2. w1=0:dw:2*pi*(0.5*sf-df);
   
3. %建立负的离散圆频率向量
   
4. w2=2*pi*(0.5*sf-df):-dw:0;
   
5. 可看到在正频率部分只有nff/2个值，应该有nff/2+1个值，对应于0~sf/2；负频率部分都不是负值，而是正值，笫1个负频率是在nff/2+2处，对应的频率是-(sf/2-df)，然后以df递增，直到nff处对应的频率为-df。所以我把该段程序改为：
   
6. %建立正的离散圆频率向量
   
7. w1=0:dw:2*pi*0.5*sf;
   
8. %建立负的离散圆频率向量
   
9. w2=-2*pi*(0.5*sf-df):dw:-dw;
   
10. %将正负圆频率向量组合成一个向量
    
11. w=[w1,w2];

复制代码

w1有2049个，w2有2047个，w还是4096个。这样求出的速度如下图。程序在附件中。

unknowno

谢谢“songzy41” 的热心帮助！
不知你是否尝试使用改正后的程序进行二次计分运算了？我试了下，一次积分是OK，不过二次积分就不OK了！
下面引用“VibrationMaster”指点的原话：“逆完FFT之后取虚部。因为速度与位移之间的关系为jw，而不是简单的w。”这样对程序进行更正结果是正确的！
不知楼上怎么看？

songzy41

在上程序中把it=2，便得到二次积分的结果，与楼主#5层的结果相符（见下图），不知楼主怎么认为“不OK”？
关于“逆完FFT之后取虚部。因为速度与位移之间的关系为jw，而不是简单的w。”这是对的。我们在把加速度积分成速度时就这样处理过了，细看一下程序，有相应的相位变换：
%进行二次积分的相位变换
   y=-a;
else
%进行一次积分的相位变换
a1=imag(a);
a2=real(a);
y=a1-a2*i;
end
二次积分时因为是-w^2，所以有 y=-a；而一次积分时考虑到jw，把实部和虚部对换。

unknowno

嘿嘿！不好意思！没看见这行程序it=1;我还是在数据文件里面该的t！

谢谢！大家的热心帮助！谢谢“songzy41 ”和“VibrationMaster”！

我是才接触振动方面的东西，好多东西只知道应用并不知道其原理，以后一定把这块给补上去！

[ 本帖最后由 unknowno 于 2007-7-2 19:22 编辑 ]

junzifei

改了之后好像还有问题
大家看一下, 那个地方有问题了又？
大家都说频域积分好，但是这里频域积分的优点没体现出来，并且还有下降的趋势（如图）
初始点也不对？ 为之奈何？哈哈

1. clear; close all; clc
   
2. sf=300; f1=6;f2=10;f3=20;f4=30;
   
3. t=0:1/sf:(1024*3-1)/sf; x=sin(f3*t)+sin(f1*t)+sin(f4*t);
   
4. yyy=-cos(f3*t)/f3-cos(f1*t)/f1-cos(f4*t)/f4+1/f3+1/f1+1/f4;
   
5. t1=1/sf;
   
6. %%%%%辛普森（simpson）算法时域积分求速度
   
7. yvs(1)=t1*(x(1)+x(2))/2; n=length(t);
   
8. for k=2:n-1
   
9. yvs(k)=yvs(k-1)+t1*(x(k-1)+4*x(k)+x(k+1))/6;
   
10. end
    
11. yvs(n)=yvs(n-1);
    
12. %%%%%辛普森（simpson）算法时域积分求速度
    
13. figure
    
14. subplot 211; title('实际信号输入信号'); plot(t,x)
    
15. subplot 212; plot(t,yvs); hold on; plot(t,yyy,'r'); title('时域积分与实际积分信号')
    
16. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
17. nfft=n; c=1;
    
18. fmin=0;%最小截止频率
    
19. fmax=150;%最大截止频率
    
20. it=1;%积分次数
    
21. y=fft(x,nfft); %FFT变换
    
22. df=sf/nfft;%计算频率间隔（Hz/s）
    
23. ni=round(fmin/df+1); na=round(fmax/df+1); %计算指定频带对应频率数组的下标
    
24. dw=2*pi*df; %计算圆频率间隔（rad/s）
    
25. w1=0:dw:2*pi*0.5*sf;%建立正的离散圆频率向量
    
26. w2=-2*pi*(0.5*sf-df):dw:-dw;%建立负的离散圆频率向量
    
27. w=[w1,w2];%将正负圆频率向量组合成一个向量
    
28. w=w.^i; %以积分次数为指数，建立圆频率变量向量
    
29. %进行积分的频域变换
    
30. a=zeros(1,nfft); a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);
    
31. if it == 2
    
32.    y=-a; %进行二次积分的相位变换
    
33. else
    
34. %进行一次积分的相位变换
    
35. a1=imag(a); a2=real(a); y=a1-a2*i;
    
36. end
    
37. a=zeros(1,nfft);
    
38. a(ni:na)=y(ni:na);%消除指定正频带外的频率成分
    
39. a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1); %消除指定负频带外的频率成分
    
40. y=ifft(a,nfft); %IFFT变换
    
41. y=real(y(1:n))*c; %取逆变换的实部n个元素并乘以单位变换系数为积分结果
    
42. %绘制积分前的时程曲线图形
    
43. figure
    
44. subplot(2,1,1); plot(t,x); title('实际时程信号'); grid on;
    
45. %绘制积分后的时程曲线图形
    
46. subplot(2,1,2); plot(t,y); hold on; plot(t,yyy,'r'); hold on; plot(t,yvs,'y')
    
47. legend('频域积分','实际积分','时域积分'); title('频域积分、时域积分与实际积分信号')
    
48. grid on;
    

复制代码

songzy41

原帖由 junzifei 于 2007-7-26 15:14 发表
改了之后好像还有问题
大家看一下
那个地方有问题了又？
大家都说频域积分好，
但是这里频域积分的优点没体现出来，
并且还有下降的趋势（如图）
初始点也不对？
为之奈何？哈哈

我实际上在本帖分析中和http://www.chinavib.com/forum/thread-48861-1-1.html帖子的分析中都发现在频域积分后有误导入一个低频振荡信号。在鸭鸭的那个帖子上很明显。仔细考虑和分析了，认为在频域积分时乘了1/jω因子，当ω很小时，1/jω就很大。同时拿 junzifei 的例子来说，信号中的最低频率是6，但由于于FFT处理中泄漏的存在，低于6的不都为0，使乘了1/jω，可能会使某些频率在频率维有较大的峰值，经IFFT后在时间域上引入了低频振荡，这便是误导入的原因。要解决的办法是引入一个低频滤波，例如在 junzifei 的程序中，把fmin=0.5，这样便得下图，可以看到频率域的积分和时间域的积分只差一个常数。

junzifei

那么常数项该怎么解决呢？
仔细看上面的的图，
还是存在下降趋势，是不是要把最小频率取得更小一点呢？
据说频域积分本身具有滤波性质，难道与相差为常数项的原因有关？

zzw1983107

麻烦问一下：
%建立正的离散圆频率向量
w1=0:dw:2*pi*0.5*sf;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%将正负圆频率向量组合成一个向量

建立负频率应该怎么理解？

fitmos

时域积分的常数项如何解决？？

zhaoch

你好，请问一下在“a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);”
中，为什么数组下标是从（2：nfft -1）开始，这是代表什么物理意义呢？

谢谢

xjtu713

学习！！！！！

volcanostar

真是太有用处了啊