новая попытка понять это.
Я пытаюсь понять, что делает функция bandpass
matlab.
Итак, копаясь в m-файлах, объект фильтра создается с помощью:
opts=signal.internal.filteringfcns.parseAndValidateInputs(x,'bandpass',var);
opts = designFilter(opts);
Затем фильтрация выполняется с помощью:
y = signal.internal.filteringfcns.filterData(x,opts);
Затем, переходя к функции filterData, фактическая фильтрация выполняется с помощью:
y = filtfilt(opts.FilterObject,x);
потому что это БИХ-фильтр. Теперь, если я загляну в функцию filtfilt
, у нее нет возможности работать с объектом фильтра, возможны только представления b-a
или sos-g
. Представление b-a
действительно плохое, поэтому я попробовал представление sos-g
, и результаты не совпадают.
Я построил три сигнала, отфильтрованные с помощью: filterData
, filtfilt
с объектом фильтра и filtfilt
с представлением sos-g
. Только первые два равны.
Далее идет PSD каждого сигнала: filtData: объект фильтра filtfilt: filtfilt sos-g: Обратите внимание на большую PSD на нижних частотах фильтра sos-g
.
Итак, что делает Matlab по-другому, когда объект фильтра используется в функции filtfilt
?
Вот пример (скачать data.txt):
load('data.txt')
srate=64;
freqrange=[0.4 3.5];
x=data(60*64:64*90);
var{1}=freqrange; var{2}=srate;
%% append at the beginning and at the end m=numel(x);
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];
%% design filter
opts=signal.internal.filteringfcns.parseAndValidateInputs(x,'bandpass',var);
opts = designFilter(opts);
%% default filtering
xx = signal.internal.filteringfcns.filterData(x,opts);
x_fil=xx(NR+1:end-NR);
%% what default filtering is supposed to do
xx = filtfilt(opts.FilterObject,x);
x_fil2=xx(NR+1:end-NR);
[z,p,k]=opts.FilterObject.zpk;
[sos, g] = zp2sos(z, p, k);
%% default filtering with sos-g method
xx = filtfilt(sos,g,x);
x_fil3=xx(NR+1:end-NR);
f=figure;
plot([x_fil x_fil2 x_fil3])
title('default vs ff_D vs ff_{sos}')
legend('filterData','filtfilt_D','filtfilt_{sos}')
grid('on')