Содержание


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

Часть 1. Моделирование и исследование сигналов на языке Perl

Comments

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

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

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

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

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

1. Роль Perl и архива CPAN в задачах компьютерного моделирования

Архив свободных программных ресурсов CPAN является огромным хранилищем ресурсов Perl и позволяет программистам найти готовые решения практически в любой области. Главным элементом CPAN являются модули, содержащие программы для решения конкретных прикладных задач, включая научные исследования и моделирование процессов. Привлекательность использования Perl в задачах математического моделирования объясняется просто - Perl унаследовал математические библиотеки языка C.

Модульность Perl обеспечивается делением программ на независимые фрагменты (пакеты), каждый из которых имеет собственное пространство имен. Это позволяет не беспокоиться о конфликтах между пакетами и вызывающей частью программы. Многократный вызов подпрограмм осуществляется с помощью модулей. Модули представляют собой пакеты, оформленные определенным образом в файлах с расширение pm. Имена файлов совпадают с именами модулей в первой строке: package NAME.

2. Назначение и состав модуля Math::FFT

Значительная доля анализа и моделирования сигналов связана с преобразованием Фурье. Алгоритм БПФ позволяет проводить вычисления в различных задачах наиболее быстро и эффективно. В этих задачах преобразование Фурье может использоваться как инструмент в определении корреляционных и передаточных функций, для вычисления функций фильтров, вычисления спектра мощности сигнала, интерполяции сигналов при моделировании аналого-цифровых преобразований в системах цифровой обработки данных. Настоящий раздел содержит краткий обзор одной из реализаций БПФ в модуле Math::FFT в виде объектно-ориентированной конструкции. В модуле реализован ряд функций для вычисления БПФ одномерной последовательности отсчетов размерностью N и кратной степени 2. В этих функциях используются алгоритмы вычисления БПФ на языке C, представленных в файле fft4g.c модуля Math::FFT.

В состав модуля входят такие функции как:

  • cdft: комплексное дискретное преобразование Фурье;
  • ddct: дискретное косинусное преобразование;
  • ddst: дискретное синусное преобразование;
  • correl: вычисление корреляции между двумя сигналами;
  • spctrm: вычисление спектра мощности,

а так же вычисление свертки двух сигналов и других полезных характеристик - вычисление среднего, средне-квадратичного отклонения, медианы и др.

Входная последовательность данных, представленная массивом $data, используется при создании объекта класса Math::FFT:

                     my $fft = new Math::FFT($data),

где $fft метод объекта. Вызов метода $fft при вычислении комплексного БПФ для набора данных x[j] задается конструкцией:

                     $coeff = $fft->cdft().

Здесь $data - ссылка на набор данных data[0... 2*N-1]:

                     data[2*j] = Re(x[j]),
                     data[2*j+1] = Im(x[j]), 0<=j<N.

Упаковка массива данных осуществляется таким образом, что вещественная часть комплексного числа размещается в четных элементах массива, а мнимая часть в нечетных элементах, соответственно. Код метода для вычисления обратного комплексного БПФ записывается следующим образом:

                     $orig_data = $fft->invcdft([$coeff]),

где $orig_data - восстановленная последовательность сигнала, а $coeff является ссылкой на массив коэффициентов Фурье coeff[0...2*N-1]:

                     coeff[2*j] = Re(x[j]),
                     coeff[2*j+1] = Im(x[j]), 0>=j<N.

Остальные методы, реализующие БПФ в модуле Math::FFT записываются аналогично рассмотренному. Лаконично записывается метод вычисления корреляции двух сигналов:

                      $fft1 = new Math::FFT($data1)
                      $fft2 = new Math::FFT($data2)
                      $corr = $fft1->correl($fft2),

где $fft1 и $fft2 - результаты преобразования Фурье сигналов $data1 и $data2 соответственно, а $corr - корреляция. Выходной массив $corr возвращается в виде круговой корреляции с положительными задержками в диапазоне отсчетов от $corr->[0] (нулевая задержка) до $corr->[$n/2-1], и с отрицательными задержками для отсчетов от $corr->[$n-1] до $corr->[$n/2]. Если сигнал $data2 задержан относительно $data1, то корреляция $corr образует пик в положительном диапазоне задержек.

Частные случаи использования методов Math::FFT удобнее показать на конкретных примерах при моделировании временных сигналов - синусоидального сигнала, импульсного прямоугольного сигнала и функции sinc, гауссова сигнала.

3. Примеры моделирования и преобразования колебаний и сигналов с помощью модуля Math::FFT

3.1. Синусоидальный сигнал

Синусоидальные колебания на Perl реализуются с помощью встроенных математических функций sin() и cos(). Пример кода для формирования смеси (для общности) синусоидальных колебаний выглядит следующим образом:

   for($k=0; $k<$N; $k++) {
   $signal->[$k] = 0.5*sin(4*$k*$PI/$N) + 0.7*cos(6*$k*$PI/$N);
      },

