Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:
AE_ME / Фильтры.rtf
Скачиваний:
16
Добавлен:
11.04.2015
Размер:
1.29 Mб
Скачать

КОНДРАТЬЕВ В.П.

Моделирование цифровых фильтров в maple

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

Усредняющие фильтры нижних частот. Пусть исходный сигнал задан дискретными отсчетами x0, x1, …, xN. Под фильтром обычно понимается диcкретная система, генерирующая выходную последовательность y0, y1, …, yN. Функция системы заключается в избирательном выделении либо подавлении некоторого диапазона частот в исходной последовательности. Наиболее просто фильтр нижних частот (сглаживающий фильтр) строится по методу усреднения отсчетов:

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

(1)

которая называется дискретной сверткой (см., например, [1, глава 4]). Здесь весовые коэффициенты hn, представляют собой отсчеты импульсной характеристики фильтра. Заметим, что в практических задачах в случае, когда известна вся входная последовательность целиком, применяется сглаживающий фильтр, использующий одинаковое количество отсчетов слева и справа от текущего отсчета (метод скользящего среднего).

Рассмотрим наиболее простой пример фильтра нижних частот с усреднением по трем отсчетам. Для усреднения по трем точкам зададим порядок фильтра M=2 и определим значения отсчетов импульсной характеристики которые являются весовыми коэффициентами в (1). Одним из способов оценки качества фильтра является изучение его комплексной частотной характеристики (ЧХ), которая вычисляется (см. [2, лекция 19]) как прямое дискретное преобразование Фурье (ДПФ) последовательности коэффициентов фильтра (то есть его импульсной характеристики):

(2)

где N=M+1 — число отсчетов импульсной характеристики, j — мнимая единица, величина соответствует шагу дискретизации T по времени, а количество коэффициентов Hm =H(m) ДПФ достаточно вычислять для Обычно переходят от дискретной последовательности номеров m=0, 1, … в частотной области к круговой частоте, возвращаясь к формулам преобразования Фурье (не дискретного), получая тем самым комплексную частотную характеристику непрерывного изменения частоты (см. [2, лекция 6]):

(3)

либо к нормированной частоте =T, где [0, ] :

(4)

а суммирование можно выполнять только по ненулевым коэффициентам импульсной характеристики, поскольку при k>N-1 имеем h(k)=0.

Модуль и аргумент частотной характеристики представляют собой

(аналогично спектрам амплитуд и фаз как функций дискретного номера гармоники), соответственно, амплитудно-частотную (АЧХ)

(5)

и фазо - частотную (ФЧХ)

(6)

характеристики фильтра, как функции непрерывного аргумента — частоты .

АЧХ импульсной характеристики можно получить из явных формул (5 - 6), либо применяя алгоритм ДПФ. Получим частотную характеристику фильтра, вычислив дискретное преобразование Фурье последовательности коэффициентов импульсной характеристики hk , дополненной нулями при k > N-1. Вычислив ДПФ этой последовательности, построим график ее АЧХ, которая, как известно (см., например, [2, лекция 11], [1, глава 5]), представляет собой дискретную периодическую (период длины N) и симметричную относительно середины последовательность коэффициентов Ck ДПФ (с точностью до постоянного нормирующего множителя).

Отметим также важную особенность: обе последовательности связаны как взаимно обратные относительно прямого ДПФ и обратного ДПФ, примененного к каждой из них, соответственно.

Приведем текст рабочего листа сеанса Maple, команды которого используют процедуры вычисления ДПФ и строят графики АЧХ и ФЧХ для различных случаев фильтров. Строки с командами рабочего листа будем снабжать номерами, на которые делаются ссылки в поясняющем тексте. Двоеточие в конце команды подавляет вывод результата ее выполнения. Большинство команд, вычисляющих вспомогательные величины, а также команды построения графиков, заканчиваются двоеточием. При самостоятельном выполнении команд рабочего листа двоеточие следует заменить на точку с запятой.

Строка 1 содержит процедуру построения ДПФ дискретной последовательности. В строках 4 – 6 описаны процедуры построения частотных характеристик обрабатываемой последовательности. Процедура ARR (строка 2) строит график АЧХ последовательности.

> restart;with(plots):with(plottools):

> 1: DTF:=proc (y,N,Y) local n,k,j,p,h; n:=N-1; h:=2*Pi/N;

1.1: for k from 0 to n do

p:=0; for j from 0 to n do

p:=p+evalf(y[j]*exp(-I*k*j*h));

end; Y[k]:=evalf(1/N*p);

end:end:;

> 2: ARR:=proc(n::integer,c) local L,H,ma,mi,k::integer,Sim::array(0..n);

ma:=c[0];mi:=c[0]; L:=line([0,c[0]],[n,c[n]],thickness=2,color=red);

for k from 1 to n do

