如何在MatLAB/Octave中创建1/3倍频程滤波器

Question

我必须设计1/3-八倍频程滤波器在MatLAB (或交替在八度)。我读过这文档的文章，我尝试过使用fdesign.octave-design二重奏，但是这种方法允许在25小时左右开始为中频创建带滤波器。我需要创建频率范围从0.5Hz到100 to的滤波器。

我知道如何计算这些滤波器的中频(fm)以及相应的上、下频段，所以我尝试使用butter或cheby1函数或使用fdatool创建它们。虽然在fdatool中创建/显示的过滤器看起来不错，但当我使用它们来使用filer函数过滤chirp信号时，我可以在下面的图像中看到如下所示：

​

顶部图是用我的自定义过滤器过滤的chirp，底部图显示相同的chirp过滤与相同的顺序、类型和中频的过滤器，但是使用fdesign.octave-design duo创建。我不知道如何摆脱这个信号增益在开始和这个实际上是一个最低的信号增益，它显示100赫兹的中频频率。越低的中频频率，我的滤波器是更坏的(实际上，低于2或3赫兹我的滤波器滤除一切)。

她的代码是我用来构建过滤器的代码：

Wn = [Fl(i) Fh(i)] / (Fs/2);
[ z, p, k ] = butter( N, Wn, 'bandpass' );
% [z,p,k] = cheby1( N, Rp, Wn, 'bandpass' );
[ sos, g ]=zp2sos( z, p, k );
f = dfilt.df2sos( sos, g );

有人能帮我设计出适合那个频率范围(0.5×100)Hz的1/3倍频程滤波器吗？

编辑:我找到了另一种方法来使用标准化的fdesign.octave方法来做我想做的事情。下面读我的答案。

Answer 1

只需手动计算出你想要的中心频率，然后计算上、下界。下面的演示使用简单的butterworth过滤器，但您可以将其与任何您喜欢的原型交换。

fs = 250;
bw = 1/3;

fMin = 0.5;
fMax = 100;

octs = log2(fMax/fMin);
bmax = ceil(octs/bw);

fc = fMin*2.^( (0:bmax) * bw ); % centre frequencies
fl = fc*2^(-bw/2); % lower cutoffs
fu = fc*2^(+bw/2); % upper cutoffs

numBands = length(fc);

b = cell(numBands,1);
a = cell(numBands,1);

figure
for nn = 1:length(fc)

    [b{nn},a{nn}] = butter(2, [fl(nn) fu(nn)]/(fs/2), 'bandpass');
    [h,f]=freqz(b{nn},a{nn},1024,fs);

    hold on;
    plot(f, 20*log10(abs(h)) );

end
set(gca, 'XScale', 'log')
ylim([-50 0])

这给了我。。。

​

​

然后，您可以使用tf2sos，如果您想要稳定性从高阶滤波器作为一个级联的二阶区段。

您可以使用下面的代码来测试chirp，看看没有什么奇怪的事情发生。。

dt=1/fs;
t = 0:dt:2000;      % 2 seconds
fo = 0.5; f1 = 120; 
x = chirp(t,fo,numel(t)*dt,f1,'logarithmic')';
y = zeros(numel(x), length(fc));
for nn = 1:length(fc)
    y(:,nn) = filter(b{nn},a{nn},x);
end

figure; plot(y)

​

​

Answer 2

虽然我接受了learnvst的回答(因为它在一开始对我有帮助)，但后来我找到了其他方法来创建正确的过滤器，这将符合ANSI S1.11 / IEC 61260，而且实际上非常容易。我的想法是基于数字滤波器的一种特殊特性，即它们可以用于不同采样率的信号，从而改变滤波器的实际通带范围。

我使用了标准化的fdesign.octave方法，但参数值不同。对于中频(25×100)赫兹范围，我做一切正常的方式，但范围从2.5Hz (高达24赫兹)，我已经指定了10倍更大的采样率和10倍的中频频率。对于0.5Hz及以上的频率(<2.5Hz)，我将这些值乘以100。因此，我获得了有效的滤波器(0.5×100)赫兹范围内的实际采样率，我的源信号。