Шумодав для рацайки. Часть 3. Фильтры.

Зачем.


В предыдущей части мы определились, что входной сигнал будет оцифровываться с частотой 128кГц в то время как ШП будет работать с сигналом оцифрованным на 8кГц.
В обратную сторону наоборот — 8кГц в 128кГц.
Плюс этого в том что увеличенная частота дискретизации в 16 раз позволит нам, согласно этому, выиграть 2 бита. Т.е. мы, после передискретизации можем считать, что входной сигнал имеет размах 14 бит. Однако для этого необходимо. чтобы на входе действовал белый шум. Ну, шума у нас достаточно. Примерно такая же ситуация на выходе.
Дальше все плюсы кончаются…



Входной фильтр.

Вообще говоря, если посмотреть на запись сигнала, то можно обойтись и без фильтра — брать каждый 16й отсчет входного сигнала, ибо спектр выхода детектора ограничен примерно 3кГц. Однако, в этом случае мы не выиграем наши 2 бита. Поэтому будем, все-таки, применять фильтр.

Чтобы определиться с выбором фильтра надо рассмотреть какими ресурсами мы располагаем: 1. Имеем плавающую точку на борту. 2. Имеем много памяти…

Нам нужен фильтр, который бы имел значительное ослабление на частоте 4кГц и пропускал бы все от 0 до, как минимум 3х кГц.
В принципе можно применить и рекурсивный фильтр (IIR) на входе, но тогда мы столкнемся с проблемой нелинейности ФЧХ (в прицепе наплевать) и с предельным циклом (блуждающими нулями). Поэтому будем использовать длинный нерекурсивный фильтр (FIR).

Запускаем матлаб, fdatool и придумываем… Ну, мне понравился фильтр Чебышева, порядок 127, частота среза 3600Гц, ослабление внеполосных сигналов -60дБ. Смотрим что получилось:



Вроде ничего себе так получилось.

Пока что пропустим детали реализации и откуда мы получим 2 лишних бита. Это позже.
Ослабление на 4кГц составляет 10дБ. Достаточно ли этого? Для того, чтобы это понять, нужно сгенерить шум с полосой от 4 до 64кГц, пропустить через этот фильтр, понизить частоту дискретизации и посмотреть что получилось.

Исходный шум + синус на 1кГц, чтобы было с чем сравнивать уровень шума, выглядит так (спектр):


Фильтрованный сигнал уже с частотой дискретизации 8кГц выглядит так:


Видим, что просачивается примерно -50дБ внеполосных сигналов интегрально. Это радует. Кстати, на нижней картинке синусоида уже усилена на 12дБ — это я учел 2 бита.

Реализация входного фильтра.

Мы на вход получаем 12 битные целочисленные отсчеты входного сигнала @128кГц. На выходе хотим получить 14 бит @8кГц. Вообще говоря, поскольку мой ШП может работать с 16и битными отсчетами, то выход фильтра мы приведем к 16и битам, помня то, что младшие 2 бита шум.
КИХ фильтр это обычная свертка и ничего интересного здесь нет.

Собственно вот код фильтра:


#define OUTBUF_SIZE				(8*16)	// 128

#define FRAME_TIME_MS			16	// mS
#define SAMPLING_RATE			128	// kHz
#define FRAME_SIZE				(FRAME_TIME_MS*SAMPLING_RATE)

extern int16_t pf[PRE_FILTER_LEN];    // 128 длина

int16_t *
pre_filter(int16_t *data, int16_t *inbuf, int16_t *outbuf)
{
	int i, j, k;
	int16_t *b = inbuf + PRE_FILTER_LEN;
	
	for (i = 0; i <FRAME_SIZE; i++)
	{
		b[i] = data[i]  - (0x1000/2);
	}

	for (i = 0, k = 1; i < OUTBUF_SIZE; i++, k += 16)
	{
		int32_t	sum = 0;
		
		for (j = 0; j < PRE_FILTER_LEN; j++)
		{
			sum += (int32_t)pf[j] * b[k-j];	// 128 sums 
		}
		
		sum += 1 << (15-0);	// Round
		sum >>= 15+1;		// shift to q1.15
		outbuf[i] = sum;	// store: Sum will be in range -16386something ... +16385something!!! However it may overload on clicks!
	}

	dmamem2move(inbuf, b + FRAME_SIZE - PRE_FILTER_LEN, PRE_FILTER_LEN);
	
	return outbuf;
}


Для начала нам нужно избавиться от постоянного смещения — это первый цикл. При этом данные будем записывать в линию задержки. Поскольку памяти у нас много, то ее пока можно не экономить.

Для фильтрации мы будем использовать целочисленную арифметику — т.е. арифметику с фиксированной точкой (где-то здесь об этом уже писалось).

Конечно же есть и float, но с ним получается несколько кривенько… не знаю почему, лень было разбираться…

