Плавающее среднее., поиск коэффициентов. Опции

[_sv_]

Добрый день всем.
Реализую простейший фильтр "скользящее среднее".
output = input + FILTER_CONST * (output – input);

Требуется адаптивно выбирать FILTER_CONST в процессе работы.
Не могу найти способ вычисления оной.

Мне, собственно, задают допустимое время накопления и порог.
Надо как-то вычислить FILTER_CONST.

Спасибо.

[_pv]

это не скользящее среднее, это БИХ фильтр НЧ первого порядка, с 20дБ/дек начиная с частоты среза (-3дБ) которая равна Fs*FILTER_CONST (Fs - частота дискретизации).

[demiurg_spb]

http://en.wikipedia.org/wiki/Low-pass_filter

[Corner]

Практически цифровая RC-цепочка...

это не скользящее среднее, это БИХ фильтр НЧ первого порядка, с 20дБ/дек начиная с частоты среза (-3дБ) которая равна Fs*FILTER_CONST (Fs - частота дискретизации).

Проще говоря если FILTER_CONST >0,5 схема превращается в генератор...

[andyp]

Вам бы надо

output(t) = alpha*input + (1-alpha)*output(t-1) , 0 < alpha < 1.

Если надо скользящее среднее на N отсчетов, то alpha~=2/(N+1) - последнее по wiki написано и наверное стоит это проверить в чем-либо типа матлаба, сравнив АЧХ получившегося фильтра и фильтра скользящего среднего c ИХ вида 1/N*ones(1,N).

Ваша FILTER_CONST = 1- alpha

Сообщение отредактировал andyp - Sep 6 2013, 18:49

[thermit]

Топикстартер, если хочешь внятных ответов, задавай вопрос внятно. Твой фильтр будет делать все. тока не усреднять.

y(n)=y(n-1)+(x(n)-y(n-1))*k
это распространенный фнч бих 1-го порядка для всяких усреднений.

k=1 - (1-L)^(1/T/sf)
L-уровень
sf-частота дискретизации
T-время достижения уровня



Сообщение отредактировал thermit - Sep 7 2013, 05:19

[_sv_]

Сразу - извиняюсь за долгое молчание.
Был в отъезде. Спасибо всем за быстрые ответы.

Теперь по делу: В первом посте досадная опечатка:
Следует читать "output = input + FILTER_CONST * (input - output)"
И это не "мой фильтр". В смысле, что я его не выдумывал.
Это первое попавшееся мне решение задачи "усреднение за время".
Мне всего-то надо усреднять данные от датчика света за некие интервалы времени.

y(n)=y(n-1)+(x(n)-y(n-1))*k

Да. Именно это и было реализовано (после устранения ошибки).

k=1 - (1-L)^(1/T/sf)

Вот. Именно эту формулу и не мог найти. Спасибо.

Вам бы надо
output(t) = alpha*input + (1-alpha)*output(t-1) , 0 < alpha < 1.

Тоже интересно. Насколько я понимаю, то это практически то же самое, что и у "thermit".
Надо бы разобраться как в этом случае считается коэффициент alpha.

Сижу с матлабом. Экспериментирую.

Еще раз - спасибо всем.
Мне кажется, что будут еще вопросы.

[thermit]

[quote name='_sv_' date='Sep 7 2013, 10:51' post='1190907']


Тоже интересно. Насколько я понимаю, то это практически то же самое, что и у "thermit".
[/qoute]

да. это то же самое. альфа=k

[andyp]

Следует читать "output = input + FILTER_CONST * (input - output)"

Да все у Вас нормально было:

out(k) = in(k) + FILTER_CONST*(out(k-1) - in(k));

после перегруппировки слагаемых получаем

out(k) = (1 - FILTER_CONST) in(k) + FILTER_CONST*out(k-1)

если FILTER_CONST = 1 - alpha, то это та форма, что приведена в wiki. Таким образом, имеем два ФНЧ один - КИХ (скользящее среднее на N отсчетов) с АЧХ вида sin(Nx)/sin(x), второй - БИХ с АЧХ вида abs(alpha*exp(j*omega)/(exp(j*omega) + (1 - alpha))). Последнее выражение можно упростить, но мне сейчас лень. Чтобы оба фильтра интегрировали одинаково по крайней мере белый шум, нужно, чтобы были равны их шумовые полосы в полосе интегрируемого сигнала. Шумовая полоса - это интеграл АЧХ, нормированный на максимальный коэффициент передачи фильтра в некотором диапазоне частот. Отсюда и получается связь между alpha и N.

[_sv_]

При более детальном изучении оказалось, что:
1) существуют два подхода к вопросу:
output = input + FILTER_CONST * (output – input);
output = input + FILTER_CONST * (input - output);
При этом и тот и другой реализуют разные вещи.
(input - output) - фильтр
(output – input) - какой-то там интегратор.

В связи с чем:

Твой фильтр будет делать все. тока не усреднять.

Ваше утверждение поставлено под сомнение.
Изучаю предмет дальше.
Действительно - пока что полученные результаты ну никак не "плавающее среднее".

[thermit]

1 можете детально изучать/подвергать сомнению итд итп сколько угодно и что угодно. только ровно до тех пор, пока вы не научитесь внятно записывать разностные уравнения, ваши (именно ваши) фильтры будут делать все что угодно.

2 к вопросу усреднения в контексте вашей задачи существует ровно 1 подход: выборочное среднее (скользящее среднее). Все остальные подходы - более общая нч-фильтрация со своими заморочками.


3 если вы думаете заменить однородный ких (скользящее среднее по N отсчетам) эквивалентным бих-ом 1-го порядка то это врядли. так что не тратьте зря время. ибо

4 строгий эквивалентный бих выглядит так:
y(n)=x(n)/N-x(n-N)/N+y(n-1)

