Спроектируем фильтр N+1-го порядка класса FIR (Finite Impulse Response или с конечной импульсной характеристикой). Документ, решающий эту задачу, представлен в файле fir.
Каждая из N ячеек временной задержки фильтра удовлетворяет следующей зависимости выходного сигнала у от входного х вида:
Подключим пакет расширения plots, нужный для графической визуализации проектирования:
> restart:with(plots):
Warning, the name changecoords has been redefined
Зададим исходные данные для проектирования полосового цифрового фильтра, выделяющего пятую гармонику из входного сигнала в виде зашумленного меандра с частотой 500 Гц:
> N := 64: # Число секций фильтра (на 1 меньше порядка фильтра)
> fs := 10000: # Частота квантования
> fl := 2300: # Нижняя граничная частота
> fh := 2700: # Верхняя граничная частота
> m := 10: # 2^m > N - число точек для анализа
Вычислим:
> Т := 2^m-1;
T := 1023
> F1 : = evalf(fl/fs);
F1 := .2300000000
> F2 := evalf(fh/fs);
F2 = .2700000000
> Dirac(0) := 1: # Функция Дирака
> fp1:=2*Pi*F1: fp2:=2*Pi*F2:
Зададим характеристику полосового фильтра:
> g : = (sin(t*fp2)-sin(t*fp1))/(t*Pi);
Вычислим FIR коэффициенты для прямоугольного окна фильтра
> С := (n) -> limit(g,t=n): h := aray(0..N): N2:=N/2:
> for n from 0 to N2 do h[N2-n]:= evalf(C(n)); n[N2+n] :=
h[M2-n]; od:
Определим массивы входного x(n) и выходного y(n) сигналов:
> х := array(-N..Т): y:= array(0..T):
Установим значение x(n) равным 0 для времени меньше 0 и 1 для времени t>=0.
> for n from -M to -1 do x[n] := 0; od:
> for n from 0 to T do x[n] := Dirac(n); od:
Вычислим временную зависимость для выходного сигнала.
> for n from 0 to Т do y[n] := sum(h[k]*x[n-k], k=0..N); od:
Построим график импульсной характеристики фильтра, отражающей его реакцию на сигнал единичной площади с бесконечно малым временем действия:
> р := [seq{[j/fs, y[j]], j=0..Т)]:
> plot(р, time=0..3*N/fs, labels=[time,output], axes=boxed,
xtickmarks=4, title=`Импульсная характеристика фильтра`, color=black);
Он показан на рис. 11.33. Нетрудно заметить, что эта характеристика свидетельствует об узкополосности фильтра, поскольку его частоты fl и fh различаются не сильно. В этом случае полосовой фильтр по своим свойствам приближается к резонансному, хотя само по себе явление резонанса не используется.
Рис. 11.33. Импульсная характеристика цифрового фильтра
Вычислим АЧХ фильтра, используя прямое преобразование Фурье. Оно, после подготовки обрабатываемых массивов, реализуется функцией FFT:
> ro := array(1..Т+1):io := array(1..Т+1):
> for n from 0 to T do ro[n+1] := y[n]; io[n+1] := 0; od:
> FFT(m,ro,io):
Построим график АЧХ фильтра:
> р := [seq( [j*fs/(Т+1),abs(ro[j + 1]+io[j + 1]*I)], j=0..T/2)]:
> plot(p, frequency=0..fs/2, labels=[frequency,gain], title=`AЧX
фильтра`, color=black);
Он представлен на рис. 11.34. Нетрудно заметить, что и впрямь АЧХ фильтра напоминает АЧХ резонансной цепи — она имеет вид узкого пика. Вы можете легко проверить, что раздвижением частот fl и fh можно получить АЧХ с довольно плоской вершиной и резкими спадами (говорят, что такая характеристика приближается к прямоугольной).
Рис. 11.34. АЧХ цифрового полосового фильтра
Теперь приступим к тестированию фильтра. Зададим входной сигнал в виде зашумленного меандра с частотой 500 Гц и размахом напряжения 2 В:
> l := round(fs/2/500):
> for n from 0 by 2*1 to T do
> for n2 from 0 to l-1 do
> if n+n2 <= T then
> x[n+n2] := evalf(-l+rand()/10^12-0.5);
> fi;
> if n+n2+1 <= T then
> x[n+n2+1] := evalf(1+rand()/10^12-0.5);
> fi;
> od:
> od:
Временная зависимость синтезированного входного сигнала представлена на рис. 11.35.
Рис. 11.35. Синтезированный входной сигнал
Вычислим реакцию фильтра на входной сигнал:
> for n from 0 to Т do
> y[n] : = sum(h[k]*х[n-k], k=0..N);
> od:
Построим график выходного сигнала.
> р := [seq([j/fs, x[j]], j=0..T)]: q:=[seq([j/fs, y[j]], j=0..T)]:
> plot(p, time=0..T/fs/4, labels=[time,volts], title= `Входной сигнал`, color=black);
> plot(q, time=0..T/fs/4, labels=[time, volts], title=`Выгодной сигнал`, color=black);
Временна́я зависимость выходного сигнала показана на рис. 11.36. Нетрудно заметить, что, в конце концов, выходной сигнал вырождается в пятую гармонику входного сигнала, но этому предшествует довольно заметный переходной процесс. Он связан с узкополосностью данного фильтра.
Рис. 11.36. Временна́я зависимость выходного сигнала цифрового фильтра
Вычислим спектры входного и выходного сигналов, подготовив массивы выборок сигналов и применив прямое преобразование Фурье с помощью функции FFT:
> ri := array(1..Т+1): ii := array(1..Т+1):
> for n from 0 to T do
> ri[n+1] := x[n]*2/T; ii[n+1] := 0;
> ro[n+1] := y[n]*2/T; io[n+1] := 0;
> od:
> FFT(m, ri, ii): FFT(m,ro,io):
Построим график спектра входного сигнала, ограничив масштаб по амплитуде значением 0.5 В:
> р := [seq([j*fs/(T+1), abs(ri[j+1]+ii[j+1]*I)], j=0..T/2)]:
> q := [seq([j*fs/(T+1), abs(ro[j+1]+io[j+1]*I)], j=0..T/2)] :
> plot(p, frequency=0..fs/2, y=0..0.5, labels=[частота, V],
titles=`Частотный спектр входного сигнала`, color=black);
Этот график представлен на рис. 11.37. Из него хорошо видно, что спектральный состав входного сигнала представлен только нечетными гармониками, амплитуда которых убывает по мере роста номера гармоники. Пятая гармоника на частоте 2500 Гц находится посередине полосы пропускания фильтра, ограниченной граничными частотами фильтра 2300 и 2700 Гц. Заметны также беспорядочные спектральные линии шума сигнала в пределах полосы прозрачности фильтра.
Рис. 11.37. Спектрограмма входного сигнала
Теперь построим график спектра выходного сигнала:
> plot(q, frequency=0..fs/2, y=0..0.5, labels=[частота, V],
title=`Частотный спектр выходного сигнала`, color=black);
Он представлен на рис. 11.38. Хорошо видно эффективное выделение пятой гармоники сигнала и прилегающей к ней узкой полосы шумового спектра.
Рис. 11.38. Спектрограмма выходного сигнала цифрового фильтра
Приведенные данные свидетельствуют, что спроектированный фильтр полностью отвечает заданным требованиям и обеспечивает уверенное выделение пятой гармоники зашумленного меандра. По образу и подобию данного документа можно выполнить проектирование и других видов цифровых фильтров.
11.3.5. Моделирование цепи на туннельном диоде
А теперь займемся моделированием явно нелинейной цепи. Выполним его для цепи, которая состоит из последовательно включенных источника напряжения Es, резистора Rs, индуктивности L и туннельного диода, имеющий N-образную вольт-амперную характеристику (ВАХ) — см. файл tdc. Туннельный диод обладает емкостью С, что имитируется конденсатором С, подключенным параллельно туннельному диоду.
Пусть ВАХ реального туннельного диода задана выражением:
> restart;
> А:=-3: а:=10: В:=1*10^(-8): b:=20:
> Id:=Ud->A*Ud*exp(-a*Ud)+B*(exp(b*Ud-1));
Id := Ud →A Udе(-aUd) + Вe(bUd-1)
Построим график ВАХ:
> plot(Id(Ud), Ud=-.02..0.76, color=black);
Этот график представлен на рис. 11.39. Нетрудно заметить, что ВАХ туннельного диода не только резко нелинейна, но и содержит протяженный участок отрицательной дифференциальной проводимости, на котором ток падает с ростом напряжения на диоде. Это является признаком того, что такая цепь способна на переменном токе отдавать энергию во внешнюю цепь и приводить к возникновению колебаний в ней различного типа.
Рис. 11.39. ВАХ туннельного диода
Работа цепи описывается системой из двух дифференциальных уравнений:
di/dt= (Es-i(t) * Rs-u(t))/L
du/dt= (i(t)-Id(u(t))/C
Пусть задано Es=0,35 В, Rs=15 Ом, C=10∙10-12, L=30∙10-9 и максимальное время моделирования tm=10∙10-9. Итак, задаем исходные данные:
> Es: = .35: Rs:=15: C:=10*10^(-12): L:=30*10^(-6): tm:=10*10^(-9):
Составим систему дифференциальных уравнений цепи и выполним ее решение с помощью функции dsolve:
> se:=diff(i(t),t) = (Es-i(t)*Rs-u(t)) /L,
difff(u(t), t) = (i(t)- Id(u(t))) / С;
> F:=dsolve({se, i(0)=0, u(0)=0}, {i(t),u(t)}, type=numeric,
method-=classiccal, stepsize=10^(-11), output=listprocedure);
F := [t = (
proc(t) …
end proc), u(t) = () i(t) = (
proc(t) …
end proc)]
Поскольку заведомо известно, что схема имеет малые значения L и С мы задали с помощью параметров достаточно малый шаг решения для функции dsolve — stepsize=10^(-11) (с). При больших шагах возможна численная неустойчивость решения, искажающая форму колебаний, получаемую при моделировании. Используя функции odeplot и display пакета plots построим графики решения в виде временных зависимостей u(t) и 10∙i(t) и линии, соответствующей спряжению Es источника питания:
