Содержание


Моделирование методов цифровой обработки сигналов на языке Perl

Часть 2. Использование математических библиотек Perl для фильтрации сигналов

Comments

Серия контента:

Этот контент является частью # из серии # статей: Моделирование методов цифровой обработки сигналов на языке Perl

Следите за выходом новых статей этой серии.

Этот контент является частью серии:Моделирование методов цифровой обработки сигналов на языке Perl

Следите за выходом новых статей этой серии.

1. Моделирование цифровой свертки на примере сглаживания сигналов.

По определению, цифровые сигналы представляют собой последовательности x(n), у которых амплитуда и время дискретны. Последовательность x(n) называется периодической с периодом N, если x(n)=x(n+N), для всех n. Фундаментальным понятием в обработке сигналов является линейная свертка. Из свойства инвариантности к сдвигу линейных систем следует, что если последовательность h(n) является импульсной характеристикой, то реакция системы определяется выражением:

 y(n) = x(k)h(n-k), (1)

которое обычно называется сверткой. Свертка двух последовательностей с периодом N, называется периодической или круговой сверткой.

В экспериментальных исследованиях одной из важнейших процедур является процедура сглаживания сигналов, которая может иметь целью как снижение уровня шума, так и его оценку. Сглаживание последовательности x(n) можно реализовать с помощью свертки с последовательностью p(n), которая представляет собой прямоугольный импульс единичной амплитуды длительностью L. На практике, при вычислении свертки сглаживания, суммирование x(n) можно осуществлять только на интервалах L, где последовательность p(n) отлична от нуля. Программная реализация на Perl свертки сглаживания "в лоб" выглядит следующим образом:

 @y = (0) x $N; $L = 6; for ($n=0; $n<$N; $n++) { $sum = 0; for ($k=0; $k<$L; $k++) { $m = $n + $k; if($m > $N-1){ last; } $sum += $x[$m]; } $y[$n] = $sum/$L; }

Как видно из кода программы, на каждом шаге, для вычисления значения выходной последовательности y(n), необходимо выполнить L сложений, что может оказаться существенным недостатком при больших окнах сглаживания L. Можно заметить, что текущие значения выходной последовательности y(n), кроме y(0), равны сумме предыдущего шага y(n-1) за вычетом первого элемента этой суммы и добавлением последнего элемента текущей суммы:

 y(n) = y(n-1) - x(n-1) + x(n+L-1). (2)

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

 @y = (0) x $N; for ($n=0; $n<$L; $n++) { $y[0] += $x[$n]; } for ($n=1; $n<$N; $n++) { $m = $n + $L - 1; if($m > $N-1){ last; } $y[$n] = $y[$n-1] - $g[$n-1] + $g[$m]; } for ($n=0; $n<$N; $n++) { $y[$n] /= $L; }

Результаты работы этой программы иллюстрируют графики сигналов на рис.1 и рис.2.

Рис. 1. Пример смеси синусоидальных колебаний и нормального шума
Рис. 1. Пример смеси синусоидальных колебаний и нормального шума
Рис. 1. Пример смеси синусоидальных колебаний и нормального шума

На рис.1 представлена последовательность смеси синусоидальных колебаний и нормального шума g(n):

 for($n=0; $n<$N; $n++) { $x[$n]=0.5*sin(4*$n*$PI/$N) + 0.7*cos(6*$n*$PI/$N) + 0.9*&smgauss($n); }

где внесистемная функция smgauss() используется для вычисления значений шума g(n) с гауссовым распределением. Результат сглаживания шума, полученный с помощью вычисления свертки с прямоугольной функцией показан на рис.2:

Рис. 2. Сглаживание с помощью вычисления свертки с прямоугольной функцией
Рис. 2. Сглаживание с помощью вычисления свертки с прямоугольной функцией
Рис. 2. Сглаживание с помощью вычисления свертки с прямоугольной функцией

При моделировании сигналов, для вычисления гауссова шума используется подпрограмма:

 sub smgauss{ local(*J) = @_; $U = 0; $AN12=12; for($i=$J;$i<$J+$AN12;$i++){ $x=2*(rand(10)/10 - .5); $U +=$x; } return $U/($AN12-1); }.

