Функция filter

@any5957 · Регистрация: 15.06.2015

Author24 — интернет-сервис помощи студентам

Здравствуйте.
В MATLAB есть встроенная функция:

Matlab M

1	y = filter(a, b, x);

Данная функция применяет фильтр к сигналу х.
Подскажите, как реализовать данную функцию filter? То есть необходимо написать свою собственную функцию, которая по функционалу была бы аналогична встроенной функции filter в MATLAB.
Входные параметры:

Matlab M

а = [0.0688 0.2750 0.4125 0.2750 0.0688], 
b = [1.0000 -0.4630 0.9703 -0.4596 0.2598],
t = 0 : 0.1 : 100;
x = sin(t/2) + sin(3*t);

Я нашла описание этой команды:

Код

 a(1)y(n) = b(1) * x(n) + b(2) * x(n - 1) + ... + b(nb + 1) * x(n - nb) 
                                           - a(2) * y(n - 1) - ... - a(na + 1) * y(n - na),

при этом входной параметр b = [b(1) b(2) ... b(nb + 1)], а параметр a = [a(1) a(2) ... ... a(na+1)].

но как это реализовать - не получается.
Спасибо.

Написала фильтр. Фильтрует, но фильтрует не как функция filter.

Matlab M

а = [0.0688 0.2750 0.4125 0.2750 0.0688], 
b = [1.0000 -0.4630 0.9703 -0.4596 0.2598],
x = sin(t/2) + sin(3*t);
 
y = zeros(1, length(x));
    y(1) = b(1)*x(1);
    for k = 2:length(x) 
        t = b(1) * x(k);
        for n= 2:length(a)
            if(k - n > 1)&&(k - n < length(a))
               y(k) = y(k) + a(n)*x(k - n) - b(n)*y(k - n);
            end            
        end 
        y(k) = y(k) + t; 
    end

В чем ошибка?

Поведение моей функции фильтрации около нуля очень странное.
Идентичности между функцией filter и моей - нет.
Переписала код, как по формуле (поменяла местами a и b в коде):

Matlab M

y = zeros(1, length(x));
    y(1) = b(1)*x(1);
    for k = 2:length(x) 
        t = b(1) * x(k);
        for n= 2:length(a)
            if(k - n > 1)&&(k - n < length(a))
               y(k) = y(k) + b(n)*x(k - n) - a(n)*y(k - n);
            end            
        end 
        y(k) = y(k) + t; 
    end

Функция filter:

Моя функция фильтрации того же сигнала:

Подскажите, где ошибка?

@Annza · 12.12.2017, 17:55

В интернетах вроде уже как написали код функции filter, попробуйте сравнить своё https://www.mathworks.com/matl... ode-filter

@any5957 · 12.12.2017, 18:10 **[ТС]**

Annza, смотрела. Там как-то и проверки на вылет из индексов массива нет и параметр z присутствует, который мне не нужен.

@Annza · 12.12.2017, 18:18

Ну в окрестности 0 переходный процесс, можно отбрасывать первые length(b)/2 отсчетов. Непонятно, почему для первого отсчета выхода особенная формула:

Matlab M

1	y(1) = b(1)*x(1);

Потом, зачем отдельно прибавлять текущий вход, умноженный на b(1):

Matlab M

1	y(k) = y(k) + t;

почему нельзя индексацию по n вести от 1.

В данном примере количество коэффициентов a и b совпадает, но если эти вектора будут разной длины, то код тоже слетит...

@any5957 · 12.12.2017, 18:31 **[ТС]**

Annza, Как учесть тот факт, что, если длины a и b разные, то код бы сохранил работоспособность? Если учесть Ваши замечания, то будет выводится следующее:

Matlab M

function y = myFilter(b, a, x)
    y = zeros(1, length(x));    
    for k = 2:length(x)   
        y(k) = b(1) * x(k);
        for n= 1:length(b)
            if(k - n > 1)&&(k - n < length(b))
               y(k) = y(k) + b(n)*x(k - n) - a(n)*y(k - n);
            end            
        end         
    end    
end

@Annza · 12.12.2017, 18:50