5 изучать эти вопросы по википедии не надо.

[_sv_]

если вы думаете заменить однородный ких ... эквивалентным бих-ом 1-го порядка то это врядли.

Вы знаете... Я несколько раз прочитал свои посты.
Не могу найти кто и откуда взял, что меня интересуют БИХ фильтры, или какие-то там замены.
Гляньте, пожалуйста и Вы.

Все, что меня действительно интересует это "скользящее среднее".
У меня все более или менее получается, кроме достаточно точного вычисления коэффициентов.
Мне хочется, что бы пользователь задал задержку в диапазоне, скажем, от 5 до 3600 секунд.
И что бы после этого фильтр на идеальных данных (в Excel или Matlab) давал ошибку в, скажем, 5%.
Пока что в диапазоне 270...420 секунд у меня ошибка почти 30%.
В других интервалах все хорошо.

... получаем ...
out(k) = (1 - FILTER_CONST) in(k) + FILTER_CONST*out(k-1)

Да. Действительно получаем. Уже добрался до этого.

Отсюда и получается связь между alpha и N.

Ваши посты продвигают меня по дебрям DSP.
Совершенно новая для меня тема.

Я получаю данные от датчика раз в секунду.
Диапазон значений - от 0 до 32767.
Мне надо достаточно надежно отличить состояния "день", "ночь", "искусственное освещение".
Уровни для каждого состояния будут определены чуток позже.
Когда появится возможность установить оборудование в реальный корпус.
Задача настолько тривиальна, что мне по простоте душевной, казалось, что ответ я получу в течении 10 минут.

[thermit]

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

вопщем мне нечего вам больше сказать.

[_sv_]

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

Я, безусловно, не специалист в DSP. И в постах и в терминах наверняка есть ошибки.
Однако, название "Скользящее Среднее" было взято из:
wiki_exponential_moving_average
Мне казалось (да и сейчас кажется), что название соответствует действительности.

вопщем мне нечего вам больше сказать.

В общем и на том спасибо.

Разбирательство с выкладками andyp привело меня на DSPLIB
Где я, похоже, нашел ответ на мои вопросы.
Проверяю в матлабе.
Спасибо ему за это.

Сообщение отредактировал _sv_ - Sep 7 2013, 21:08

[andyp]

вот небольшой скрипт для матлаб - октав иллюстрирующий выравнивание шумовых полос для значения alpha = 2/M (т.е. фильтр скользящего среднего на М отсчетов и экспоненциальный ФНЧ при этом значении alpha имеют примерно равные шумовые полосы.) Все становится не очень хорошо при малом M = 5 (ошибка до 20%), но уже для М = 36 она не превышает одного процента. Для малых M можно еще пошаманить с подбором alpha.

Код

Ms = [5, 36, 360, 3600];  %количество отсчетов, участвующее в усреднении
show_plot = 1;
for  ii = 1 : length(Ms)
    M = Ms(ii);
    freq_ananlysis_len = 4096;
    f = linspace(0,freq_ananlysis_len-1,freq_ananlysis_len) * pi/(freq_ananlysis_len-1) * min(1, 20/M);
    s1 = abs(freqz(1/M*ones(1,M),1,f));
    alpha = 2/(M);
    s2 =  abs(freqz(alpha,[1 , -(1 -alpha)],f));
    ss(ii) = sum(s1.^2)/sum(s2.^2);
    aa = ss(ii)
    if show_plot
        mm = 10*log10(s1);
        ix = find(mm < -30); mm(ix) = -30;    
        plot(mm,'r');
        hold on;
        mm = 10*log10(s2);
        ix = find(mm < -30); mm(ix) = -30;    
        plot(mm,'g');
        hold off
        grid on
        pause(5)
    end
end

err = abs(ss - 1)

[_sv_]

вот небольшой скрипт для матлаб...
Все ... не очень хорошо при малом M = 5 (ошибка до 20%), но уже для М = 36 она не превышает одного процента.

У меня результаты сильно отличаются от Ваших.
Почему - буду разбираться дальше.
Правда, ошибки у меня уже в пределах, для меня допустимых.
Продолжу завтра.

Спасибо.

[Corner]

Использую функцию выход=(выход+вход)/2. Если ее каскадировать, то получается хороший сплайн для прямоугольных сигналов.
Правильный эквивалент все-таки: выход=выход+(вход-выход)/2.
Или в записи _sv_: output = output + FILTER_CONST * (input - output), где FILTER_CONST <= 0,5 иначе схема - генератор. Это интегрирующая цепочка, только вот частота среза далеко не Fs*FILTER_CONST, так как при FILTER_CONST = 0,5 она равна Fs*0,113.

Сообщение отредактировал Corner - Sep 9 2013, 06:42

[_sv_]

Правильный эквивалент ... output = output + FILTER_CONST * (input - output),

Та вроде - да.

интегрирующая цепочка ... частота среза ... равна Fs*0,113.

И это тоже, похоже, верно.


Математические эксперименты мною пока прекращены. Обнаружились проблемы на плате.
Установлены два двух ваттных ИК светодида. Жарят так, что забивают мне датчик.
Это при том, что приняты все (как мне казалось) меры для предотвращения этого эффекта.
Фильтры, механические пергородки и т.д.

К теме вернусь чуток позже.

[_sv_]

К теме вернусь чуток позже.

Вернулись к вопросу.
Оказалась проблема не совсем там, где искали.
Фильтры и все остальное работают точно так, как должны.

Вылезла проблема перегрева платы при включенной ИК подсветке.
Грелся датчик вместе с платой. Давал медленно ползущие ввех значения.
Я считал, что ошибка в матемалике фильтра.

Спасибо всем участовавшим и отдельное andyp