Подпрограмма smgauss() возвращает значения нормального шума с дисперсией около 0.033 и со средним порядка 10-3 - 10-4, при выборке равной 10000. Исходные данные для расчета нормального шума формируются с помощью встроенной функции rand() генератора шума. Значения переменной $x в подпрограмме smgauss() представляют собой равномерно-распределенный шум с диапазоном значений в интервале: -1 < $x <1.

2. Расчет коэффициентов фильтра нижних частот Баттеруорта

Синусный фильтр нижних частот Баттеруорта является фильтром второго порядка

 y(n) = - A1y(n-1) - A2y(n-2) + B0x(n)

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

 |H(n)|2 = 1/[1 + (sin(π/T)/sin(πBT))2M]. (3)

В этой формуле, символ M означает число полюсов фильтра, T означает интервал времени, с которым производится выборка, символом B обозначена частота отсечки фильтра. Решение, с помощью которого можно определить коэффициенты A1 и A2 этого фильтра, представляется довольно сложным и не является целью настоящей статьи. В работе1 весьма подробно представлены способы для вычисления весовых коэффициентов A1 и A2 и приведена программа на Фортране. Эта программа на Perl имеет следующий вид:

 $pi = 3.1415926539; $FACT = sin($pi*$T*$BW); $F = 1; $M1 = int($M/2); $SECTOR = $pi/$M; $WEDGE = $SECTOR/2; for ($i=0;$i<$M1;$i++) { $B = $FACT*sin($i*$SECTOR + $WEDGE); $C = 1 - $FACT*$FACT; $D = .5*(-$C + sqrt($C*$C + 4*$B*$B)); $E = sqrt($D + 1) + sqrt($D); $G = 2*((2*$B*$B/$D) - 1)/($E*$E); $H = -1/($E**4); $F = $F*(1 - $G - $H); $A1[$i] = -$G; $A2[$i] = -$H; } $B0=$F**(1/$M1);

Результаты расчета, весовых коэффициентов фильтра нижних частот Баттеруорта представлены ниже:

 A1=-1.86338251 A2=0.88656726 B0=0.02222324 A1=-1.72709408 A2=0.74839568 B0=0.02222324

При вычислении были приняты следующие значения констант: порядок фильтра $M=4, период $T=0.005 секунд, частота отсечки $BW=5 герц (в формуле обозначена символом B). Расчет частотной характеристики фильтра, в соответствии с формулой (3) представлен на рис.3:

Рис. 3. Частотная характеристика синусного фильтра нижних частот Баттеруорта
Рис. 3. Частотная характеристика синусного фильтра нижних частот Баттеруорта
Рис. 3. Частотная характеристика синусного фильтра нижних частот Баттеруорта

Отличительная черта фильтров Баттерворта состоит в том, что их передаточная функция является монотонно плоской в полосе пропускания и монотонно спадающей в полосе непропускания, т.е. монотонной по всей частотной области.

3. Фильтры Чебышева и моделирование фильтрации сигналов в частотной области.

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

 |H(n)|2 = 1/[1 + e2Tp(v)2],

где модуль |H(n)| - амплитудно-частотная характеристика фильтра, Tp(v) - полином Чебышева степени p, v = n/nv, nv - частота среза фильтра, e - коэффициент пульсаций. Применением аппроксимации по Чебышеву достигается улучшение аппроксимации идеальной прямоугольной характеристики фильтра нижних частот. Коэффициент e < 1 вводиться для ограничения пульсаций фильтра в полосе пропускания. Чем меньше коэффициент e, тем лучше аппроксимация в полосе пропускания, но одновременно снижается крутизна спада характеристики в области частоты среза nv. Размах пульсаций определяется выражением

 U = 1 - 1/(1 + e2)1/2.

Приемлемый вариант передаточной характеристики достигается подбором коэффициента e и степени полинома p. Полиномы Чебышева имеют вид

 Tp(n) = cos(p)arccos(n), p = 0,1,2,3,... . (4)

Замечательным свойством полиномов Чебышева является то, что при любом целом положительном p, полиномы (4) представляют собой обыкновенные степенные полиномы:

 T0(n) = 1, T1(n) = n, T2(n) = 2n2 - 1, (5) T3(n) = 4n3 - 3n, T4(n) = 8n4 - 8n2 + 1, ... .

