Простой полосовой фильтр

Вот мигалочка из 4 полосовых «фильтров» на тиньке.



Как работает эта штука?

Теория


Представим себе такую конструкцию.

Маятник

Грузик висит на конце пружинки, прикреплённой к неподвижному предмету. Если слегка оттянуть шарик влево или вправо, он начнёт колебаться с некоторым периодом T0. Как мы знаем из школьного курса физики из гугла, период этот равняется:

T_0 = 2 %pi sqrt{m over k}

Где m — масса грузика, а k — жёсткость пружинки.

А собственная частота колебаний, соответственно, будет равна 1 / T0. Думаю, ни для кого не секрет, что если подталкивать грузик ровно с этой частотой, амплитуда колебаний будет максимальна. Причём, чем слабее подталкивать, тем сильнее амплитуда колебаний будет зависеть от частоты подталкивания.

Получается такой механический фильтрик. Попробуем сделать его модель.

Из физики мы знаем ещё несколько формул. Сила, с которой пружинка действует на грузик, равняется:

F = -kx,

где k — жёсткость пружинки, а x — растяжение пружинки. Эта сила придаёт грузику ускорение, равное:

a = F/m,

где m — масса грузика. Отсюда получается, что:

a = -xk/m,

Заменим k/m на T. Из самой первой формулы получаем, что:

T = 4π2F2

Теперь к самой модели. Нам понадобятся 2 переменных — скорость движения грузика v и положение грузика x. Также понадобится шаг моделирования t.

На каждом шаге скорость грузика будет меняться по такой формуле:

v1 = v + at,

или

v1 = v — xTt,

а положение по формуле

x1 = x + vt.

В нашей системе t — в секундах, F — в герцах, а x — в попугаях. Стоит перейти к милисекундам и килогерцам, чтобы не получать страшные числа. Все формулы при этом волшебным образом остаются теми же. :)

Реализация


Для начала, попробуем просто получить колебательный процесс.

#include <avr/io.h>
#include <avr/interrupt.h>

// T = (2*pi*Fkhz)^2
#define T_1000HZ		39

volatile int v = 0, x = 500;

// 16 kHz
ISR(TIMER1_COMPA_vect)
{
	v -= (x * T_1000HZ) >> 4;
	x += v >> 4;
}

int main()
{
	// Timer1 CTC 16 kHz
	TCCR1B = (1<<WGM12)|(1<<CS10);
	OCR1A = F_CPU/16000-1;
	TIMSK |= 1<<OCIE1A;

	sei();

	DDRB = 0xff;
	DDRD = 0xff;

	while(1) {
		PORTB = (x>>3) + 128;
		PORTD = (v>>5) + 128;
	}

	return 0;
}


Моделируем колебания грузика с шагом в 1/16 мс (умножение на 1/16 мс заменяем на сдвиг вправо на 4 бита). Данные отправляем в порты B и D чтобы понимать что у нас там твориться. :)

Схемка

К портам подключаем ЦАПы, а к ним осциллограф.

Осциллограммка

Вроде работает. Положение и скорость меняются по синусоидам, сдвинутым на 90 градусов, прям как ток и напряжение в колебательном контуре. Частота ровно 1 кГц, как и планировалось. :)

Теперь добавим входной сигнал.

#include <avr/io.h>
#include <avr/interrupt.h>

// T = (2*pi*Fkhz)^2
#define T_1000HZ		39

volatile int v,x;

ISR(TIMER1_COMPA_vect)
{
	v -= (x * T_1000HZ) >> 4;
	x += v >> 4;
}

ISR(ADC_vect)
{
	int sample = ADC - 512;
	x = (x * 31 + sample) >> 5;
}

int main()
{
	// Timer1 CTC 16 kHz
	TCCR1B = (1<<WGM12)|(1<<CS10);
	OCR1A = F_CPU/16000-1;
	TIMSK |= 1<<OCIE1A;

	// ADC enable (Fcpu/32)
	ADCSRA = (1<<ADEN)|(1<<ADSC)|(1<<ADFR)|(1<<ADIE)|(1<<ADPS2)|(1<<ADPS0);

	sei();

	DDRB = 0xff;

	while(1) {
		PORTB = (x>>2) + 128;
	}

	return 0;
}


Здесь мы немного «подталкиваем грузик» сэмплами из АЦП. Вот так подключено.

Схемка

А вот и результат.

Работа фильтра

Синий график — сигнал с генератора, жёлтый — с ЦАПа.

Тест


Вот такую мигалочку я и собрал чтобы проверить идею. :)

Схемка

Платко.
Платка

Вот окончательный код.

#include <avr/io.h>
#include <avr/interrupt.h>
#include <util/delay.h>
#include <stdlib.h>