где $N - размерность дискретной последовательности сигнала, $PI=3.1415926539. График функции сигнала $signal представлен на рис.1 и построен с помощью подпрограммы rasp2() для символьно-графического представления сигналов:

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

Спектр смеси синусоидальных колебаний $signal представлен на рис.2 и получен с помощью дискретного преобразования Фурье:

                     $fft = new Math::FFT($signal);
                     $coeff = $fft->rdft();,

где функция rdft() выполняет вычисление реальной части дискретного преобразования Фурье, а $coeff – коэффициенты ряда Фурье.

Рис. 2. Спектр Фурье смеси синусоидальных колебаний.
Рис. 2. Спектр Фурье смеси синусоидальных колебаний.
Рис. 2. Спектр Фурье смеси синусоидальных колебаний.

Следующий пример позволит убедиться в правильности вычислений алгоритма БПФ с помощью Math::FFT. Известно, что спектр синусоиды представляет собой сумму двух дельта-функций. Выполним обратную процедуру - вычислим обратное комплексное преобразование Фурье суммы дельта-функций. Код для формирования дельта функции выглядит следующим образом:

                     for($i=0;$i<$N;$i++){
                                        $coeff->[2*$i]=0;
                                        $coeff->[2*$i+1]=0;
                     }
                     $coeff->[2*3]=$N/2;
                     $coeff->[2*$N-2*3]=$N/2;

где [3] номер трансформанты (гармоники) дискретной последовательности косинуса. Четные трансформанты [2*$i] соответствуют реальной части комплексного числа, а нечетные - мнимой части. Размерность комплексного массива составляет 2N значений. Выполним обратное комплексное преобразование Фурье массива коэффициентов $coeff:

                     $fft = new Math::FFT($coeff);					
                     $cos_from_delta = $fft->invcdft($coeff);
                     for($i=0;$i<$N;$i++){					
                                      $g[$i]=$cos_from_delta->[2*$i];    
                     }
                     &rasp2(\30,\@g,\$N);

Здесь функция invcdft() является обратным комплексным дискретным Фурье преобразованием, массив @g - выходной массив данных, представленных на рис.3. Функция rasp2() обеспечивает символьно-графическое представление одномерных сигналов.

Рис. 3. Результат преобразования суммы двух дельта-функций.
Рис. 3. Результат преобразования суммы двух дельта-функций.
Рис. 3. Результат преобразования суммы двух дельта-функций.

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

3.2. Прямоугольный импульс и функция sinc

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

                     for($i=0;$i<$N;$i++){
                                        $data->[2*$i]=0;
                                        $data->[2*$i+1]=0;
                      }
                     for($i=0;$i<$size;$i++){
                                        $data->[2*$i]=1;
                     }
                     for($i=0;$i<$size-1;$i++){
                                        $data->[2*$N-2-2*$i]=1;
                     }

где длительность (число отсчетов) импульса составляет $L=2*$size-1. При длительности $L=5, график сигнала будет выглядеть следующим образом:

Рис. 4. Пример построения симметричного прямоугольного импульса.
Рис. 4. Пример построения симметричного прямоугольного импульса.
Рис. 4. Пример построения симметричного прямоугольного импульса.

Комплексное преобразование Фурье задается кодом:

                     $fft = new Math::FFT($data);					
                     $coeff = $fft->cdft();

Преобразование Фурье симметричного сигнала будет содержать только вещественную составляющую, то есть только отсчеты coeff[2*j]=Re(x[j]). Мнимая часть coeff[2*j+1]=Im(x[j]) комплексной функции спектра сигнала в этом случае будет равна нулю. Для графического представление спектра будет достаточно учитывать только вещественную компоненту:

                     for($i=0;$i<$N;$i++){					
                                      $g[$i]=$coeff->[2*$i];    
                     }
                     &rasp2(\25,\@g,\$N);.
Рис. 5. Спектр симметричного прямоугольного импульса.
Рис. 5. Спектр симметричного прямоугольного импульса.
Рис. 5. Спектр симметричного прямоугольного импульса.

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

                     $Tx=9;
                     for($i=0;$i<$L;$i++){ 
                                          $arg=$Tx*$i*$PI/$L;
                                          $g[$i]=&dsinc(\$arg);
                     }
                     &rasp2(\25,\@g,\$L);,

где параметр $Tx соответствует половине длительности прямоугольного импульса, $L - размерность последовательности дискретных данных, а функция dsinc() вычисляется с помощью подпрограммы:

                     sub  dsinc{
                                local(*D)=@_;
                                my $scc=1.;
                                if($D!=0.) { $scc=sin($D)/($D); }
                                return $scc;
                     }.