Коэффициенты фильтра из матлаба сохраняем в формате Q1.31 — 32 битные целые. При этом это все будут числа от -1 до примерно 1, а точка находится между 31м и 30м разрядом числа. Бит 31 знаковый.
Произведение 2х чисел с фиксированной точкой делается в 2 (или даже в 3) этапа — умножение, округление и приведение к нужной размерности. При этом результат умножения должен иметь размер равный сумме количества дробных бит первого и второго числа.
У наших коэффициентов из 31 и у данных 11. Т.е. промежуточный результат должен иметь как минимум 31+11 + знак == 43 бита. Это не очень хорошо, поэтому мы понизим разрядность коэффициентов до 20 бит — сдвинем все коэффициенты вправо на 12. Тогда, вроде бы, все нормально получается для одного числа. Однако приведение одного произведения к требуемой размерности с последующим суммированием может усилить шум «обработки».
Мы поступим по-другому — мы будем сначала умножать и складывать, а потом результат приводить к требуемой размерности. Но нам нужно просуммировать аж 128 раз! Это требует лишних 7 бит. Однако заметим, что самый большой коэффициент занимает 27 — 12 == 16 бит. Тогда 11+15 + 7 == 32 бита. Т.е. нам бы еще 1 бит надо найти. Но! Ситуация, когда требуется вся разрядная сетка возможна лишь в том случае, если на входе будет меандр полного размаха с частотой 4 кГц. Во всех остальных случаях нужно меньше. Поэтому будем считать, что переполнения не случится.

Теперь, собственно, что делается в коде. Во втором цикле для каждого выходного отсчета вычисляется произведение с накоплением 128и предыдущих входных отсчетов и коэффициентов фильтра. Затем результат приводится к формату Q1.15. Правда, мне тут удобно использовать не весь диапазон 16и битного числа, а лишь его половину. За сим сдвигать будем не на 15, а на 16.

Ну и потом нам надо сохранить последние 128 отсчетов входного сигнала в буфере для следующего кадра. Это делает dmamem2move(). Вот ее код:



#define DMA_STREAM               DMA2_Stream7
#define DMA_CHANNEL              DMA_Channel_1
#define DMA_STREAM_CLOCK         RCC_AHB1Periph_DMA2 

void
releasedmamemmove(void)
{
    while (DMA_GetCmdStatus(DMA_STREAM) != DISABLE) {     }
}

void
dmamem2move(int16_t *dst, int16_t *src, int size)
{
	DMA_InitTypeDef  DMA_InitStructure;
	
	DMA_DeInit(DMA_STREAM);
	
	while (DMA_GetCmdStatus(DMA_STREAM) != DISABLE)  {  }
	
	  /* Configure DMA Stream */
  DMA_InitStructure.DMA_Channel = DMA_CHANNEL;  
  DMA_InitStructure.DMA_PeripheralBaseAddr = (uint32_t)src;
  DMA_InitStructure.DMA_Memory0BaseAddr = (uint32_t)dst;
  DMA_InitStructure.DMA_DIR = DMA_DIR_MemoryToMemory;
  DMA_InitStructure.DMA_BufferSize = (uint32_t)size;
  DMA_InitStructure.DMA_PeripheralInc = DMA_PeripheralInc_Enable;
  DMA_InitStructure.DMA_MemoryInc = DMA_MemoryInc_Enable;
  DMA_InitStructure.DMA_PeripheralDataSize = DMA_PeripheralDataSize_HalfWord;
  DMA_InitStructure.DMA_MemoryDataSize = DMA_PeripheralDataSize_HalfWord;
  DMA_InitStructure.DMA_Mode = DMA_Mode_Normal;
  DMA_InitStructure.DMA_Priority = DMA_Priority_High;
  DMA_InitStructure.DMA_FIFOMode = DMA_FIFOMode_Disable;         
  DMA_InitStructure.DMA_FIFOThreshold = DMA_FIFOThreshold_Full;
  DMA_InitStructure.DMA_MemoryBurst = DMA_MemoryBurst_Single;
  DMA_InitStructure.DMA_PeripheralBurst = DMA_PeripheralBurst_Single;
  DMA_Init(DMA_STREAM, &DMA_InitStructure);
    
   /* DMA Stream enable */
  DMA_Cmd(DMA_STREAM, ENABLE);
}


Смотрим, сколько времени занимает эта процедура… 1мс! Много… коэффициенты в RAM — 820мкс. Уже лучше.

Выходной фильтр.

После обработки надо передескретизировать сигнал с 8кГц до 128кГц. Это делается так — исходный сигнал дополняется нулями так, чтобы на один отчсет выходного сигнала приходилось 1 входной + 15 нулей. 1+15 == 16 == 128/8.
Потом фильтруется выходной сигнал фильтром, имеющим частоту среза примерно равную макс. частоте входного сигнала.

При этом реально нам не надо добивать вход нулями — мы будем использовать просто каждый 16й коэффициент фильтра. Иногда это называют полифазным фильтром.

Будем использовать тот же фильтр. Вот, собственно, код:



#define POSTINBUF_SIZE				(8*16)	// 128
#define SATU16(x)				(((x) > 65535)?(65535):(((x)<0)?(0):(x)))

extern int16_t pf[PRE_FILTER_LEN];

