声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 16060|回复: 24

[综合讨论] 请教功率频谱如何转成1/3 octave 频谱

[复制链接]
发表于 2007-8-3 01:04 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我有一笔加速规撷取数据单位为电压值经FFT转成功率频谱后, 要再化成1/3 octave band与 加速度RMS(dB)的频谱下面该如何作?
下面是目前程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
x=xlsread('data.xls','a1:a8192');
fs=1/0.002; N=length(x); freqStep=fs/N; freq=freqStep*(-(N-1)/2:(N-1)/2);
y=fft(y); y1=abs(y); magx1=fftshift(y1); T=16;
fl=[0.891,1.122,1.412,1.778,2.238,2.817,3.547,4.465,5.621,7.077,8.909,11.216,14.120,17.776,22.378,28.173,35.467,44.651,56.212,70.767,89.090,112.158,141.198,177.758,223.784];
fu=[1.122,1.413,1.779,2.240,2.819,3.550,4.469,5.626,7.082,8.916,11.225,14.131,17.790,22.396,28.195,35.495,44.686,56.256,70.823,89.160,112.246,141.310,177.898,223.961,281.950]; Fc=[1,1.25,1.6,2,2.5,3.15,4,5,6.3,8,10,12.5,16,20,25,31.5,40,50,63,80,100,125,160,200,250];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fl为频宽下限,  Fu为频宽上线,  Fc为中心频率
?????????????

谢谢:@D

[ 本帖最后由 ChaChing 于 2010-4-20 13:48 编辑 ]
回复
分享到:

使用道具 举报

发表于 2007-8-3 08:03 | 显示全部楼层
不懂"1/3 octave band与 加速度RMS(dB)的频谱",先搜搜论坛吧,加速度、速度、位移积分的问题论坛有的。建议发帖之前读读本版的置顶帖子,下次再越(遇)到这样的帖子我就删了,念你刚注册这个帖子就算了

[ 本帖最后由 无水1324 于 2007-8-3 08:37 编辑 ]
发表于 2007-8-3 09:11 | 显示全部楼层
1/3倍频程的程序在MATHWORKS上有:
http://www.mathworks.com/matlabc ... ile&objectId=69
FFT后再进行1/3倍频程分析,在王济和胡晓编 “MATLAB在振动信号处理中的应用” (中国水利水电出版社)一书中有一节用介绍1/3倍频程分析,它是在FFT之后用1/3倍频程滤波器对信号进行分析处理,求出1/3倍频程滤波器输出的均方根值,并提供了MATLAB程序,楼主可参考。

评分

1

查看全部评分

发表于 2007-8-3 10:00 | 显示全部楼层

回复 #3 songzy41 的帖子

呵呵,觉得1/3倍频没有3分频直观容易理解
 楼主| 发表于 2007-8-4 14:03 | 显示全部楼层
感谢响应 看了MATHWORKS上对于octave的范例,filtband主要是从25Hz开始,而我想求的是较低频的人体振动量测1~80HzISO2631规范比较,不知有无CPB的程序范例?
发表于 2007-8-4 15:53 | 显示全部楼层
如果要用MATHWORKS的octave程序包,可参照它的方法,把采样频率、频带等参数改为适合自已的参数,然后使用。可能楼主更适合采用我提到的王济那书中的程序,其中1/3倍频程的中心频率从1Hz开始。
 楼主| 发表于 2007-8-6 00:57 | 显示全部楼层
感谢songzy41推荐的“MATLAB在振动信号处理中的应用一书
我参照书上程序来编程 想请教一些观念厘清

