使用过滤器对象和 b,a 系数的结果非常不同
Very different results using filter object and b,a coefficients
我正在比较使用 MATLAB filter
对象的数字滤波输出和使用 tf
函数的 b-a 系数,它们确实不同。如果有人能指出我正确的方向,我将不胜感激。这里有我的结果和测试脚本(下载data.txt
):
srate=64;
freqrange=[0.4 3.5];
file='data.txt';
load(file);
file=split(file,'.');
file=file{1};
data=eval(file);
st=srate*60;
ed=srate*60*2;
data=data(st:ed);%1 minute data
m=numel(data);
x=data;
R=0.1;%10% of signal
Nr=50;
NR=min(round(m*R),Nr);%At most 50 points
x1=2*x(1)-flipud(x(2:NR+1));%maintain continuity in level and slope
x2=2*x(end)-flipud(x(end-NR:end-1));
x=[x1;x;x2];
[xx,dbp]=bandpass(x,freqrange,srate);%same result if use filter(dbp,xx)
data_fil=xx(NR+1:end-NR);
[b,a]=dbp.tf;
xx=filter(b,a,x);
data_fil_ba=xx(NR+1:end-NR);
f=figure;
s(1)=subplot(2,1,1,'Parent',f);
plot(s(1),[data data_fil])
title(s(1),'dbp')
s(2)=subplot(2,1,2,'Parent',f);
plot(s(2),[data data_fil_ba])
title(s(2),'ba')
linkaxes(s,'x');
这里是我得到的结果:
这是由 [b,a]
滤波器表示中的数值不稳定引起的,来自 Matlab docs:
Numerical Instability of Transfer Function Syntax
In general, use the [z,p,k]
syntax to design IIR filters. To analyze or implement your filter, you can then use the [z,p,k]
output with zp2sos
. If you design the filter using the [b,a]
syntax, you might encounter numerical problems. These problems are due to round-off errors and can occur for n as low as 4.
用 fvtool
查看你的过滤器,你会看到第一个图(digitalFilter
对象)是平滑的,而第二个图是用传递函数表示生成的非常平滑在较低的频率范围内不稳定。该滤波器的幅度甚至超过 0 dB,这意味着该滤波器不稳定并将放大这些频率。这就是为什么你的结果爆炸了,它变成了 ~10^92。
解决办法就是不使用传递函数([b,a]
)表示,而是使用digitalFilter
或[z,p,k]
表示。有关为什么要小心切换 filter/model 表示的更多信息,请参阅 here。
代码
srate=64;
freqrange=[0.4 3.5];
[~,dbp]=bandpass(rand(100,1),freqrange,srate);
[b,a]=dbp.tf;
fvtool(dbp)
fvtool(b,a)
我正在比较使用 MATLAB filter
对象的数字滤波输出和使用 tf
函数的 b-a 系数,它们确实不同。如果有人能指出我正确的方向,我将不胜感激。这里有我的结果和测试脚本(下载data.txt
):
srate=64;
freqrange=[0.4 3.5];
file='data.txt';
load(file);
file=split(file,'.');
file=file{1};
data=eval(file);
st=srate*60;
ed=srate*60*2;
data=data(st:ed);%1 minute data
m=numel(data);
x=data;
R=0.1;%10% of signal
Nr=50;
NR=min(round(m*R),Nr);%At most 50 points
x1=2*x(1)-flipud(x(2:NR+1));%maintain continuity in level and slope
x2=2*x(end)-flipud(x(end-NR:end-1));
x=[x1;x;x2];
[xx,dbp]=bandpass(x,freqrange,srate);%same result if use filter(dbp,xx)
data_fil=xx(NR+1:end-NR);
[b,a]=dbp.tf;
xx=filter(b,a,x);
data_fil_ba=xx(NR+1:end-NR);
f=figure;
s(1)=subplot(2,1,1,'Parent',f);
plot(s(1),[data data_fil])
title(s(1),'dbp')
s(2)=subplot(2,1,2,'Parent',f);
plot(s(2),[data data_fil_ba])
title(s(2),'ba')
linkaxes(s,'x');
这里是我得到的结果:
这是由 [b,a]
滤波器表示中的数值不稳定引起的,来自 Matlab docs:
Numerical Instability of Transfer Function Syntax
In general, use the
[z,p,k]
syntax to design IIR filters. To analyze or implement your filter, you can then use the[z,p,k]
output withzp2sos
. If you design the filter using the[b,a]
syntax, you might encounter numerical problems. These problems are due to round-off errors and can occur for n as low as 4.
用 fvtool
查看你的过滤器,你会看到第一个图(digitalFilter
对象)是平滑的,而第二个图是用传递函数表示生成的非常平滑在较低的频率范围内不稳定。该滤波器的幅度甚至超过 0 dB,这意味着该滤波器不稳定并将放大这些频率。这就是为什么你的结果爆炸了,它变成了 ~10^92。
解决办法就是不使用传递函数([b,a]
)表示,而是使用digitalFilter
或[z,p,k]
表示。有关为什么要小心切换 filter/model 表示的更多信息,请参阅 here。
代码
srate=64;
freqrange=[0.4 3.5];
[~,dbp]=bandpass(rand(100,1),freqrange,srate);
[b,a]=dbp.tf;
fvtool(dbp)
fvtool(b,a)