Разработка цифрового полосового КИХ-фильтра

       H(e^jω) – действительная частотная характеристика,

       W(e^jω) – положительная  весовая функция, которая позволяет выбирать разный вес ошибки для  различных участков полосы частот.

Значение  р =2  соответствует квадратичной   норме ошибки.  При p → ∞ соответствующая норма дает максимальное  по модулю абсолютное  отклонение  желаемой  частотной характеристики от действительной, т.е.   . 

Цель оптимального чебышевского метода синтеза – найти коэффициенты фильтра,  т.е. значения импульсной характеристики h[n], при которых значение  максимальной ошибки  E(ω) будет минимальным,       т.е.

Рис. 2. Графическое объяснение метода

 

2. Проектирование фильтра

 

Спецификация  цифрового фильтра, в соответствии с техническим заданием:

 

Fp1=400;         %Гц, нижняя граничная частота полосы пропускания

Fp2=800;         %Гц, верхняя граничная частота полосы пропускания

Fs1=300;         % Гц, нижняя граничная частота полосы задерживания

Fs2=900;        % Гц, верхняя граничная частота полосы задерживания

Ft=8000;       % Гц, частота дискретизации

Ap=1;              % Дб, неравномерность передачи в полосе пропускания

As=60;             % Дб, минимальное ослабление в полосе задерживания

Tt=1/Ft;           % сек, период дискретизации

Fc1=(Fp1+Fs1)/2=350;   %Гц, нижняя частота среза

Fc2=( Fp2+Fs2)/2=850;   %Гц, верхняя частота среза

 

С помощью  следующего script-файла проведём сравнение различных методов разработки КИХ-фильтров:

 

% Программа сравнения  различных методов проектирования

Ft=8000;

% фильтр с окном Кайзера

fcuts=[300  400 800 900];

mags=[0  1  0];

defs=[0.001 0.01  0.001];

[n, Wn, beta,ftype]=kaiserord(fcuts,mags,defs,Ft);

hk=fir1(n, Wn, ftype, kaiser(n+1,beta), 'noscale');

fvtool(hk,1);

legend('Фильтр с окном Кайзера')

% Оптимальный фильтр по  алгоритму Паркса- МакКлелана

As=60; Ap=1;

f=[300  400  800  900];

a=[0  1  0];

dev=[10^(-As/20),  1-10^(-Ap/20),   10^(-As/20)];

[n, fn, a0, w]= firpmord(f, a, dev, Ft);

[b error opt]=firpm(n,fn,a0,w);

fvtool(b,1) 

legend('Оптимальный фильтр')

 

Результат выполнения:

 

Рис. 3. АЧХ полученного фильтра

а) Фильтр с  окном Кайзера

б) Оптимальный  фильтр

 

 

При этом порядок фильтра, полученный при  помощи функции kaiserord(), составляет 291. При помощи функции firpmord() – всего лишь 141.

 

Рассмотрим  ближе АЧХ фильтра с окном  Кайзера.

Рис. 4. Увеличенный  фрагмент ПП АЧХ фильтра с окном  Кайзера

 

Рис. 5. Увеличенный  фрагмент ПЗ АЧХ фильтра с окном  Кайзера

 

Из рис. 5 видно, что минимальное ослабление в полосе задерживания не соответствует  требуемой спецификации.

 

 

Рис. 6. Увеличенный фрагмент полосы пропускания АЧХ оптимального фильтра

 

Рис. 7. Увеличенный фрагмент полосы задерживания АЧХ оптимального фильтра

 

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

Так как  оба спроектированных фильтра нам  не подходят, то попробуем увеличить порядок до 151 по сравнению  с вычисленным функцией     firpmord() и заново пересчитаем фильтр:

 

Рис. 8. АЧХ оптимального фильтра с порядком 150

 

Рис. 9. Увеличенные фрагменты АЧХ оптимального фильтра порядка 150

 

Из увеличенных  фрагментов (рис. 9) видно, что фильтр удовлетворяет условиям неравномерности в полосе пропускания и минимальному ослаблению в полосе задерживания. Увеличим фрагменты АЧХ в линейном масштабе, чтобы наглядно проверить соответствие ПП и ПЗ.

 