Сделайте два внутренних цикла, один по индексам коэффициентов числителя передаточной функции b, другой по - a и добавляйте прямые и обратные связи по отдельности к сигналу. В остальном код соответствует приведенной в 1 сообщении формуле

@any5957 · 12.12.2017, 18:59 **[ТС]**

Annza, кстати, внешний цикл по k от 1?
Если стандартной функцией фильтровать, то сигнал от -1.5 до 1.5, а с помощью моего фильтра от -2 до 2. По-моему, это важное отличие. Почему такое происходит?
Получилось следующее:

Matlab M

function y = myFilter(b, a, x)
    y = zeros(1, length(x));    
    for k = 1:length(x)   
        y(k) = b(1) * x(k);
        for n = 1:length(b)
            if(k - n > 1)&&(k - n < length(b))
               y(k) = y(k) + b(n)*x(k - n);
            end            
        end  
        for n = 1:length(a)
            if(k - n > 1)&&(k - n < length(b))
               y(k) = y(k) - a(n)*y(k - n);
            end            
        end 
    end    
end

Не фильтрованный сигнал также находится в интервале от -2 до 2. И получается, что сигнал в итоге и не фильтруется вовсе.
Если применить стандартную функцию, то сигнал имеет такой вид:

@Krasme · 12.12.2017, 19:47

у меня такой не получается...
видимо, другой код выкладываете

@any5957 · 12.12.2017, 19:53 **[ТС]**

Krasme, Дико извиняюсь. Массивы a и b имеют другие значения. Вот, со всеми входными данными:

Matlab M

а = [0.0000595 0.0002381 0.0003572 0.0002381 0.0000595], 
b = [1.0000 -3.7087 5.2313 -3.3234 0.8020],
x = sin(t/2) + sin(3*t);
y = zeros(1, length(x));    
    for k = 1:length(x)   
        y(k) = b(1) * x(k);
        for n = 1:length(b)
            if(k - n > 1)&&(k - n < length(b))
               y(k) = y(k) + b(n)*x(k - n);
            end            
        end  
        for n = 1:length(a)
            if(k - n > 1)&&(k - n < length(b))
               y(k) = y(k) - a(n)*y(k - n);
            end            
        end 
    end

Krasme, если я правильно поняла, то у массива а - это умножение значений на 10 в минус третьей степени.

Krasme, но даже при измененных входных параметрах, общая картина с графиками остается той же. И следствие из этого - фильтр не рабочий.

@any5957 · 12.12.2017, 21:14 **[ТС]**

Сравнила исходный массив сигнала и массив сигнала после применения моего фильтра. Они абсолютно идентичны (кроме первых 10 значений). Значит, с кодом точно что-то не порядок...В чем может быть ошибка?

Добавлено через 10 минут
Подозреваю, что ошибка где-то в условии if.

@Krasme · 12.12.2017, 21:21

вектор а другой, чтобы разные размеры (a,b) были.

Кликните здесь для просмотра всего текста

Matlab M

function lklklk
clc; clear all; clf;
a = [0.0000595 0.0000595 0.0000595 0.0002381 0.0003572 0.0002381 0.0000595];
b = [1.0000 -3.7087 5.2313 -3.3234 0.8020];
t = 0 : 0.1 : 25;
x = sin(t/2) + sin(3*t);
y = filter(a, b, x);
z=myFilter(b,a,x);
plot(t,x,'g',t,y,'.b',t,z,'or')
legend('Orig','filter','MyFilt',0)
end
 
function y = myFilter(a, b, x)
na=length(a);
nb=length(b);
y = zeros(1,length(x));
for n=max(na,nb):length(x);
    y(n)=b(1)*x(n);
    NB=0; NA=0;
    for k=2:nb
        NB=NB+b(k)*x(n-k+1);
    end
    for k=2:na
        NA=NA+a(k)*y(n-k+1);
    end    
    y(n)=y(n)+NB-NA;
     y(n)=y(n)*a(1);
end
end
 
%   a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
%                         - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

@any5957 · 12.12.2017, 21:29 **[ТС]**

Krasme, ого. Спасибо!
Два вопроса:
Почему в цикле по n мы идем от максимальной длины векторов a, b?

