Главная > DSP > Цифровые фильтры (часть 2). Практика.

Цифровые фильтры (часть 2). Практика.

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

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

Итак, в предыдущей статье мы выяснили, что цифровой фильтр представляет собой свертку соответствующей импульсной характеристики со входным сигналом. Очевидно, спроектировать КИХ-фильтр означает как-то найти его импульсную характеристику, которая, к слову, в этом случае называется его ядром (kernel).

Кроме того, в предыдущей статье упоминалось, что импульсная характеристика и АЧХ связаны преобразованием Фурье — с помощью него первое получается из второго; соответственно, обратным преобразованием достигается переход в противоположном направлении.

Фильтр цифровой — какое ядро нарисуем, такие характеристики и будут. Никаких ограничений, казалось бы. В свете этого совершенно естественным представляется желание сделать совершенно идеальный фильтр с характеристикой-прямоугольником:

Сказано — сделано. Рисуем эту характеристику, делаем преобразование Фурье. Это можно было бы сделать с помощью FFT, однако в данном случае это будет несколько нерационально, поскольку результат преобразования Фурье для прямоугольника известен с незапамятных времен — это функция sinc(x):

Собственно, на этом почти все. Дальше надо просто посчитать значения этой функции в точках дискретизации и получится готовое ядро фильтра. Почти — потому что один нюанс все же есть. Как упоминалось в предыдущей статье, мир устроен таким образом, что преобразование Фурье конечных функций есть бесконечные функции, и наоборот. Прямоугольник, описывающей идеальную АЧХ, конечен, потому, говоря, что его Фурье-образом является функция sinc(x), мы, разумеется, полагаем ее продолжающейся по осям координат бесконечно в обе стороны.

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


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

Вот на этом точно все.

По мотивам соответствующих глав книги, упомянутой в самом начале, я написал (на Lua, разумеется) программу-генератор ядра фильтра для заданной частоты дискретизации, частоты среза и длины ядра (от нее зависит скорость спада/возрастания характеристики). На выходе получается CSV-файл с ядром и файл с массивом, записанным по правилам синтаксиса Lua. Собственно, в программе можно подсмотреть параметры окна, конкретный вид sinc-функции и прочие детали вроде алгоритма нормализации.

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

Опять же, мне лень заливать код файлом, потому я приведу его листингом в конце статьи. Но сначала покажу картинки, относящиеся к спроектированному с помощью этого кода фильтру.

Фильтр низкой частоты с частотой среза 500 Гц, для частоты дискретизации 8 кГц, единичное усиление, длина ядра — 301 выборка. Кстати, важный момент — для обеспечения линейности фазо-частотной характеристики фильтра ядро должно быть симметричным, соответственно, состоять из нечетного количества выборок. Для этого длина ядра, указываемая в программе ниже, должна быть четной — в ходе генерации добавится еще одна выборка.

Ядро:

filter_kernelТеоретическая АЧХ, полученная обратным преобразованием Фурье ядра (MathCAD):

fr

 

Результат прогона белого шума через фильтр (синий — исходный спектр реализации шума, красный — спектр на выходе фильтра):

sp

 

По осям X отложены номера выборок; самая правая точка графика соответствует половине частоты дискретизации.

Код:

 

SAMPLE_FREQ		= 8000;
CUTOFF_FREQ		= 500;
FILTER_GAIN		= 1;
KERNEL_LEN		= 300;
HIGHPASS_FILTER	        = false;
 
kernel = {};
 
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--                             Basic functions
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-- The heart of a filter, a sinc function representing impulse response
function sinc(n,fc)
  if n~=0 then
    return (math.sin(2*math.pi*fc*n)/n);
  else
    return 2*math.pi*fc;
  end;
end;
 
-- Blackman smoothing window
function Blackman_window(n,M)
  return (0.42 - 0.5*math.cos((2*math.pi*n)/M) + 0.08*math.cos((4*math.pi*n)/M));
end;
 
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--                       Filter response computation
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-- Load kernel with a shifted sinc function
 
k=0;
while k<=KERNEL_LEN do
  kernel[k]=sinc(k - KERNEL_LEN/2,CUTOFF_FREQ / SAMPLE_FREQ);
  k=k+1;
