关于三分之一倍频程中有效值的求法问题

小海豚zc

第一个问题：《MATLAB在振动信号处理中的应用》一书中有三分之一倍频程的程序（程序在下面），请高手帮我看看程序有没有问题，主要是中间求有效值那里

1. %三分之一倍频程处理
   
2. clear
   
3. %clc
   
4. %close all hidden
   
5. format long
   
6. %fni=input('三分之一倍频程处理－输入数据文件名','s');
   
7. %fid=fopen(fni,'r')
   
8. %sf=fscanf(fid,'%f',1); %读入采样频率值
   
9. %fno=fscanf(fid,'%d',1); %读入输出数据文件名
   
10. %x=fscanf(fid,'%f',[1 inf]); %读入输入数据存成列向量
    
11. %status=fclose(fid);
    
12. sf=200;
    
13. %fno='out6_4.mat';
    
14. %load y
    
15. %x=y;
    
16. N=20000;
    
17. nn=1:N;
    
18. t=nn./sf;
    
19. x=3*sin(2*pi*2.25*t)+3*sin(2*pi*15*t)+3*sin(2*pi*23*t);
    
20. %定义三分之一倍频程的中心频率
    
21. f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
    
22. fc=[f,10*f,100*f,1000*f,10000*f];
    
23. %中心频率与下限频率的比值
    
24. oc6=2^(1/6);
    
25. nc=length(fc);
    
26. n=length(x);
    
27. nfft=2^nextpow2(n);
    
28. a=fft(x,nfft);
    
29. for j=1:nc
    
30. fl=fc(j)/oc6;
    
31. fu=fc(j)*oc6;
    
32. nl=round(fl*nfft/sf+1);
    
33. nu=round(fu*nfft/sf+1);
    
34. if fu>sf/2
    
35. m=j-1;break
    
36. end
    
37. b=zeros(1,nfft);
    
38. b(nl:nu)=a(nl:nu);
    
39. b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    
40. c=ifft(b,nfft);
    
41. yc(j)=sqrt(var(real(b(1:n))));%计算每个中心频率段的有效值
    
42. end
    
43. figure
    
44. subplot(2,1,1);
    
45. t=0:1/sf:(n-1)/sf;
    
46. plot(t,x);
    
47. xlabel('时间(s)');
    
48. ylabel('加速度(g)');
    
49. grid on;
    
50. subplot(2,1,2);
    
51. plot(fc(1:m),yc(1:m));
    
52. xlabel('频率(Hz)');
    
53. ylabel('有效值');
    
54. xlim([0 100])
    
55. grid on;
    
56. %fid=fopen(fno,'w');
    
57. %for k=1:m
    
58. % fprintf(fid,'%f%f\n',fc(k),yc(k));
    
59. %end
    
60. %status=fclose(fid);

复制代码

xinglong-liu

好像第41行“yc(j)=sqrt(var(real(b(1:n))));%计算每个中心频率段的有效值”改为“yc(j)=std(c)”即可

小海豚zc

xinglong-liu 发表于 2015-7-16 22:43
好像第41行“yc(j)=sqrt(var(real(b(1:n))));%计算每个中心频率段的有效值”改为“yc(j)=std(c)”即可

首先谢谢你，还有几个问题想请教下
1.c是时域的，b是频域的，看本论坛的帖子说原理上都可以实现，那用b怎么做呢？
2.如果让输入是一个正弦信号，那均方根（有效值）应该是幅值除根2，可是用上面的程序却算不对（接近幅值除根2），误差是因为什么呢？？  是因为FFT的能量泄露吗？？
3.有这样一种想法：既然信号可以做傅里叶级数展开，那我可以先把信号做FFT幅值谱，然后幅值谱除根2，再对每个频带做均方根。。这个做法可行吗？？

TestGuru

非整周期采样都会有泄漏问题，正弦信号用FFT做出来的谱线能量总会有泄漏，泄漏多少看所用的窗函数，能量总小于理论值。

小海豚zc

TestGuru 发表于 2015-7-18 10:12
非整周期采样都会有泄漏问题，正弦信号用FFT做出来的谱线能量总会有泄漏，泄漏多少看所用的窗函数，能量总 ...

恩，有道理
能帮我看看楼上的第三个问题吗？？用这种除根2求有效值的方式可行吗？

rossbin

三分之一倍频程的求法在频域内是相应octave内的均方根值，频域内的计算是平方和开根号即可。

rossbin

小海豚zc 发表于 2015-7-17 09:42
首先谢谢你，还有几个问题想请教下
1.c是时域的，b是频域的，看本论坛的帖子说原理上都可以实现，那用b ...