Представление полиномов Чебышева в виде (5) может непосредственно использоваться для практической реализации фильтра с помощью программирования на Perl:

 @H = (0) x $N; $N2 = $N/2 + 1; for ($n=0; $n<$N2; $n++) { $U=$n/$CN; $T[0][$n]=$U; $T[1][$k]=2.*$U**2-1; } for ($k=2; $k<$LP; $k++) { for ($n=0; $n<$N2; $n++) { $U=$n/$CN; $T[$k][$n]=2.*$U*$T[$k-1][$n] - $T[$k-2][$n]; } } for ($n=0; $n<$N2; $n++) { $H[$n]=1./(1.+$E2*$T[$LP-1][$n]**2); }

Программный вариант для вычисления квадрата модуля амплитудно-частотной характеристики удобно выполнить в три этапа. На первом - вычислить полиномы степени 0 и 1, на втором определить остальные полиномы T, а на третьем этапе вычислить значения цифрового фильтра H. На рис.4 показан результат расчета цифрового фильтра Чебышева при следующих значениях параметров фильтра: период $N = 64, трансформанта среза $CN = 9, степень полинома $LP = 5, коэффициент пульсаций $E = 0.08.

Рис. 4. Частотная характеристика фильтра нижних частот Чебышева
Рис. 4. Частотная характеристика фильтра нижних частот Чебышева
Рис. 4. Частотная характеристика фильтра нижних частот Чебышева

Для выполнения фильтрации сигнала на рис.1, с помощью цифрового фильтра Чебышева, будем использовать программный модуль CPAN Math::FFT быстрого преобразования Фурье. Для этого, сначала выполним преобразование Фурье исходного сигнала (рис.1), затем результат этого преобразования умножим на передаточную характеристику фильтра (рис.4) и выполним обратное преобразование Фурье:

 use Math::FFT; $N = 64; $N2 = $N/2 + 1; for($i=1;$i<$N2-1;$i++){ $H[$N-$i] = $H[$i]; } use Math::FFT; for($i=0;$i<$N;$i++){ $data->[2*$i]=0.5*sin(4*$i*$PI/$N)+0.7*cos (6*$i*$PI/$N)+0.9*&smgauss($i); $data->[2*$i+1]=0; } $fft = new Math::FFT($data); $coeff = $fft->cdft(); for($i=0;$i<$N;$i++){ $coeff->[2*$i] = $H[$i]*$coeff->[2*$i]; $coeff->[2*$i+1] = $H[$i]*$coeff->[2*$i+1]; } $result = $fft->invcdft($coeff); for($i=0;$i<$N;$i++){ $B[$i]=$result->[2*$i]; }.

В четвертой строке, приведенного выше кода, выполняется формирование фильтра на весь период N. В последнем цикле for() нам потребуется только вещественная составляющая выходной последовательности $result->[2*$i]. На рис.5 показан результат низкочастотной фильтрации с помощью фильтра Чебышева с частотной характеристикой, представленной на рис.4.

Рис. 5. Результат низкочастотной фильтрации с помощью фильтра Чебышева
Рис. 5. Результат низкочастотной фильтрации с помощью фильтра Чебышева
Рис. 5. Результат низкочастотной фильтрации с помощью фильтра Чебышева

4. Согласованная фильтрация сигналов.

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

 y(n) = s(k)s(k+n).

Импульсный отклик фильтра будет иметь вид h(n)=s(k-n), а частотная характеристика H(n)=F*(n), где F*(n) - комплексно-сопряженная функция преобразования Фурье искомого сигнала. Для оптимального фильтра, если требуется обнаружить сигнал s(n) в шумах g(n):

 y(n) = s(n) + g(n), (6)

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

 H(n) = F*(n)/G(n),

где G(n) - спектральная плотность мощности шума. Программная реализация согласованной фильтрации будет состоять из следующих шагов:

  • создание эталона сигнала s(n);
  • преобразование Фурье эталона s(n) для получения комплексного фильтра F(n);
  • преобразование Фурье входной последовательности y(n) для получения комплексного спектра Y(n);
  • перемножение спектра Y(n) и комплексно-сопряженного фильтра F*(n);
  • вычисление обратного преобразования Фурье от произведения Y(n)F*(n).