#define T1		16L		// ~600 Hz
#define T2		64L		// ~1200 Hz
#define T3		128L	// ~1800 Hz
#define T4		512L	// ~3600 Hz

#define led1_on(ms) { led1 = ms; PORTB &= ~(1<<PB3); }
#define led2_on(ms) { led2 = ms; PORTB &= ~(1<<PB4); }
#define led3_on(ms) { led3 = ms; PORTB &= ~(1<<PB5); }
#define led4_on(ms) { led4 = ms; PORTB &= ~(1<<PB6); }

volatile uint8_t led1, led2, led3, led4;

volatile int x1,v1;
volatile int x2,v2;
volatile int x3,v3;
volatile int x4,v4;

// ovf ~2ms
ISR(TIMER0_OVF0_vect)
{
	if(led1 > 2) led1 -= 2;
	else PORTB |= 1<<PB3;

	if(led2 > 2) led2 -= 2;
	else PORTB |= 1<<PB4;

	if(led3 > 2) led3 -= 2;
	else PORTB |= 1<<PB5;

	if(led4 > 2) led4 -= 2;
	else PORTB |= 1<<PB6;
}

// ctc 16 kHz
ISR(TIMER1_CMPA_vect)
{
	v1 -= (x1 * T1) >> 4;
	x1 += v1 >> 4;

	v2 -= (x2 * T2) >> 4;
	x2 += v2 >> 4;

	v3 -= (x3 * T3) >> 4;
	x3 += v3 >> 4;

	v4 -= (x4 * T4) >> 4;
	x4 += v4 >> 4;
}

ISR(ADC_vect)
{
	int sample = ADC - 512;

	x1 = (x1 * 63 + sample) / 64;
	x2 = (x2 * 63 + sample) / 64;
	x3 = (x3 * 63 + sample) / 64;
	x4 = (x4 * 63 + sample) / 64;
}

int main()
{
	// I/O
	PORTA = 0xff & ~((1<<PA3)|(1<<PA0));
	PORTB = 0xff;
	DDRB = (1<<PB3)|(1<<PB4)|(1<<PB5)|(1<<PB6);

	// Timer0 Fcpu/64 (ovf ~2ms)
	TCCR0 = (1<<CS01)|(1<<CS00);
	TIMSK |= 1<<TOIE0;

	// Timer1 CTC 16 kHz
	TCCR1B = (1<<CTC1)|(1<<CS11);
	OCR1C = F_CPU/2/16000-1;
	TIMSK |= 1<<OCIE1A;

	// ADC Fcpu/32
	ADCSR = (1<<ADEN)|(1<<ADSC)|(1<<ADFR)|(1<<ADIE)|(1<<ADPS2)|(1<<ADPS0);

	sei();

	while(1) {
		if(x1 > 250) led1_on(20);
		if(x2 > 250) led2_on(20);
		if(x3 > 250) led3_on(20);
		if(x4 > 250) led4_on(20);
	}

	return 0;
}


На коллекторе транзистора должна быть примерно половина питания (если что, подстраивается номиналом R1). Тинька затактована от внутренней RC-цепочки на 8 МГц.

Хм, ну вот и всё…

Файлы в топике: filter-tests.zip, blinker.zip

Комментарии (36)

RSS свернуть / развернуть
Вбрось лучше сюда:
we.easyelectronics.ru/blog/Soft/
0
Да ладно. Может кто-нить напишет как по людски тоже самое сделать.
0
Автор, я сейчас отложил небольшой кирпич, когда включил видео работы сей мигалки. Звук у меня на всю стоял. :-)
+1
собери лучше ультразвуковую ванночку на базе 10-литрового тазика)
хочется нечто поздоровее того, что продают в китаях (промывать детали от мотора), но опять же 200 долларов на большой формат бюджета нету.
0
В ультразвуковой ванночке есть две проблемы.
1) УЗ излучатель ватт на 100-200 (для 10-литрового тазика-то)
2) Качающий его генератор (не то чтоб сильно сложно, но на выходе там приличный дроссель и несколько киловольт амплитуды)
И если второе собрать можно, то первое на каждом углу не валяется.
0
А если несколько излучателей поставить? )
0
Тогда нужно купить несколько излучателей! Опять же, там уже может потребоваться и расчет, чтобы они внезапно не оказались работающими в противофазе.
0
Наверно это всё же проще чем найти таблетку на 100-200вт.:)
Хм, правда сама ванночка, как я понимаю, требует рассчёта.
0
Честно говоря, я и на 50-то ватт таблеток не видел. Только в ванне.
0
что-нибудь из этого? :-)
www.voltmaster.ru/cgi-bin/qwery_i.pl?code=%F3%EB%FC%F2%F0%E0
0
Эммм что например?
0
там есть раздел с датчиками (пр/пер) — подумал, может оно)))
а еще вот здесь чо-то обсуждают похожее
flyback.org.ru/viewtopic.php?t=1178&postdays=0&postorder=asc&start=0
0
Ну вон что-то из раздела «излучатели звука» подойти может. Это те, которые примерно по 1.2кр.
0
Забавненько. А чем не устроили формулы из ЦОС? Хотя, ЕМНИП, там получается все примерно то же (что и неудивительно).
0
  • avatar
  • Vga
  • 06 июня 2011, 02:20