1、频域计算三分之一倍频程是相应octave内的rms值；
2、对简单的正弦信号，整周期采用，FFT到频谱后，频域的幅值需要scale，貌似要除以N，单边谱的话要除以N再乘以2。需要注意的是时域和频域的均方根值得算法。均方根值就是rms，时域先s，再m，最好r。但频域先s，然后直接r，不需要m，用matlab可以直接验证的。最近不做这个了，应当没有记错。
3、周期性信号才可展开傅里叶级数，非周期的要用傅里叶变换。个人认为不可行，你需要了解时域到频域信号的处理过程。

小海豚zc

rossbin 发表于 2015-7-22 20:06
1、频域计算三分之一倍频程是相应octave内的rms值；
2、对简单的正弦信号，整周期采用，FFT到频谱后，频 ...

1.你说的octive是专门的一个软件吗？还是哪个概念？不好意思，我知道的太少了
2.你说的整周期采样，可以用上面的程序试验下，我自己试过了：一个正弦信号，即使整周期采样，还是达不到幅值除以根2（一个正选信号的有效值不就是幅值除以根2吗？），总是缺少一些，不清楚原因？？？还请麻烦您运行试下，帮我找下到底是什么问题

3.第三问题，其实我也一直有疑惑，有人说“功率信号才能做功率谱，能量信号才能做幅值谱”，可是功率谱有一种做法（Cooley-Turkey法）不就是在幅值谱的基础上调整幅值得到的吗？从程序算法的角度讲，这两种谱的求法一定存在关系，只不过从物理意义上讲就讲不通了。。     
然后回到我的问题，你说的非周期信号不能用fft，那非周期信号怎么用频域法？还有就是有限长度的非周期信号不是也可以拓展成周期信号吗？
不知道我说的对不对，还请您指正

rossbin

小海豚zc 发表于 2015-7-24 19:48
1.你说的octive是专门的一个软件吗？还是哪个概念？不好意思，我知道的太少了
2.你说的整周期采样，可以 ...

你需要多看书。

小海豚zc

rossbin 发表于 2015-7-24 19:54
你需要多看书。

哈哈，听君一席话，胜读十年书啊

我又想了下第三个问题：对一个信号，我没法保证对关注的频率都是整周期采样，所以，这样直接除根2 的算法确实是有问题的，

但是第二个问题还是考虑不清楚，原理上是对的

TestGuru

整周期采样，算出来的正弦波的幅值一定会对的。否则检测程序是否有错误。

非整周期采样，需要采用适当的窗函数把往外泄漏的能量集中到主谱线附近，算出来的正弦波的幅度也一定会对的，否则检测程序是否有错误。

小海豚zc

TestGuru 发表于 2015-7-25 17:08
整周期采样，算出来的正弦波的幅值一定会对的。否则检测程序是否有错误。

非整周期采样，需要采用适当的 ...

您帮忙看一下，我找不到问题

1. %三分之一倍频程处理
   
2. clear
   
3. %clc
   
4. %close all hidden
   
5. format long
   
6. 
   
7. sf=2000;
   
8. N=100000;
   
9. nn=1:N;
   
10. t=nn./sf;
    
11. x=2*sin(2*pi*2*t)+2*sqrt(2)*sin(2*pi*10*t)+2*sin(2*pi*30*t);
    
12. %定义三分之一倍频程的中心频率
    
13. f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
    
14. fc=[f,10*f,100*f,1000*f,10000*f];
    
15. %中心频率与下限频率的比值
    
16. oc6=2^(1/6);
    
17. nc=length(fc);
    
18. n=length(x);
    
19. nfft=2^nextpow2(n);
    
20. a=fft(x,nfft);
    
21. for j=1:nc
    
22.     fl=fc(j)/oc6;
    
23.     fu=fc(j)*oc6;
    
24.     nl=round(fl*nfft/sf+1);
    
25.     nu=round(fu*nfft/sf+1);
    
26.     if fu>sf/2
    
27.         m=j-1;break
    
28.     end
    
29.     b=zeros(1,nfft);
    
30.     b(nl:nu)=a(nl:nu);
    
31.     b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    
32.     c=ifft(b,nfft);
    
33.     yc(j)=sqrt(var(real(c)));
    
34. end
    
35. figure
    
36. subplot(2,1,1);
    
37. t=0:1/sf:(n-1)/sf;
    
38. plot(t,x);
    
39. xlabel('时间(s)');
    
40. ylabel('加速度(g)');
    
41. grid on;
    
42. subplot(2,1,2);
    
43. plot(fc(1:m),yc(1:m));
    
44. xlabel('频率(Hz)');
    
45. ylabel('有效值');
    
46. xlim([0 100])
    
47. grid on;
    

复制代码

GuoCL

学习了，谢谢。

wwhblue

fc最大8*10000，sf=200 ??? 不符合采用定理吧！？