请教功率频谱如何转成1/3 octave 频谱
我有一笔加速规撷取数据单位为电压值经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=;
fu=; Fc=;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fl为频宽下限,Fu为频宽上线,Fc为中心频率
?????????????
谢谢:@D
[ 本帖最后由 ChaChing 于 2010-4-20 13:48 编辑 ] 不懂"1/3 octave band与 加速度RMS(dB)的频谱",先搜搜论坛吧,加速度、速度、位移积分的问题论坛有的。建议发帖之前读读本版的置顶帖子，下次再越(遇)到这样的帖子我就删了，念你刚注册这个帖子就算了
[ 本帖最后由 无水1324 于 2007-8-3 08:37 编辑 ] 1/3倍频程的程序在MATHWORKS上有：
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectType=file&objectId=69
FFT后再进行1/3倍频程分析，在王济和胡晓编 “MATLAB在振动信号处理中的应用” （中国水利水电出版社）一书中有一节用介绍1/3倍频程分析，它是在FFT之后用1/3倍频程滤波器对信号进行分析处理，求出1/3倍频程滤波器输出的均方根值，并提供了MATLAB程序，楼主可参考。
回复 #3 songzy41 的帖子
呵呵，觉得1/3倍频没有3分频直观容易理解 感谢响应 看了MATHWORKS上对于octave的范例,其filtband主要是从25Hz开始,而我想求的是较低频的人体振动量测1~80Hz与ISO2631规范比较,不知有无CPB的程序范例?如果要用MATHWORKS的octave程序包，可参照它的方法，把采样频率、频带等参数改为适合自已的参数，然后使用。可能楼主更适合采用我提到的王济那书中的程序，其中1/3倍频程的中心频率从1Hz开始。 感谢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=;
fc=; 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(); grid on;
subplot(2,2,4); plot(freq,magxl1,'r');
xlabel('频率(Hz)'); ylabel('功率'); axis(); 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 编辑 ] 问一个很弱的问题：为什么要将功率频谱转成1/3 octave 频谱，是他有什么特殊的用途和意义吗？
回复 #9 无水1324 的帖子
标题想说明的意思有错误,其实是想做1/3 octave spectrum,用意是要与ISO2631的规范做比较,其横轴是1/3 octave center frequency,纵轴是RMS(dB)值,所以不晓得我上面程序是否有问题,请各位老师给一点建议 看了楼主的程序，也看了原书中的程序，但感觉原书中的这语句似乎有错：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');
便会出错。
感谢songzy41回应
由于excel档案太大,数据改用txt贴上
yc(j)=sqrt(var(real( b(1:n))));
yc(j)=sqrt(var(real( c(1:n))));
这两句的语法不清楚差异在何处 能否稍作解释, 作RMS?
另外下面建议部分已改正
感谢指导 从上面程序中可看出，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）运算。 RMS应该是norm(x)/sqrt(length(x))吧.
yc(j)=sqrt(var(real( c(1:n)))); 这只是标准差吧~~ 请教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 编辑 ] 因为原始数据长只是n，而多余的nff-n个点是由于在FFT时用nfft来变换， 在n个数据后补了nff-n个0值。所以在IFFT后只有前n个数据与原始数据相对应。
[ 本帖最后由 ChaChing 于 2010-5-5 23:44 编辑 ]
页:
[1]
2