Забавно видеть на плате штыри, к которым припаяны провода. Ладно, если было бы какое-то устройство, которое делал кто-то другой. Но сам же поставил штыри, и потом передумал, и решил припаяться к ним… ;)
Порадовал тег «ненужное». Это из разряда «чтобы продать что-нибудь ненужное, нужно сначала купить что-нибудь ненужное»?
0
А меня больша порадовала связка тегов «интересное+ненужное» :)
0
Ооо ты тоже юзаешь православные кт315 =)
0
Хм, классная статья. Вывод формул порадовал, спасибо.
0
Прикольная статья, интересный алгоритм. А как с быстродействием, по сравнению, например, с тем же БПФ?
ЗЫ: а что за трек играет в видеоролике?
0
Ну вот так на си, с оптимизацией по скорости выполнения — сравнимо в принципе с реализацией от chan'а на асме. Однако памяти жрёт поменьше.
Играет lifelover — lethargy
+1
Вчера намучался… полез разбираться в коде — долго циклился на формуле T=(2*pi*Fkhz)^2 — пока сам не вывел. Т — не самое удачное обозначение, вводит в ступор. Ну и костлява реализация мигания светодиодов в посднем куске кода:) Второй пример промоделировать в протеусе не удалось…
Попутно возник вопрос: а что, если входной сигнал будет сдвинут по фазе, тогда ведь модель маятника нормально не сработает, или я не прав?
0
Неа, маятник «подстраивается» под фазу.
+1
подскажите, какой потолок у данного алгоритма? ну скажем, 15-20 кГц сигнал можно будет таким методом обработать? и еще один вопрос, на схеме нарисован модуль от In до ADC, где транзистор и прочее — он какую роль играет? (извиняюсь за вопросы, просто начинающий еще совсем)
0
Частота моделирования, думаю, раза в 4 хотя бы должна превышать частоту выделяемого сигнала. Т.е., скажем, 64 кГц. При тактовой частоте 20000 кГц получается ~125 тактов на шаг. Так что, думаю, можно упихнуться. Другое дело, что АЦП у AVR работает примерно до 16 кГц. А желательно, чтобы частота дискретезации в пару раз превашала частоту сигнала.

Транзистор это тупо усилитель. Можно его выкинуть, или заменить чем-нить поприличнее. И low-pass фильтр, вот он очень желателен.
0
Ой, при тактовой частоте 8000 кГц, от внутренней RC-цепочки.
0
ну а если внешний кварц повесить? По-шустрее. Может получиться или нет? Вот именно на 20 кГц сигнал.
0
Что-то проект под мегу 8 не компилится, а 16 тини в протеусе нету, чтобы проверить :(
0
В тиньке по другому называются вектора прерываний. Под межку надо немного подправить исходник.
ну а если внешний кварц повесить?
Почему бы и нет? Просто я юзал тиньку 26L, которая до 8 МГц.
0
можете помочь под мегу 8 или 16 подправить исходник? хочется промоделировать как в работе будет.
0
Хм, мона. Помоему быстрее написать заново. Там же одна формула)
0
да я в сях пока еще как свинья в апельсинах. Пытаюсь разобраться, но пока идет туговато. Вот с паскалем — милое дело ;)
0
Вообще-то эту штуку лучше на асме написать) Куда меньше ресурсов будет жрать)
На Си приходится включать оптимизацию по скорости.
0
Ну нету разницы между паскалем и С в плане работы с железом. Точно так же меняешь имена ISR-ок и т.п. на соответствующие меге. Базовый же синтаксис С можно выучить за вечер.
0
Обычный усилок. Довольно неплохо их тема рассмотрена в цикле статей «Усилитель» от mzw.
0
Давно вот хочу поинтересоваться. А как реагирует фильтр на сигнал вида sin(w1*t)+sin(w2*t), где w1 — собственная частота «колебаний маятника», а w2 — какая-то «левая». А при w2=n+*w1 (n=2,3,4...)
0
  • avatar
  • Deer
  • 26 августа 2011, 21:08
Всё, что не совпадает с собственной частотой, оказывает минимальное влияние)
int sample = ADC - 512;
x = (x * 31 + sample) >> 5;

Причём, чем меньше влияние тут нового сэмпла, тем меньше влияют посторонние частоты (впрочем, и тем медленнее оно раскачивается на собственной частоте).
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.