int16_t *
post_filter(int16_t *data, int16_t *outbuf, int16_t *inbuf)
{
	int i, j, k;
	int16_t *b = inbuf + POSTINBUF_SIZE;
	int32_t	sum;
	
	dmamem2move(b, data, POSTINBUF_SIZE);
	releasedmamemmove();

	for (i = 0; i < FRAME_SIZE; i += 16)
	{
		for (k = 0; k < 16; k++)
		{
			sum = 0;
			for (j = 0; j < 8; j++)
			{
				sum += (int32_t)pf[(j<<4) + k] * b[(i>>4)-j];			// 
			}

			sum += 1 << (15-1);	// Round
			sum >>= 15+0;				// shift to q1.15
			outbuf[i+k] = SATU16(sum + 32767);		// store
		}
	}
	dmamem2move(inbuf, b, POSTINBUF_SIZE);
	return outbuf;
}


Строка

outbuf[i+k] = SATU16(sum + 32767);

добавляет к выходному отсчету постоянную составляющую, ибо наш ЦАП на вход принимает числа от 0 до 4095 прижатые влево к границе 16и бит. Ну и, на всякий случай, насыщает…

Чтобы посмотреть сколько «мусора» попадет к нам в уши, подадим на вход фильтра белый шум @8кГц и посмотрим спектр выходного сигнала:



Т.е. все то, что выше 4кГц есть зеркально отраженный спектр нашего исходного сигнала т.е. алиас. Из картинки видно. что его не много.

По времени эта процедура занимает 1.1мС. Приемлемо.

Осталось это правильно вызвать:

frame process.



struct bufs_s {
	int16_t in1[FRAME_SIZE + PRE_FILTER_LEN];
	int16_t in2[FRAME_SIZE + PRE_FILTER_LEN];
	int16_t out1[OUTBUF_SIZE];
	int16_t out2[OUTBUF_SIZE];
	int16_t postbuf[POSTINBUF_SIZE * 2];
};

struct bufs_s buf;

void
frame_process(void *mem, int16_t *in_data, int16_t *out_data)
{
	int i;
	int16_t *fr_data;
	struct supl_s *supl = get_supinfo();
	
	deb_port_on();		
////////////////////////////////////	

	fr_data = pre_filter(in_data, buf.in1, buf.out1);	// returns buf.out1
	
	nr_set_params(mem,supl);

	frame_process_nr(mem,fr_data,fr_data);

	sql_gate_control(get_sql(mem));

	post_filter(fr_data, out_data, buf.postbuf);

////////////////////////////////
	deb_port_off();	
}



Собственно вот… Если исключить из рассмотрения 3 неописанные функции, то у нас получится копирование входа на выход.

Функция nr_set_params(mem,supl); устанавливает текущий уровень шумоподавления и уровень squelsh.
frame_process_nr(mem,fr_data,fr_data); собственно давит шум.
sql_gate_control(get_sql(mem)); дергает ногу процессора — есть или нет голоса.

— Сам шумодав в следующей части… если интересно… Замечу, что пре и пост обработка заняла 2мс. Т.е. на шумодав нам осталось 14мс или же 14/16*128 == 112МГц.
  • +4
  • 11 января 2013, 20:41
  • diwil

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

RSS свернуть / развернуть
Статья нормальная. Только теги поправь.
0
По-прежнему очень интересно!

Впрочем, по перечислению того, что валяется на столе в одном из комментариев к прошлым постам, можно догадаться, что автор не просто эмбеддер, а DSP'шник.

Спасибо!
0
А вообще любопытно было бы посмотреть реализацию этого же на TMS320F28xx и сравнение его с STM32F4 на этой задаче.
0
Так, а где заявленный в предыдущих статьях спектральный метод? Пока на сколько я понял просто тупо по фильтровали и по уровню вкл, выкл шумоподавитель… или я чето не так понял?
0
Пока мясо нам не показали. Пока нам показали всё обёртку — подготовку сигнала к обработке и вывод обратно. Мясо, я так понимаю, будет в следующей статье. Пока — каркас-скелет.
0
Запускаем матлаб, fdatool и придумываем…
А есть какие нибудь простые утилиты для расчета коэффициентов для цифровых фильтров? А то не у всех есть матлаб.
0
WinFilter
0
А не целесообразнее ли разбить фильтр по каскадам?? 128-й порядок многовато как-то. Я сейчас пишу программу для SSB трансивера на STM32f407, ВЧ частота 96 кГц, звуковая частота 8 кГц. Передискретизация в 12 раз получается. Так я решил остановится на трех каскадах интерполяции/децимации (3/2/2). Притом применил еще разбитие каждого фильтра еще на два каскада, в матлабе в fdatool это реализуется с помощью Interpolate FIR. В итоге при пересчете это становится эквивалентно фильтру ~40 порядка, при том что полоса фильтра 1700 Гц, а подавление 90 Дб. А вот по поводу шумоподавителя решил не заморачиваться и сделать пороговым.
0
можно и так, но:
1. 128 — это фигня
2. по каскадам — отжор памяти будет + возрастет задержка.

а для трансивера шумодав вообще — зло!
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.