На рис.6 представлен график функции sinc(t) с размерностью последовательности отсчетов равной 35.

Рис. 6. График функции sinc(t).
Рис. 6. График функции sinc(t).
Рис. 6. График функции sinc(t).

3.3. Гауссов сигнал

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

                     $sigma=1.22; 
                     $a=2.*$sigma*$sigma;
                     $sht=$N/2-1;
                     for($i=0;$i<$N;$i++){ 
                                        $t=$i - $sht; 
                                        $g[$i]=exp(-$t*$t/$a);
                     },

где параметр $sigma эквивалентен половине длительности гауссова импульса, определяемой на уровне 0,606 от амплитуды импульса; параметр $sht обеспечивает смещение импульса. График синтезированного сигнала представлен на рис. 7.

Рис. 7. Гауссов импульс со смещением во временной области.
Рис. 7. Гауссов импульс со смещением во временной области.
Рис. 7. Гауссов импульс со смещением во временной области.

После выполнения комплексного преобразования Фурье, спектр сигнала определяет как спектр мощности гауссова импульса:

    for($i=0;$i<$N;$i++){					
    $g[$i]=$coeff->[2*$i]*$coeff->[2*$i] + $coeff->[2*$i+1]*$coeff->[2*$i+1]; 
      }

График спектра мощности гауссова импульса показан на рис.8.

Рис. 8. Спектр мощности гауссова импульса.
Рис. 8. Спектр мощности гауссова импульса.
Рис. 8. Спектр мощности гауссова импульса.

4. Конструирование собственных модулей в Perl

Как отмечалось выше, модули это пакеты, оформленные таким образом, что имя файла модуля совпадает с именем в объявлении пакета в первой строке:

                     package mymodul.pm;
                     BEGIN{
                     require Exporter;  
                     @ISA= qw(Exporter);  
                     @EXPORT=qw(&rasp2 &dsinc);
                     }
                     sub rasp2{
                     local(*L,*B,*N)=@_;
                     $AMAX = -65536;
                     $AMIN =  65536;
                     for ($i=0;$i<$N;$i++){
                                        if($B[$i]< $AMIN){$AMIN=$B[$i];}
                                        if($B[$i]> $AMAX){$AMAX=$B[$i];}
        	         }
                     for ($i=0;$i<$N;$i++){$X[$i]=$B[$i]-$AMIN;}
                     $AM = $AMAX-$AMIN;
                     $AMAX = sprintf "%12.3f",$AMAX;
                     $AMIN = sprintf "%14.5f",$AMIN;
                     print "MAX=$AMAX MIN=$AMIN\n";        
                     if($AM == 0){$AM=1;}
                     @hr = (".") x ($L+4);
                     print "@hr\n";
                     for ($i=0;$i<$N;$i++){$X[$i]=$L*$X[$i]/$AM;}
                     for ($j=0;$j<$N;$j++){
                                               $bt = sprintf "%14.5f",$B[$j];
                                               $num = sprintf "%4d",$j+1;
                                               @H3 = (" ") x $L;
                     		           $K=int($X[$j]+.5);
                                               if($K > $L-1){$K=$L-1;}
                                               $H3[$K]="*";
                                               print "$num|@H3|$bt\n"; 
                     }
                     print "@hr\n";
                     }
                     sub  dsinc{
                     local(*D)=@_;
                     my $scc=1.;
                     if($D!=0.) { $scc=sin($D)/($D); }
                     return $scc;
                     }
                     return 1;
                     END {  }

В строке 3 загружается модуль Exporter, который управляет внешними интерфейсами модуля. Строка 4 инициализирует специальный массив @ISA, который существует на уровне пакетов. В строке 5 список процедур (подпрограмм) заноситься в массив @EXPORT, который так же существует на уровне пакетов. Далее идут тексты собственных подпрограмм rasp2() и dsinc(), которые использовались в примерах по моделированию.

В заключение приведем обобщенный пример загрузки собственного модуля mymodul.pm и модуля Math::FFT для вычисления спектра Фурье входной последовательности данных:

                     #!/usr/bin/perl
                     print "Content-type: text/plain\n\n";
                     use Math::FFT;
                     use mymodul;
                     $PI = 3.1415926539;
                     my $N = 32;
                     for($k=0; $k<$N; $k++) {
   $g[$k] = 0.5*sin(4*$k*$PI/$N) + 0.7*cos(6*$k*$PI/$N);
        $series->[$k] = $g[$k];
                     }
                     &rasp2(\25,\@g,\$N);
                     my $fft = new Math::FFT($series);
                     my $coeff = $fft->rdft();
                     for (my $k=0; $k<$N; $k++) {
                                                   $g[$k] = $coeff->[$k];
                     }
                     &rasp2(\25,\@g,\$N);
                     exit;

5. Выводы

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


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


Комментарии

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

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