if c[k]>ma then ma:=c[k];end if; if c[k]<mi then mi:=c[k];end if;

end do;

H:=ma-mi; if H=0 then RETURN(L) end;

for k from 0 to n do

if abs(c[k])<H/1000 then

Sim[k]:=ellipse([k,c[k]],0.2,0.01*H,color=blue);

else Sim[k]:=arrow([k,0],[k,c[k]],0.2,0.2,0,color=black);

end if;

end do; convert(Sim,list);

end:

> 3:Spectr:=proc(n,C,A,phi) local k,R;global Risphi;

6.1: for k from 0 to n do A[k]:=evalf(abs(C[k])):

phi[k]:=evalf(argument(C[k])); R[k]:=[eval(k),eval(phi[k])];

end:;

Risphi:=plot(convert(R,list),thickness=2,color=blue,style=point,

symbol=box):

end:

> 4: H:=proc(omega,h,N) local n,z;

z:=sum(h[n]*exp(-I*omega*n),n=0..N-1);

end:

> 5: Amp:=proc(omega,h,N) local z; global H;

z:=evalf(abs(H(omega,h,N)));

end:

> 6: phaza:=proc(omega,h,N) local z; global H;

z:=evalf(argument(H(omega,h,N)));

end:

> 7: frect:=proc(t) piecewise(t<2*Pi*Fc,1,0) end proc:

Зададим (строки 8 – 9) частоту среза (называемую также частотой Найквиста или граничной частотой полосы пропускания), а также частоту дискретизации, которая (в согласии с теоремой Котельникова) должна быть по крайней мере в два раза выше частоты среза. Затем зададим шаг дискретизации, называемый иначе периодом дискретизации. Поскольку нас интересует математическая сторона процесса цифровой фильтрации, физические размерности этих величин (секунды — для временной оси, герцы — в частотной области) мы не будем принимать во внимание. Определим порядок фильтра M=2 (строка 12) и зададим коэффициенты импульсной характеристики для усреднения по трем отсчетам (строки 14 – 16), дополнив последовательность hk нулями. Вычислим ДПФ этой последовательности и построим (рисунок 1) ее частотные характеристики (строки 20 – 22).

> 8: Fc:=8:Fd:=4*Fc:#частота среза и частота дискретизации

> 9: T:=1/Fd;omega[d]:=2*Pi*Fd:# шаг и круговая частота

10:Risrect:=plot(frect(f*2*Pi),f=0..Fd/2,color=brown, thickness=2):

> 11:Risimp:=proc(N,himp) local Gra,GrafF,Graflin;

Gra:=zip((x,y)->[x,y],[seq(k,k=0..N)],convert(himp,list)):

GrafF:=plot(Gra,style=point,symbol=box,symbolsize=14,

color=blue):

Graflin:=plot(Gra,color=blue): display(GrafF,Graflin):

end:

> 12: M:=2:

> 13: h:=array(0..M):

> 14: h[0]:=1/3:h[1]:=1/3:h[2]:=1/3:

> 15: N:=23:himp:=array(0..N,sparse):

> 16: for n from 0 to M do himp[n]:=h[n];

end do:

> 17: Risimp(N,himp):

> 18:n:=N-1:

> 19: C:=array(0..n):A:=array(0..n):;

> 20: DTF(himp,N,C):

> 21: Spectr(n,C,A,phi):

> 22:display(ARR(n,A));

Рисунок 1. АЧХ и ФЧХ трехточечного фильтра.

Построение графика АЧХ (рисунок 2) как функции непрерывной частоты (строки 24 – 27) показывает, что она далека от "идеальной" АЧХ фильтра нижних частот (строка 7), имеющей прямоугольную форму. ФЧХ (рисунок 1) имеет нечетную симметрию относительно середины числа отсчетов и испытывает скачки на величину в полосах задерживания и переходных полосах ([2, лекция 18]). На рисунках в целях наглядности по оси абсцисс графика откладываются абсолютные значения частот, что позволяет оценить величину полосы пропускания сравнительно с заданной частотой среза Fc=8. Обычно используют шкалы нормированных частот f [0, 0.5] либо [0, ], позволяющие абстрагироваться от значений частоты и оценивать относительные размеры величин области пропускания, переходной области и полосы подавления.

Увеличим порядок фильтра и построим график АЧХ, сравнивая его с предыдущим графиком. Строки команд, выполняющие аналогичные предыдущим вычисления, снабдим номерами, отличающимися на 100.

> 112: M:=6:

> 113: h:=array(0..M):

> 114: for k from 0 to M do h[k]:=1/(M+1):end do:;

> 115: N:=35:himp:=array(0..N,sparse):

> 116: for j from 0 to M do himp[j]:=h[j]; end do:

> 117:Risimp(N,himp):

> 126: Ris2:= plot(Amp(2*Pi*f/Fd,himp,N),f=0..Fd/2,color=black,

thickness=3,linestyle=4):;