Рис. 10. АЧХ  спроектированного фильтра в  линейном масштабе

 

Рис. 11. Увеличенные  фрагменты граничных частот ПЗ и  ПП

 

По представленным выше графикам, можно сделать вывод  о соответствии спроектированного  фильтра требуемой спецификации. По графику АЧХ видно соответствие граничных частот: Fp1=300 Гц, Fp2=800 Гц, Fs1=400 Гц, Fs2=900 Гц (рис.11). Графики ЛАЧХ  показывают, что неравномерность передачи в полосе пропускания составляет не более 1 дБ, а ослабление в полосе задерживания составляет не менее 60 дБ, что говорит о правильности расчета данного фильтра.

Фильтр с окном Кайзера не  полностью удовлетворяет спецификации по затуханию в полосе задерживания.  Оптимальный фильтр с порядком, полученным при помощи функции firpmord(),  не подходит по условиям ПП и ПЗ. Наилучшую результирующую  АЧХ имеет чебышевский  фильтр с порядком 150, он полностью соответствует требованиям.

 

3. Уравнение фильтра

 

Окончанием  проектирования фильтра можно считать  получение его уравнения. Получим  же уравнение нашего спроектированного  фильтра:

 

       Уравнение полосового КИХ-фильтра имеет вид:

,     N -1 – порядок, h[n] - коэффициенты фильтра. Уравнение фильтра –  дискретная свертка.

 

 

Рис. 12. Структурная  схема   КИХ -  фильтра для порядка 150

 

 

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

 

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

 

4. Тестирование фильтра

 

Фильтр  спроектирован. Проверим его работу.

А) Подаём сигнал с частотой, из полосы пропускания, 600 Гц

 

n=0:1:400;

t=1/8000;

subplot(2,1,1)

x=cos(2*pi*600*n*t);

plot(n,x)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал с частотой 600 Гц')

subplot(2,1,2)

F=filter(b,1,x);

plot(n,F)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал на выходе фильтра ')

axis([0,400,-1,1]);

 

Рис. 13. Фильтрация сигнала из полосы пропускания

 

По графикам, представленным, на рис. 13 видно, что  такой сигнал фильтр пропускает без  искажений.

 

Б) Подаём сигнал с частотой, из полосы задерживания, 2000 Гц

 

n=0:1:400;

t=1/8000;

subplot(2,1,1)

x=cos(2*pi*2000*n*t);

plot(n,x)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал с частотой 2000 Гц')

subplot(2,1,2)

F=filter(b,1,x);

plot(n,F)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал на выходе фильтра ')

axis([0,400,-1,1]);

 

Рис. 14. Фильтрация сигнала из полосы задерживания

 

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

 

В) С помощью следующего script-файла сгенерируем сигнал, равный сумме гармоник с частотами 600 Гц и 2000 Гц, и осуществим фильтрацию данного сигнала с помощью спроектированного фильтра:

 

n=0:1:400;

t=1/8000;

subplot(4,1,1)

x=cos(2*pi*600*n*t);

plot(n,x)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал с частотой 600 Гц')

subplot(4,1,2)

y=cos(2*pi*2000*n*t);

plot(n,y)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал с частотой 2000 Гц')

subplot(4,1,3)

z=x+y;

plot(n,z)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сумма двух сигналов')

subplot(4,1,4)

F=filter(b,1,z);

plot(n,F)

set(gca, 'FontName', 'Times New Roman Cyr', 'FontSize' ,10);

legend('Сигнал на выходе фильтра ')

axis([0,400,-1,1]);

 

Результат выполнения:

 

Рис. 15. Графики результирующие применение полученного фильтра

 

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

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

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

 

5. Программирование и анализ быстродействия фильтра

 

Реализуем программно наш цифровой КИХ-фильтр нижних частот. Программирование фильтра осуществляется в среде Borland C++ (Приложение 3).

 

