请教功率频谱如何转成1/3 octave 频谱

gogo885

我有一笔加速规撷取数据单位为电压值经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 编辑 ]

花如月

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

[ 本帖最后由 无水1324 于 2007-8-3 08:37 编辑 ]

songzy41

1/3倍频程的程序在MATHWORKS上有：
http://www.mathworks.com/matlabc ... ile&objectId=69
FFT后再进行1/3倍频程分析，在王济和胡晓编 “MATLAB在振动信号处理中的应用” （中国水利水电出版社）一书中有一节用介绍1/3倍频程分析，它是在FFT之后用1/3倍频程滤波器对信号进行分析处理，求出1/3倍频程滤波器输出的均方根值，并提供了MATLAB程序，楼主可参考。

花如月

呵呵，觉得1/3倍频没有3分频直观容易理解

gogo885

感谢响应 看了MATHWORKS上对于octave的范例,其filtband主要是从25Hz开始,而我想求的是较低频的人体振动量测1~80Hz与ISO2631规范比较,不知有无CPB的程序范例?

songzy41

如果要用MATHWORKS的octave程序包，可参照它的方法，把采样频率、频带等参数改为适合自已的参数，然后使用。可能楼主更适合采用我提到的王济那书中的程序，其中1/3倍频程的中心频率从1Hz开始。

gogo885

感谢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 编辑 ]

无水1324

问一个很弱的问题：为什么要将功率频谱转成1/3 octave 频谱，是他有什么特殊的用途和意义吗？

gogo885

标题想说明的意思有错误,其实是想做1/3 octave spectrum,用意是要与ISO2631的规范做比较,其横轴是1/3 octave center frequency,纵轴是RMS(dB)值,所以不晓得我上面程序是否有问题,请各位老师给一点建议

songzy41

看了楼主的程序，也看了原书中的程序，但感觉原书中的这语句似乎有错：
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');
便会出错。

gogo885

data.txt (142.26 KB, 下载次数: 68)

2007-8-7 01:34 上传

点击文件名下载附件

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

yc(j)=sqrt(var(real( b(1:n))));
yc(j)=sqrt(var(real( c(1:n))));
这两句的语法不清楚差异在何处 能否稍作解释, 作RMS?
另外下面建议部分已改正
感谢指导

songzy41

从上面程序中可看出，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）运算。

zhj0231984

RMS应该是norm(x)/sqrt(length(x))吧.  
yc(j)=sqrt(var(real( c(1:n)))); 这只是标准差吧~~

hotman007

请教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 编辑 ]

songzy41

因为原始数据长只是n，而多余的nff-n个点是由于在FFT时用nfft来变换， 在n个数据后补了nff-n个0值。所以在IFFT后只有前n个数据与原始数据相对应。

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