end;
 
-- Now apply smoothing
 
k=0;
while k<=KERNEL_LEN do
  kernel[k]=kernel[k] * Blackman_window(k,KERNEL_LEN);
  k=k+1;
end;
 
-- Normalize kernel for it to have an unity sum of all elements
 
sum=0;
 
k=0;
while k<=KERNEL_LEN do
	sum=sum + kernel[k];
	k=k+1;
end;
 
k=0;
while k<=KERNEL_LEN do
	kernel[k]=kernel[k] / sum;
	k=k+1;
end;
 
-- Transform low-pass kernel into a high-pass if necessary
 
if HIGHPASS_FILTER==true then
 
  k=0;
  while k<=KERNEL_LEN do
 
	if k%2==0 then
		kernel[k]=-1.0*kernel[k];
	end;
 
    k=k+1;
  end;
 
end;
 
-- Write ready kernel into files
 
kfile = io.open("kernel.csv","w");
cfile = io.open("kernel.lua","w");
 
cfile:write("\nkernel_len = " .. KERNEL_LEN + 1 .. "\n\nkernel = {\n");
 
k=0;
while kernel[k]~=nil do
  kfile:write(k,";",kernel[k],"\n");
  cfile:write(kernel[k],",\n");
  k=k+1;
end;
 
cfile:write("};\n");
 
kfile:close();
cfile:close();
 
-- Print some info
 
print("Transition bandwidth estimation: ",(4/(KERNEL_LEN-1))*SAMPLE_FREQ," Hz");
 
if ((4/(KERNEL_LEN-1))*SAMPLE_FREQ)>CUTOFF_FREQ then
	print("WARNING: Transition bandwidth is wider than cutoff bandwidth - consider using longer kernel.");
end;

 

Видно, что в программе есть переменная-переключатель, отвечающая за генерацию фильтра высоких частот (по умолчанию генерируется ФНЧ). Как это достигается? Как вообще делаются все остальные фильтры?

Преобразование ядра для ФНЧ в ядро для ФВЧ достигается умножением каждой второй его выборки на -1, что эквивалентно умножению его на синусоиду с частотой, равной половине частоты дискретизации. Это приводит к спектральному сдвигу, выглядящему как отражение АЧХ фильтра слева направо. В результате из ФНЧ получается ФВЧ.

Полосовой фильтр получается из комбинации ФНЧ и ФВЧ с соответствующими частотами среза. Для получения полосового фильтра их импульсные характеристики надо перемножить. Если нужен режекторный фильтр — характеристики складывают. Разумеется, результирующую импульсную характеристику надо заново нормировать.

Да, совсем забыл сказать. Рассмотренный в этой статье фильтр так и называется — windowed sinc filter, что отражает вид его импульсной характеристики и применение оконной функции в процессе расчета ядра оного.

Рубрики:DSP
  1. Роман
    03/02/2015 в 02:46

    код на ruby?

    • YS
      03/02/2015 в 15:56

      Боже упаси. : ) Код написан на Lua. Чудесный язык, мне очень нравится. Я его применяю как язык для скриптов, облегчающих жизнь, несложных расчетов и просто пишу на нем для души.

      Кстати, Lua применяется во встроенных системах. Его интерпретатор, в частности, портирован на TL-MR3020 и некоторые радиомодемы.

  2. Тимофей
    20/01/2015 в 14:24

    А как книга называется? Там ссылка не работает…

    • YS
      20/01/2015 в 22:28

      Хм, странно. Я проверил — у меня работает. dspguide.com, Steven W. Smith «The Scientist and Engineer’s Guide to Digital Signal Processing»

  3. Алексей
    12/06/2014 в 10:26

    Спасибо!

    • YS
      12/06/2014 в 11:47

      Не за что.🙂 На самом деле эта статья — скорее аннотация к главе упомянутой книги, в которой рассказывается про sinc-фильтры. Очень рекомендую саму книгу, там все объясняется гораздо полнее и, что важнее всего, очень доступно.

  1. No trackbacks yet.

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s