Результат вычислений B(n) будет представлять собой оценку взаимной корреляции сигналов s(n) и y(n). Последнее обстоятельство позволяет значительно упростить программный код, если воспользоваться подпрограммой correl() из модуля CPAN Math::FFT, рассмотренного в первой статье цикла. Программная реализация для обнаружения гауссова сигнала s(n) представлена ниже:

 use Math::FFT; $N = 64; $N2 = $N/2 + 1; $sigma = 1.44; $a = 2.*$sigma*$sigma; for ($i=0; $i<$N2; $i++) { $s[$i] = exp(-$i*$i/$a); } for($i=1;$i<$N2-1;$i++) { $s[$N-$i] = $s[$i]; } for($i=0;$i<$N;$i++) { $signal->[2*$i] = $s[$i]; $signal->[2*$i+1]=0; } $cm = 20; for ($i=0; $i<$N; $i++) { $t=$i - $cm; $data->[2*$i] = exp(-$t*$t/$a) + &smgauss($i) + 1; $data->[2*$i+1] = 0; } $fft = new Math::FFT($data); $corr = $fft->correl($signal);

Выходной результат возвращается в виде круговой корреляции $corr (ссылка на массив @corr) с задержками во времени. Для получения нулевой задержки, центр гауссова импульса совмещается с нулевым отсчетом, как это показано на рис.6:

Рис. 6. Формирование эталонного сигнала для получения нулевой задержки
Рис. 6. Формирование эталонного сигнала для получения нулевой задержки
Рис. 6. Формирование эталонного сигнала для получения нулевой задержки

На рис.7 показана случайная последовательность сигнала, представляющая собой аддитивную смесь искомого импульса (с произвольной задержкой) и нормального шума:

Рис. 7. Пример формирования сигнала в шумах на входе согласованного фильтра
Рис. 7. Пример формирования сигнала в шумах на входе согласованного фильтра
Рис. 7. Пример формирования сигнала в шумах на входе согласованного фильтра

На рис.8 показан результат согласованной фильтрации.

Рис. 8. Результирующий сигнал на выходе согласованного фильтра
Рис. 8. Результирующий сигнал на выходе согласованного фильтра
Рис. 8. Результирующий сигнал на выходе согласованного фильтра

5. Оценка спектров мощности сигналов.

Необходимость оценки спектра мощности возникает во многих приложениях, например при создании оптимальных линейных фильтров, включая оптимальный согласованный фильтр рассмотренный выше. Результат преобразования Фурье, примененный к автокорреляционной функции случайной последовательности y(n), называется спектральной плотность мощности функции y(n). Спектр мощности можно вычислить непосредственно как квадрат комплексного спектра, как было показано в первой статье цикла. Так же можно воспользоваться подпрограммой spctrm() модуля CPAN Math::FFT:

 $fft = new Math::FFT($signal); $spectrum = $fft->spctrm;

Первое значение массива $spectrum необходимо удвоить. Для построения графика спектра мощности целесообразно создать дополнительный массив P:

 for($i=0;$i<$N;$i++) { $P[$i] = $spectrum->[$i]; } $P[0] *=2;

На рис.9 показан график спектральной плотности мощности функции y(n), в соответствии с формулой (6).

Рис. 9. Спектр мощности случайной последовательности y(n)
Рис. 9. Спектр мощности случайной последовательности y(n)
Рис. 9. Спектр мощности случайной последовательности y(n)

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

Выводы

Представленные примеры моделирования методов цифровой фильтрации сигналов показали широкие возможности применения математических библиотек Perl. Использование подпрограмм модуля CPAN Math::FFT позволяет в значительной степени упростить программный код и создавать компактные программные решения.

Литература

1 Р. Отнес, Л. Эноксон. Прикладной анализ временных рядов.

2 А. Оппенгейм, Р. Шафер. Цифровая обработка сигналов.

3 Р. Дуда, П. Харт. Распознавание образов и анализ сцен.


Ресурсы для скачивания


Комментарии

Войдите или зарегистрируйтесь для того чтобы оставлять комментарии или подписаться на них.

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=40
Zone=Linux, Open source
ArticleID=517931
ArticleTitle=Моделирование методов цифровой обработки сигналов на языке Perl: Часть 2. Использование математических библиотек Perl для фильтрации сигналов
publish-date=09092010