> 127:display(Risrect,Ris1,Ris2);

Сравнение графиков АЧХ, выведенных командой строки 127, показывает, что фильтр большего порядка (пунктирная линия на рисунке 3) имеет меньшую полосу пропускания, а величина боковых лепестков на верхних частотах существенно не уменьшается. Тем не менее, простое перераспределение величин коэффициентов импульсной характеристики (рисунок 4, выведенный командой строки 217) уже улучшает качество фильтра (строки 214 -- 227).

> 214: h[0]:=1/20:h[1]:=3/20:h[2]:=1/5:h[3]:=1/5:h[4]:=h[2]:h[5]:=h[1]:

h[6]:=h[0]:

> 215: N:=35:himp:=array(0..N,sparse):s:=0:

> 216: for n from 0 to M do

himp[n]:=h[n];s:=s+h[n];

end do: s;

> 217:Risimp(N,himp):

> 226: Ris3:= plot(Amp(2*Pi*f/Fd,himp,N),f=0..Fd/2,color=black,thickness=2):

> 227:display(Risrect,Ris1,Ris2,Ris3);

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

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

Рисунок 2.

Рисунок 3.

Рисунок 4.

Оконные функции. Для получения АЧХ формы, близкой к прямоугольной, возьмем в качестве значений коэффициентов импульсной характеристики дискретные отсчеты функции , являющейся Фурье-преобразованием прямоугольного сигнала. Повторим вычисления строк 8 – 27 для этого случая, приводя только измененные команды соответствующих строк.

> 8: Fc:=50:Fd:=4*Fc:#частота среза и частота дискретизации

> 9: T:=1/Fd:;omega[d]:=2*Pi*Fd:# шаг и круговая частота

> 12: M:=18:

> 13: h:=array(0..M):c:=array(0..M):;

> omega_n:=2*Pi*Fc:

> c[0]:=omega_n*T/Pi;h[M/2]:=c[0]:K:=M/2:

> 14: for n from 1 to M/2 do

c[n]:=evalf(sin(omega_n*T*n)/(Pi*n));

end:;

> s:=c[0]:for n from 0 to M/2-1 do

h[n]:=c[K-n];h[M-n]:=h[n]; s:=s+2*h[n];

end:;print(s);

> 15: N:=M:

> 17: Risimp(N,h);

> 18:n:=N-1:

> 19: C:=array(0..n):A:=array(0..n):;

> 20: DTF(h,N,C):

> 21: Spectr(n,C,A,phi):

> 22:display(ARR(N/2,A)):;

> 23: plot(Amp(omega,h,N),omega=0..Pi):

> 24: plot(Amp(2*Pi*f/Fd,h,N),f=0..Fd/2):

> 27: display(Ris1,Risrect);

Команда строки 17 выводит график импульсной характеристики (рисунок 5), а в строке 27 командой display строятся графики АЧХ фильтра в сравнении с прямоугольной АЧХ (рисунок 6).

Рисунок 5.

Рисунок 6.

Низкий уровень боковых лепестков АЧХ позволяет подавить мощность верхних частот, но вблизи точки разрыва прямоугольной АЧХ наблюдается эффект Гиббса, связанный с усечением ряда Фурье до конечного числа слагаемых. Увеличение порядка фильтра не устраняет этот эффект (рисунки 7 – 8, построенные для новых параметров командами строк 12 – 27, которые полностью аналогичны приведенным выше командам). Изменив порядок фильтра (строка 12), выполним вычисления строк 12 – 27 и построим аналогичные графики импульсной характеристики (рисунок 7) и АЧХ фильтра (рисунок 8) при новых значениях параметров задачи.

> 12: M:=60:

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

> 400: ThetaH:=proc(n) local z,alpha;alpha:=0.54;

z:=piecewise(abs(n)>M/2,0,alpha+(1-alpha)*cos(2*Pi*n/M));

evalf(z);end:;

> 500:ThetaB:=proc(n) local z;

z:=piecewise(abs(n)>M/2,0,0.42+0.5*cos(2*Pi*n/M)+

0.08*cos(4*Pi*n/M));

evalf(z);end:

> 600:ThetaK:=proc(n) local k,z,beta,J0;beta:=1.2;

J0:=proc(x) local y;y:=1+sum((1/k!*(x/2)^2)^k,k=1..7);evalf(y);end;

z:=piecewise(abs(n)>M/2,0,J0(beta*sqrt(1-(2*n/M)^2))/J0(beta));

evalf(z);end:;

Рисунок 7.

Рисунок 8.

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

Рисунок 9.

Рисунок 10.

Строка 14 имеет вид

> 14: for n from 1 to M/2 do

c[n]:=evalf(sin(omega_n*T*n)/(Pi*n))*ThetaH(n);

end:;

Соседние файлы в папке AE_ME