和模拟,其中“复制模拟”会产生更好的解决方案。
这有点漏掉了重点。并不是人们很关心匹配或复制“模拟”,而是数字 IIR 滤波器具有一些非常好的和有用的特性。
例如在音频中,IIR 滤波器非常常见,我每天都使用 Butterworths。
最直接的原因只是过滤器长度。当采样率为 48 kHz 时,要在 40 Hz 下做任何有意义的事情,FIR 滤波器需要有数千个抽头长度。在某些情况下,您可以使用多速率过滤器来非常有效地实现这一点,但我不确定这是否普遍适用。另一种可能的选择是基于 FFT 的过滤器(重叠相加等),但这在延迟和效率之间提出了一个重要的权衡。
如果你坚持使用线性相位滤波器,延迟很快就会变得非常大,并且你会失去因果关系。从技术上讲,您可以通过简单地反转单位圆外的零点来将任何线性相位 FIR 转换为最小相位 FIR,但这在高阶数字上是一个尴尬的过程。
人耳对单耳相位失真(例如最小相位)也相当不敏感,只要其原因。然而,它对预振铃非常敏感,因此线性相位滤波器在感知上可能存在问题。
一个更“哲学的原因”是人类听觉或多或少地使用对数频率轴,而 IIR 滤波器更适合于此。例如,8kHz 的 IIR 倍频程带通滤波器看起来与 125Hz 的滤波器完全相同(中心频率除外)。FIR 滤波器将大不相同。
IIR 非常适合分频:奇阶巴特沃斯低通和高通滤波器总和为平坦的频率响应(尽管它通常是全通的)。对于偶数订单,您可以使用 Linkwitz Riley 过滤器(基于 Butterworths)。您可以改用线性相位 FIR,但同样会遇到延迟和因果关系问题。
像实时可调高通和低通滤波器这样的操作也很容易:动态设计巴特沃斯滤波器或简单地列出极点位置并进行插值非常容易。过滤器顺序始终保持不变,并且保证您具有相同的斜率和形状。
换句话说:打开你家中的任何扬声器,你会听到很多巴特沃思和朋友在行动:-)
示例 1
这只是一个非常简单的高通,您可以在任何花园品种有源扬声器中使用。以 48 kHz 采样的 40 Hz。我使用firls()
. 匹配是手动完成的,所以它们看起来一样。我需要在 FIR 滤波器上进行大约 6000 个抽头才能使其达到预期效果
% IIR
[z,p,k] = butter(6,40*2/48000,'high');
%FIR
h=firls(6000,[0 10 65 fs/2]*2/fs,[0 0 1 1])';
传递函数看起来像这样。
对我来说,FIR 在延迟、内存占用方面似乎非常逊色。一个 6 阶 BW 只有 13 个系数和 6-8 个状态变量。您可以通过利用零都在z=1.
示例 2:
这是对前一个示例的扩展,它使高通的截止频率可以实时调整。这会创建一个“滑动高通”滤波器,这在密封的智能扬声器中很常见。我试图让需求变得现实,就像你在真实产品中发现的那样。
Butterworth 允许非常高效的实现,没有可听伪影、没有固有延迟、非常低的内存占用和非常低的 CPU 消耗。
%% Script to implement a sliding high pass filter, that can be adjusted on the fly
% This type of "sliding highpass" is typically used in smart speakers with a
% "closed box" topology to control trasnducer excursion at high output
% volume
%
% Requirements:
% - sample rate: 48 kHz
% - frame size: 128 samples
% - cutoff frequency varies from 40 Hz to 100 Hz
% - highpass slope: 36 dB/octave. In other words level at half the cutoff
% frequency should be <= -36 dB
% - Level at cutoff frequency not less than -3dB
% - passband ripple: < .01 dB above 200 Hz (where spectral perception is
% more sensitive)
% - cutoff frequency is updated once per frame
% - no additional latency
% - smooth updating filter, no pops or clicks
%
% Implementation
% - 6th order butterworth
% - pole location and filter gain as a function of frequency are
% approximated as a polynomial fit. At these low frequencies, a linear
% fit (1rst order polynomial) works perfectly.
% - The filter coefficients are calcuated once per frame based on the
% current cutoff frequency
% - Filter is implemented as cascaded second order sections in Direct Form
% I. For Direct Form I, the state variable are simply the inputs and
% outputs, so updateing the filter deosn't create a discontinuity in the
% state variables.
%
% Test
% - input signal is a 50 Hz sine wave, low frequency sine waves are very
% sensitive to artifacts
% - cutoff frequency input is a mixture of step function, up and down sweeps
% and uniformly distributed random numbers between 40Hz and 100Hz
%
%% Table up the poles and the gain of the filters, perform polynomial fit
ord = 6;
fs = 48000;
fr = (40:.1:100)';
nfr = length(fr);
p0 = cell(nfr,1); % pole locations
k0 = zeros(nfr,1); % filter gains
for i = 1:nfr
[z,p,k] = butter(ord,2*fr(i)/fs,'high');
p0{i} = p(imag(p)>0);
k0(i) = k;
end
p0 = [p0{:}].'; % convert cell array to regular array
%% do a simple linear for poles and gains
% we fit the real part and the gain in (1-x) since they are very close to 1
ppPoles = cell(3,2);
for i = 1:3
ppPoles{i,1} = polyfit(fr,1-real(p0(:,i)),1); % 1 minus real part
ppPoles{i,2} = polyfit(fr,imag(p0(:,i)),1); % imginary part
end
% and the gains in 1-k
ppGain = polyfit(fr,1-k0,2);
%% test the polyfit, calculate the resulting transfer functions at a few
% frequencies
frTest = 40:10:100;
nfrTest = length(frTest);
nFFT = 16384*2;
d0 = zeros(nFFT,1); d0(1) = 1;
fy = zeros(nFFT,nfrTest);
sos = zeros(ord/2,6);
sos(:,1) = 1; sos(:,2) = -2; sos(:,3) = 1;
for i = 1:nfrTest
f = frTest(i);
for ip = 1:3
x = 1-polyval(ppPoles{ip,1},f); % real part
y = polyval(ppPoles{ip,2},f); % imginary part
sos(ip,4:6) = [1 -2*x x^2+y^2];
end
k = 1-polyval(ppGain,f);
fy(:,i) = fft(k*sosfilt(sos,d0));
end
% this checks out fine
%% Real time part, preperation
frameSize = 128; % frame size
frameRate = fs/frameSize; % number of frames per second
% let's do 10 seconds but an integer number of frames
nx = 10*fs;
numFrames = floor(nx/frameSize);
nx = numFrames*frameSize;
% build a test signal: 50 hz sine wave
xin = sqrt(.5)*sin(2*pi*(0:nx-1)'*50/fs); % 50 Hz sine wave at -3 dB
% build control frequency signal, one frequency per frame
frInput = 40*ones(numFrames,1);
% let's do a bit of rectangular switching plus some random stuff
t = (1:frameRate);
frInput(frameRate+t) = 100;
frInput(3*frameRate+t) = 70;
% up and down sweep
frInput(5*frameRate+t) = linspace(40,100,frameRate);
frInput(6*frameRate+t) = flip(linspace(40,100,frameRate));
% some random numbers for good measure
frInput(7*frameRate+1:end) = 40+60*rand(3*frameRate,1);
%
% in order to smooth the very aprupt frequency transitions in the test
% vector, we smooth the frequency input with a time constant to 30ms
timeConstant = 0.03;
frSmooth = exp(-1./(timeConstant*frameRate));
frCurrent =40;
%% Real time over all frames
% we implement this as driect form I so the states are always guaranteed to
% be continous
signalState = zeros(2,4); % filter state, we need total of 8
freqState = frInput(1); % cutoff frequency states
xout = 0*xin; % initialize output
t = 1:frameSize; % time vector
t0 = 0; % current time
yy = [xout xout xout];
for iFrame = 1:numFrames
% get frequency input and apply smoothing
frCurrent = frSmooth*frCurrent + (1-frSmooth)*frInput(iFrame);
% calculate filter gain, grab input and scale it
k = 1-polyval(ppGain,frCurrent);
y = k*xin(t0+t);
% over all biquads
x1 = signalState(1,1); % grab input state for the first biquad
x2 = signalState(2,1);
for iPole = 1:3
% calculate the filter coeefficents
pReal = 1-polyval(ppPoles{iPole,1},frCurrent); % real part
pImag = polyval(ppPoles{iPole,2},frCurrent); % imaginary part
a1 = -2*pReal; % filter coefficient "a1"
a2 = pReal*pReal+pImag*pImag; % biquad coefficient "a2"
% grab the output state
y1 = signalState(1,iPole+1);
y2 = signalState(2,iPole+1);
% inner loop. Here efficieny is the most important
for i = t
x0 = y(i); % get input sample
y0 = x0-2*x1+x2-a1*y1-a2*y2; % DF1 Butterworth stage
% update state
x2 = x1; x1 = x0;
y2 = y1; y1 = y0;
y(i) = y0; % write output
end % end sample loop
% save the output state of this stage
signalState(1,iPole) = x1;
signalState(2,iPole) = x2;
% grab the input state for the next stage
x1 = signalState(1,iPole+1);
x2 = signalState(2,iPole+1);
% now write the output state
signalState(1,iPole+1) = y1;
signalState(2,iPole+1) = y2;
yy(t0+t,iPole) = y;
end % end poles/stages
% write output
xout(t0+t) = y;
t0 = t0 + frameSize;
end % end frame
plot(xout);