程序:
x=xlsread('data.xls','a1:a8192'); %实际量测数据 单位电压(V)
x=(x-2.488)/1.005);   % 转成加速度值 单位(g)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%预处理
L=length(x); %t1=1:L;
amean=sum(x)/L; a=x-amean;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sf=500;
f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
fc=[f,10*f,100*f]; oc6=2^(1/6); nc=length(fc); n=length(x);
nfft=2^nextpow2(n);
freqStep=sf/n;
freq=freqStep*(-(n-1)/2:(n-1)/2);
a1= fft(a,nfft); a2=abs(a1);
a3=fft(a,nfft)*2/n;   %<--------幅值
a4=abs(a3);
magxl=fftshift(a4);
magxl1=magxl.^2; %<--------功率
for j=1:nc
fl=fc(j)/oc6;
fu=fc(j)*oc6;
nl=round(fl*nfft/sf+1);
nu=round(fu*nfft/sf+1);
if fu > sf/2, m=j-1;break; end
b=zeros(1,nfft);
b(nl:nu)=a1(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a1(nfft-nu+1:nfft-nl+1);

c=ifft(b,nfft);

yc(j)=sqrt(var(real(b(1:n))));  %<------
不懂??
yc(j)=20*log10(yc(j)/0.000001); %<-------
dB
end
subplot(2,2,1); t=0:1/sf:(n-1)/sf; plot (t,x);
xlabel('时间(s)'); ylabel('加速度(g)'); grid on;
subplot(2,2,3); plot(freq,magxl,'r');
xlabel('(Hz)'); ylabel('振幅'); axis([0,10,0,0.5]); grid on;
subplot(2,2,4); plot(freq,magxl1,'r');
xlabel('频率(Hz)'); ylabel('功率'); axis([0,10,0,0.1]); grid on;
subplot(2,2,2); plot(fc(1:m),yc(1:m),'g');
xlabel('频率(Hz)'); ylabel('RMS(dB)'); grid on;

想请教这样做法有问题吗?
还有分贝谱的部分是否正确?

补图

补图


[ 本帖最后由 ChaChing 于 2010-5-5 23:40 编辑 ]
发表于 2007-8-6 08:02 | 显示全部楼层
问一个很弱的问题:为什么要将功率频谱转成1/3 octave 频谱,是他有什么特殊的用途和意义吗?
 楼主| 发表于 2007-8-6 13:30 | 显示全部楼层

回复 #9 无水1324 的帖子

标题想说明的意思有错误,其实是想做1/3 octave spectrum,用意是要与ISO2631的规范做比较,其横轴是1/3 octave center frequency,纵轴是RMS(dB)值,所以不晓得我上面程序是否有问题,请各位老师给一点建议
发表于 2007-8-6 16:20 | 显示全部楼层
看了楼主的程序,也看了原书中的程序,但感觉原书中的这语句似乎有错:
yc(j)=sqrt(var(real(b(1:n))));  %<------不懂??
应改为
yc(j)=sqrt(var(real( c(1:n))));
另楼主的程序没有给出数据,最好把数据附上,便于对照。
freqStep=sf/n;
freq=freqStep*(-(n-1)/2:(n-1)/2);
改为
freqStep=sf/nfft;
freq=freqStep*(-(nfft-1)/2:(nfft-1)/2);
否则n≠nfft时
plot(freq,magxl,'r');
plot(freq,magxl1,'r');
便会出错。

评分

1

查看全部评分

 楼主| 发表于 2007-8-7 01:34 | 显示全部楼层
data.txt (142.26 KB, 下载次数: 68)


感谢songzy41回应
由于excel档案太大,数据改用txt贴上

yc(j)=sqrt(var(real( b(1:n))));
yc(j)=sqrt(var(real( c(1:n))));
这两句的语法不清楚差异在何处 能否稍作解释, RMS?
另外下面建议部分已改正
感谢指导
发表于 2007-8-7 09:34 | 显示全部楼层
从上面程序中可看出,b是信号在某一个1/3滤波器内的频谱值:
b(nl:nu)=a1(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a1(nfft-nu+1:nfft-nl+1);
它是在本身是复数值,如果取其实部只取了其频谱值的一部分。
c是由b经IFFT得到的,是信号经某一个1/3滤波器后输出的时域数据,它原应该为实数值,但由于计算误差的存在,使虚部也有值,但很小。对c只取实部即得到了滤波器的输出数据。var函数和sqrt函数便是对滤波器的输出作了方均根(RMS)运算。
发表于 2008-1-18 14:38 | 显示全部楼层
RMS应该是norm(x)/sqrt(length(x))吧.  
yc(j)=sqrt(var(real( c(1:n)))); 这只是标准差吧~~
发表于 2008-4-10 14:33 | 显示全部楼层
请教songzy41老师:
关于程序 yc(j)=sqrt(var(real( c(1:n)))) 一行,我个人认为也应该是c,但是应该是c(1:n)还是c(1:nfft)?? 因为 c=ifft(b,nfft)变换得到的是一个nfft个点的信号,而程序中为什么只截取到n呢??
急切盼望songzy41以及各位1/3倍频程的高手指点!

怎么没有回答?是我的问题问错了还是太幼稚了??:@Q

[ 本帖最后由 ChaChing 于 2010-5-5 23:41 编辑 ]
发表于 2008-4-17 20:13 | 显示全部楼层
因为原始数据长只是n,而多余的nff-n个点是由于在FFT时用nfft来变换, 在n个数据后补了nff-n个0值。所以在IFFT后只有前n个数据与原始数据相对应。

[ 本帖最后由 ChaChing 于 2010-5-5 23:44 编辑 ]

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-2 23:35 , Processed in 0.090082 second(s), 27 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表