关于1/3倍频程的讨论和请教

hotman007

本人刚接触信号处理方面的知识，遇到的问题求解一组检测信号的是1/3倍频程。因为主要关心振源产生的振动对精密仪器的影响，因此主要关心的频率段为1～80Hz范围，也就是相对较低的部分。
此前认真搜索了本论坛这方面的知识，遇到了几个问题，还希望各位高人能够予以指点。
首先是在gogo885的帖子，
http://forum.vibunion.com/thread-49623-1-1.html
其中程序的核心部分也根据songzy41老师的建议，部分程序如下：
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);
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(c(1:n))));  <---疑问
end

主要的疑问是，上面的c(1:n)是否是c(1:nfft)呢？？

此外，songzy41老师还推荐1/3倍频程的程序在MATHWORKS上给的函数,经简要改造后如下：
******************************************
function [B,A] = oct3dsgn(Fc,Fs,N);
if (nargin > 3) | (nargin < 2) %nargin变量的个数
  error('Invalide number of arguments.');
end
if (nargin == 2)
  N = 3;
end
if (Fc > 0.88*(Fs/2))
  error('Design not possible. Check frequencies.');
end
  
% Design Butterworth 2Nth-order one-third-octave filter
pi = 3.14159265358979;
f1 = Fc/(2^(1/6)); %带通的起始频率
f2 = Fc*(2^(1/6)); %带通的截止频率
Qr = Fc/(f2-f1);
Qd = (pi/2/N)/(sin(pi/2/N))*Qr;
alpha = (1 + sqrt(1+4*Qd^2))/2/Qd;
W1 = Fc/(Fs/2)/alpha;
W2 = Fc/(Fs/2)*alpha;
[B,A] = butter(N,[W1,W2]);
********************************************
function [p,f] = oct3bank(x);
pi = 3.14159265358979;
Fs = 200;     % Sampling Frequency
N = 3;      % Order of analysis filters.
ff = ((2^(1/3)).^[0:19]);  % Exact center freq.  
p = zeros(1,20);
% Design filters and compute RMS powers in 1/3-oct. bands
for i = 1:20
   [B,A] = oct3dsgn(ff(i),Fs,N);
   y = filter(B,A,x);              %<-----滤波器
   p(i) = sqrt(var(y));        %得到rms
end
*****************************************
　　但是另外一个问题是，采用以上两种方法得到的结果，却总是和仪器自带软件的处理结果有一定差异，特别是在低频阶段差异较明显。仪器的说明上说是采用FIR滤波器，因此，在这里请教各位老师前辈：
采用滤波器对1/3倍频程的影响表现在那些方面？
建议采用哪种滤波器能适合这种较低频率的1/3倍频程求解？？？

[ 本帖最后由 hotman007 于 2008-4-17 16:24 编辑 ]

songzy41

我曾有几个帖子介绍1/3倍频程滤波器，主要是介绍1/3倍频程滤波器的设计和使用方法。在http://forum.vibunion.com/forum/viewthread.php?tid=49623&extra=&page=1
中的程序主要取自王济和胡晓编 “MATLAB在振动信号处理中的应用”一书，在该方法中主要是利用FFT和IFFT，来获得1/3倍频程的滤波。 它仅是一种方法，在实际仪器设备上可能早已不采用了。
对楼主的问题“主要的疑问是，上面的c(1:n)是否是c(1:nfft)呢？？”应取c(1:n)。数据序列长n，但n不一定是2的幂次，在程序中把n扩展到2的幂次nfft，即在FFT时通过补0使数据序列长nfft，所以在滤波IFFT后，为了和原 数据长n对应，所以取c(1:n）。

在MATHWORK中的方法是用IIR的Butterworth滤波器设计成1/3倍频程滤波器。

　　楼主认为“采用以上两种方法得到的结果，却总是和仪器自带软件的处理结果有一定差异”，这和采用什么方法及怎么处理有关。滤波器用什么样的滤波器，IIR或FIR，即使用IIR，是用什么类型的，用FIR，又是用什么方法，用多少阶数;又在处理上每帧用多长(即多长求一次均方值)，等等这些因素都影响1/3倍频程滤波的结果。所以仪器和自带软件的处理结果有一定差异也是预料之中的。

hotman007

非常感谢song老师的耐心指导和回答，大部分疑惑已得以解决。
正如song老师所说，[ 滤波器用什么样的滤波器，IIR或FIR，即使用IIR，是用什么类型的，用FIR，又是用什么方法，用多少阶数；又在处理上每帧用多长（即多长求一次均方值），等等这些因素都影响1/3倍频程滤波的结果 ]。小弟对同一条信号采MATHWORK的3阶IIR的Butterworth滤波器在1～8Hz范围内的结果明显比仪器自带的FIR滤波器的结果小，因此向song老师以及各位在这方面有经验的老师请教：
为能得到较低频率(1~100Hz)的倍频程正确结果，应该选用哪种滤波器并如何设定其合理的参数呢？？希望大家不吝赐教，再次向song老师及各位关心该话题的朋友表示感谢！

[ 本帖最后由 hotman007 于 2008-4-18 11:40 编辑 ]

shashishan

songzy41 发表于 2008-4-17 18:19

                               
登录/注册后可看大图

我曾有几个帖子介绍1/3倍频程滤波器，主要是介绍1/3倍频程滤波器的设计和使用方法。在http://www.chinavib. ...

song老师你好，关于倍频程的mathwork方法有些疑问。设计的IIR滤波器的方法是：
f1 = Fc/(2^(1/6)); %带通的起始频率
f2 = Fc*(2^(1/6)); %带通的截止频率
Qr = Fc/(f2-f1);
Qd = (pi/2/N)/(sin(pi/2/N))*Qr;
alpha = (1 + sqrt(1+4*Qd^2))/2/Qd;
W1 = Fc/(Fs/2)/alpha;
W2 = Fc/(Fs/2)*alpha;
[B,A] = butter(N,[W1,W2]);
请问这里计算Qd，alpha起什么作用。同时，这里设计的是模拟滤波器，而输入信号是离散的，是不是应该设计成数字滤波器？希望您能给与解答。