сравнение отфильтрованных данных3: filtData (полоса пропускания) и filtfilt с объектом фильтра и filfilt с матрицей sos

новая попытка понять это.

Я пытаюсь понять, что делает функция 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: filtData объект фильтра filtfilt: filtfilt filter object filtfilt sos-g: 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')

person kurokirasama    schedule 20.11.2019    source источник
comment
вы пробовали более простой сигнал, который вы можете предсказать результаты?   -  person Gideon Kogan    schedule 24.11.2019
comment
... еще не об этом. Однако дизайн фильтра зависит от данных, поэтому у меня может быть другое поведение...   -  person kurokirasama    schedule 25.11.2019