Matlab M

1	for n=max(na,nb):length(x);

Почему мы в коэффициенте у х прибавляем еще единицу?

Matlab M

1	NB=NB+b(k)*x(n-k+1);

@Krasme · 12.12.2017, 21:35

Сообщение от any5957

коэффициенте у х прибавляем еще единицу?

смотрите логику в цикле

Matlab M

1	b(2)*x(n-1)

Сообщение от any5957

цикле по n мы идем от максимальной длины векторов a, b?

надо учитывать отступ по массиву

Matlab M

    x(n-k+1);

если взять n меньше, чем na или nb, то будет ошибка индекса массива x

@any5957 · 12.12.2017, 21:41 **[ТС]**

Krasme, поняла. Спасибо

@any5957 · 12.12.2017, 21:45 **[ТС]**

Krasme, фото результата работы нашего фильтра. Он такой и должен получиться (умножение на 10 в минус 7)?

Встроенный фильтр:

@Krasme · 12.12.2017, 21:47

Кликните здесь для просмотра всего текста

@any5957 · 12.12.2017, 21:48 **[ТС]**

Krasme, это для Ваших входных данных.
А для моих что-то вот так выводит...

@Krasme · 12.12.2017, 21:50

а вы повторили код, который я выложила?

@any5957 · 12.12.2017, 21:51 **[ТС]**

Krasme, да, конечно. Я им заменила код в своей функции. Только входные данные остались свои и отрисовка графика.

@Krasme · 12.12.2017, 21:52

ну давайте свой новый код полностью... посмотрю

@Krasme 6830 / 4890 / 2065 Регистрация: 02.02.2014 Сообщений: 13,048
	12.12.2017, 21:47	16
	Кликните здесь для просмотра всего текста 0

@Annza 79 / 61 / 25 Регистрация: 07.04.2013 Сообщений: 204
	12.12.2017, 17:55	2
	В интернетах вроде уже как написали код функции filter, попробуйте сравнить своё https://www.mathworks.com/matl... ode-filter 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 18:10 [ТС]	3
	Annza, смотрела. Там как-то и проверки на вылет из индексов массива нет и параметр z присутствует, который мне не нужен. 0

@Annza 79 / 61 / 25 Регистрация: 07.04.2013 Сообщений: 204
	12.12.2017, 18:50	6
	Сделайте два внутренних цикла, один по индексам коэффициентов числителя передаточной функции b, другой по - a и добавляйте прямые и обратные связи по отдельности к сигналу. В остальном код соответствует приведенной в 1 сообщении формуле 0

@Krasme 6830 / 4890 / 2065 Регистрация: 02.02.2014 Сообщений: 13,048
	12.12.2017, 19:47	8
	у меня такой не получается... видимо, другой код выкладываете 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 21:14 [ТС]	10
	Сравнила исходный массив сигнала и массив сигнала после применения моего фильтра. Они абсолютно идентичны (кроме первых 10 значений). Значит, с кодом точно что-то не порядок...В чем может быть ошибка? Добавлено через 10 минут Подозреваю, что ошибка где-то в условии if. 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 21:41 [ТС]	14
	Krasme, поняла. Спасибо 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 21:45 [ТС]	15
	Krasme, фото результата работы нашего фильтра. Он такой и должен получиться (умножение на 10 в минус 7)? Встроенный фильтр: 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 21:48 [ТС]	17
	Krasme, это для Ваших входных данных. А для моих что-то вот так выводит... 0

@Krasme 6830 / 4890 / 2065 Регистрация: 02.02.2014 Сообщений: 13,048
	12.12.2017, 21:50	18
	а вы повторили код, который я выложила? 0

@any5957 0 / 0 / 0 Регистрация: 15.06.2015 Сообщений: 64
	12.12.2017, 21:51 [ТС]	19
	Krasme, да, конечно. Я им заменила код в своей функции. Только входные данные остались свои и отрисовка графика. 0

@Krasme 6830 / 4890 / 2065 Регистрация: 02.02.2014 Сообщений: 13,048
	12.12.2017, 21:52	20
	ну давайте свой новый код полностью... посмотрю 0

Функция filter

Решение