namespace S {

 

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

 

const int    N  = 150;

const int    Fc1= 350;

const int    Fc2=850;

const double Fd = 8000;

const double PI = 3.14159;

 

double w[150] =

{

0,-8.548e-6,-1.953e-5,5.258e-5,1.406e-4,-1.374e-5,-3.423e-4,-2.659e-4,4.151e-4,7.893e-4,-6.093e-5,-1.275e-3,-8.656e-4,1.214e-3,2.116e-3,-1.519e-4,-2.991e-3,-1.928e-3,2.587e-3,4.339e-3,-3.014e-4,-5.764e-3,-3.624e-3,4.756e-3,-2.991e-3,-1.928e-3,2.587e-3,4.339e-3,-3.014e-4,-5.764e-3,-3.624e-3,4.756e-3,7.829e-3,-5.352e-4,-0.01,-6.282e-3,8.179e-3,0.013,-9.125e-4,-0.017,-0.011,0.014,0.023,-1.613e-3,-0.031,-0.02,0.027,0.048,-3.559e-3,-0.076,-0.058,0.103,0.31,04,0.294,0.084,-0.65,-0.072,3.419e-3,0.049,0.023,-0.023,-0.03,1.595e-3,0.024,0.012,-0.012,-0.016,9.036e-4,0.014,6.909e-3,-7.285e-3,-9637e-3,5.299-4,8.045e-3,4.021e-3,-4.098e-3,-5.493e-3,2.89e-4,4.451e-3,2.18e-3,-2.229e-3,-2.844e-3,1.498e-4,2.164e-3,1.021e-3,-9.963e-3,-1.206e-3,5.972e-5,8.013e-4,3.458e-4,-3.027e-4,-3.193e-4,1.322-5,1.391e-4,4.215e-5,-2.08e-5,-6.793e-6,0   

};

 

double * GetSignal() {

    int n;

    double * p = (double*)malloc(sizeof(double)*N);

    for ( n = 0; n < N ; ++n ) {

        p[n] = cos(2*PI*600*n*1/Fd);

    }

    return(p);

}

 

double * CBEPTKA() {

    int n, k;

    double s = 0, *p, *r;

 

    r = (double*)malloc(sizeof(double)*N);

    p = GetSignal();

    if (!r || !p) {perror(NULL);return NULL;}

 

    for ( n = 0; n < N ; ++n) {

        for (k = 0; k < n; ++k) {

            s += p[n-k]*w[k];

        }

        r[n] = s;

    }

    free(p);

 

    return(r);

}

} 

 

В данном исходном коде мы задаём длину импульсной характеристики N = 150, частоту среза Fc1 = 350 Гц и Fc2=850 Гц, частоту дискретизации Fd = 8000 Гц, то есть задаём исходные требования. Далее мы определяем все 150 значений импульсной характеристики нашего фильтра, и определяем некий входной сигнал, в данной программе за входной сигнал принят уже используемый ранее нами сигнал cos(2*pi*600*n*t), частота которого входит в полосу пропускания фильтра. Для простоты, программная реализация задана в виде дискретной свёртки во времени, а не в виде вычисления свёртки в частотной области с помощью алгоритмов быстрого преобразования Фурье (БПФ). Результат программы – сигнал на выходе фильтра, выводится в окне Windows, также сразу появляется окошко, сообщающее о времени данной реализации фильтра. Код вывода результата на экран и подсчёт быстродействия имеет следующий вид:

 

WM_PAINT:

{

hdc = BeginPaint(hWnd, &ps);

RECT rt;

GetClientRect(hWnd, &rt);

SetMapMode(hdc, MM_HIMETRIC);

 

SetViewportOrgEx(hdc, 0, rt.bottom/2, NULL);

 

DWORD dwCurTime = GetTickCount(); // получаем текущее время

double * p = S::CBEPTKA();

MoveToEx(hdc,0,0,NULL);

 

for (int i = 0; i < LiS::N; ++i) { // строим график

LineTo(hdc, i*70, 100*(int)p[i]);    }

 

// заканчиваем отсчёт времени

dwCurTime = GetTickCount() - dwCurTime; // Время - время(старое)

char buf[10];

 

sprintf(buf, "%d ms", dwCurTime);

MessageBox(NULL,buf," Время вычисления", MB_OK);

// Выводим  сообщение в диалоге

 

EndPaint(hWnd, &ps);

break; 

} 

 

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