Почему фильтр некорректно работает при частотах фильтруемой гармоники ниже 95 рад/с?

Question

presence @presence

Обработка сигналов

Почему фильтр некорректно работает при частотах фильтруемой гармоники ниже 95 рад/с?

Задача: реализовать НЧФ на основе ряда Фурье, используя сглаживающие множители Кайзера, частота пропускания 21рад/с, Т дискретизации 0.1с, ширина переходной полосы 1р/с, неравномерность в полосе задержки/пропускания 0.05, коэф. усиления 1.

function del=bess(k) % Вычисление значения модифицированной функции Бесселя первого рода нулевого порядка 
F=1;
s=1;
S=0;
i=1;
while (s>10^(-10)) 
    F=F*i;
    s=((k/2)^i/F)^2;
    S=S+s;
    i=i+1;
end
del=1+S;
end

function Aw=frequency_response(w,L,wl,dt,wk) % Вычисление отсчетов АЧХ
S=0;
for i=1:1:L
   S=S+wk(L-i+1)*cos(i*dt*w);
end
Aw=wl+2*S;
end

delta=0.05; % Ввод исходных данных
wp=1;
dt=0.1; 
wv=3.14/0.1;
dw=wv/74;
w1=7;
w2=95;%Не может быть выше 125.6 рад/с по теор. котельникова
R=((0.34^2-4*1.610*(19.6+20*log10(delta)))^(0.5)-0.34)/(2*1.610); % Вычисляем параметр R множителей Кайзера
L=fix(1+2*(R^2+pi^2)^(0.5)/(wp*dt)); % Определяеv число членов аппроксимирующего ряда Фурье
I0=bess(R);  % Вычисление сглаживающих множителей Кайзера
for k=1:1:L 
dk(k)=bess(R*(1-(k/L)^2)^(0.5))/I0;
end
for k=1:1:L % Вычисление отсчетов импульсной характеристики
wk(k)=dk(L-k+1)/((L-k+1)*pi)*sin((L-k+1)*((21+wp/2)*dt));
wk(2*L-k+1)=wk(k);
t=[1:1:(2*L)]; %Построение графика ИХ
plot(t,wk)
end
wl=(21+wp/2)/(pi/dt);
m=1;
for i=0:0.1:(pi/dt) %Вычисление отсчетов АЧХ
I(m)=i;
w(m)=frequency_response(i,L,wl,dt,wk);
f(m)=-L*dt*i; % Вычисление отсчетов ФЧХ
m=m+1;
end
figure; plot(I,w) % Построение графика АЧХ
grid on
set(gca, 'XTick',0:1:(pi/dt)+2)
set(gca, 'YTick',-0.2:0.05:1.2)
figure; plot(I,f) % Построение графича ФЧХ
grid on

q1=fix(w1/dw);
q2=fix(w2/dw);

pp=0:0.1:50; 
x=sin(w1*pp)+sin(w2*pp);
figure; plot(pp, x) % Построение входного полигармонического сигнала
grid on

pe=0:0.05:50;
x=w(q1)*sin(w1*dt*pe+f(q1))+w(q2)*sin(w2*dt*pe+f(q2));
figure; plot(pe, x) % Построение выходного сигнала
grid on

pu=0:0.05:50;
x=w(q1)*sin(w1*dt*pu)+w(q2)*sin(w2*dt*pu);
figure; plot(pu, x) % Построение выходного скорректированого сигнала
grid on

Так код реализован на матлабе, 2 функции bess.m и frequency_responce.m при частотах второй гармоники (не входящей в полосу пропускания) более 95 рад/с, но меньше 125 фильтр работает корректно и выявляет гармонику, входящую в полосу пропускания, но при более низких частотах (но выше 21 рад/с (полоса пропускания) сигнал выявляется некорректно. С чем это связанно?