[_sv_ 0]

Добрый день всем.

Реализую простейший фильтр "скользящее среднее".

output = input + FILTER_CONST * (output – input);

 

Требуется адаптивно выбирать FILTER_CONST в процессе работы.

Не могу найти способ вычисления оной.

 

Мне, собственно, задают допустимое время накопления и порог.

Надо как-то вычислить FILTER_CONST.

 

Спасибо.

 

 

 

[_pv 52]

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

[demiurg_spb 0]

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

[Corner 0]

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

 

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

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

[andyp 9]

Вам бы надо

 

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

Изменено 6 сентября, 2013 пользователем andyp

[thermit 1]

 

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

 

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

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

 

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

L-уровень

sf-частота дискретизации

T-время достижения уровня

 

 

Изменено 7 сентября, 2013 пользователем thermit

[_sv_ 0]

Сразу - извиняюсь за долгое молчание.

Был в отъезде. Спасибо всем за быстрые ответы.

 

Теперь по делу: В первом посте досадная опечатка:

Следует читать "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 1]

 

 

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

[/qoute]

 

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

 

 

 

 

 

 

 

[andyp 9]

Следует читать "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_ 0]

При более детальном изучении оказалось, что:

1) существуют два подхода к вопросу:

output = input + FILTER_CONST * (output – input);

output = input + FILTER_CONST * (input - output);

При этом и тот и другой реализуют разные вещи.

(input - output) - фильтр

(output – input) - какой-то там интегратор.

 

В связи с чем:

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

Ваше утверждение поставлено под сомнение.

Изучаю предмет дальше.

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

 

 

 

 

 

[thermit 1]

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

 

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

 

 

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

 

4 строгий эквивалентный бих выглядит так:

y(n)=x(n)/N-x(n-N)/N+y(n-1)

 

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

[_sv_ 0]

если вы думаете заменить однородный ких ... эквивалентным бих-ом 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]

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

 

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

 

 

[_sv_ 0]

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

Я, безусловно, не специалист в DSP. И в постах и в терминах наверняка есть ошибки.

Однако, название "Скользящее Среднее" было взято из:

wiki_exponential_moving_average

Мне казалось (да и сейчас кажется), что название соответствует действительности.

 

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

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

 

Разбирательство с выкладками andyp привело меня на DSPLIB

Где я, похоже, нашел ответ на мои вопросы.

Проверяю в матлабе.

Спасибо ему за это.

Изменено 7 сентября, 2013 пользователем _sv_

[andyp 9]

вот небольшой скрипт для матлаб - октав иллюстрирующий выравнивание шумовых полос для значения 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)