Глава 2. Специальные методы

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

• Гармонический анализ

• Фильтрация

• Амплитудно-частотная характеристика

2.1. Гармонический анализ

Основным инструментом современного гармонического анализа является преобразование Фурье. В MATLAB'е этому соответствуют функции fft, ifft, fft2, ifft2, fftn и ifftn. Напомним, что по определению, дискретное преобразование Фурье последовательности вычисляется по формуле:

Вот перевод справочного экрана, получаемого при вызове help fft:

FFT Discrete Fourier transform.

Дискретное преобразование Фурье ДПФ

FFT(X) is the discrete Fourier transform (DFT) of vector X. If the length of X is a power of two, a fast radix-2 fast-Fourier transform algorithm is used. If the length of X is not a power of two, a slower non-power-of-two algorithm is employed. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton dimension.

FFT(X) - дискретное преобразование Фурье (DFT) вектора X. Если длина X – степень числа 2, то применяется быстрый алгоритм преобразования Фурье. Если длина X - не степень числа 2, более медленный алгоритм используется. Для матриц, операция FFT применяется к каждому столбцу. Для многомерных массивов операция FFT оперирует на первом нетривиальном измерении.

**********************************

FFT(X,N) is the N-point FFT, padded with zeros if X has less than N points and truncated if it has more.

FFT(X,N) - N-точечное ДПФ, дополняемое нулями, если X имеет меньше чем N точкек, и обрезанное, если X имеет больше отсчетов.

FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM.

FFT(X,[],DIM) или FFT(X,N,DIM) применяет операцию FFT к размерности DIM.

For length N input vector x, the DFT is a length N vector X, with elements

Для входного вектора x длины N, X=DFT(x) имеет N комплексных координат

Built-in function

Встроенная функция системы

**********************************

Замечание 1. "FFT" – аббревиатура от " Fast Fourie Transform". Соответствующие русские эквиваленты: "БПФ" - от " Быстрое Преобразование Фурье ". Успех вычислительной процедуры является именно следствием применения специального, "быстрого" варианта алгоритма.

Замечание 2. Как отмечено в [Ахо и Ульман], вектор fq=fft(p,n), представляет собой массив значений полинома с коэффициентами, составленными из координат вектора p, полученных при подстановке в него чисел .

Обратным к дискретному преобразованию Фурье последовательности {X(n), n=0,..,N-1} является вектор x, координаты которого вычисляются по формуле:

Функция ifft используется при обращении ДПФ и в коде содержит обращение к встроенной функции fft:

**************************************

function y = ifft(x,nfft,dim)

if nargin == 3 % Выполняем операцию над размернрстью dim

y = conj(fft(conj(x),nfft,dim));

elseif nargin == 2 %

y = conj(fft(conj(x),nfft)); % Работаем с первой нетривиальной размерностью

% с данными длиной nfft

else

y = conj(fft(conj(x))); % Работа без модификации длины блока

end

% Присвоим переменной dim номер первой нетривиальной размерности

if nargin<3,

dim = min(find(size(x)>1));

if isempty(dim), dim = 1; end

end

% Посчитаем число данных в обрабатываемом блоке

m = size(y,dim);

y = y/m;

**************************************

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

Замечание 4. Чем длиннее , тем более частой будет сетка для оси частот в приближении непрерывного преобразования Фурье с помощью дискретных данных. Однако, следует иметь в виду эффект симметрии ДПФ в случае реального входного вектора. Он состоит в том, что, если все x(i) - вещественны, то X(N-i)=conj(X(i)), где conj(a+i*b)=a-i*b. Поэтому имеет смысл рассматривать только первые N/2 точек ДПФ.

 

2.2. Встроенная функция filter

Выполнение операций с данными "на-лету", т.е. по мере поступления в порт ЭВМ последовательных блоков данных требует строгой оптимизации вычислительных алгоритмов по скорости выполнения и памяти. Можно использовать преобразования Фурье, требующее в своем быстром варианте C*N*lnN операций с блоком длины N (C – постоянная, не зависящая от N). Однако, при решении задач обработки данных с известными спектральными характеристиками выгоднее использовать фильтры. В зависимости от того, какой идеальный фильтр используется при проектировании фильтра и какими дополнительными свойстваим должен обладать полученнй реальный фильтр, используют средства MATLAB'a для проектирования (см. раздел 2.4.) и визуализации характеристик фильтра (см. раздел 2.5). Для стендовых испытаний разработанных фильтров используется filter – встроенная функция, реализация которой на Вашем любимом языке может быть выполнена на основе следующего описания алгоритма. Значение n-ой компоненты y(n) выходного вектора y происходит согласно схеме:

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

Вызов Y = FILTER(B,A,X) решает разностное уравнение

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

Если a(1) не равно 1, FILTER нормализует фильтр делением на a(1).

Когда X - матрица, FILTER выполняет операции по колоннам X.

Вызов [Y,Zf] = FILTER(B,A,X,Zi) создает доступ к начальным значениям Zi и конечным значениям Zf - задержек. Zi это вектор длины MAX(LENGTH(A),LENGTH(B))-1 или массив таких векторов для каждой колонки X.

 

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

Функции MATLABbutter, bessel, cheby1, cheby2, ellip адвтоматически проектируют фильтр, удовлетворяющие требованиям, спецификаций, сформулирпванных в терминах зависимости коэффициента усиления соответствующей линейной системы от частотыю Так опредекляются низкочастотные, высокочастотные, полоснозаграждающие ( lowpass, highpass, bandstop), и полосовые (bandpass) фильтры, Как они могут быть разработаны, как показано на примерах. Частоты среза должны быть определена так, чтобы находились в пределах [0, 1], причем 1 соответствует частоте Найквиста fc=1/(2T). Чтобы спроектировать цифровой низкочастотный фильтр Баттерворта, используйте команду:

[b,a] = butter(n,Omegac)

где n - число полюсов и нулей в передаточной функции фильтра, и Omegac - нормированная цифровая частота среза. Чтобы спроектировать высокочастотный фильтр с частотой среза Omegac, используйте команду

[b,a] = butter(n,Omegac,'high')

Чтобы спроектировать полосовой фильтр с полосой пропускания от Omega1 до Omega2, определите Omega = [Omega1, Omega2] и используйте команду

[b,a] = butter(n,Omega)

Чтобы спроектировать полоснозаграждающий фильтр Omega1 до Omega2 определите как выше и используйте команду

[b,an] = butter(n,Omega,'stop')

Проектирование чебышевский фильтр n-го порядка первого типа, выполнен, используя те же самые методы что butter, за исключением того, что "butter" заменено на "cheby1".

Примеры обращения к описанным функциям рассмотрены в Главах 4-5.

 

2.4. Амплитудно-частотные характеристики фильтра.

Чтобы наглядно представлять воздействие на сигнал линейной системы, реализованной с помощью функции filter, мы будем использовать функцию FREQZ – (requency response of digital filters) частотный отклик цифрового фильтра, уделив внимания как способам обращения, так и реализации.

Описание

[h,w]=freqz(b,a,n) возвращает комплексный частотный отклик H(ej)

цифрового фильтра, заданного векторами b,a , вычисленный в n равномерно расположенных точках на отрезке между 0 и , возвращаемых в векторе w. Лучше, хотя и не обязательно, выбрать для n точную степень 2, что позволит применить БПФ-алгоритм. По умолчанию n=512.

[h,f] = freqz(b,a,n,Fs) специфицируеи отсчетность Fs, Гц. По умолчанию Fs=1. Тогда вектор f содержит актуальные частоты между 0 и Fs/2 (частотой Найквиста) в которых будет вычислен частотный отклик.

h = freqz(b,a,w) возвращает частотный отклик в точках, указанных в векторе w (измеренных в рад/отсчет).

freqz(...) без выходных аргументов строит АЧХ – амплитудно-частотную характеристику и ФЧХ – фазо-частотную характеристику на указанном частотном диапазоне в текущем окне.

Пример 1. Вызов

[bu1,au1]=freqz([1 2], [2 -1], 3)

приводит к результату

bu1 =

3.0000

0.5000 - 1.4434i

-0.2143 - 0.6186i

au1 =

0

1.0472

2.0944

Пример 2. Используем значения, устанавливаемые по умолчанию, чтобы посмотреть на графики:

freqz([1 2], [2 -1])

Алгоритм

freqz использует алгоритм преобразования Фурье (вызывая fft), когда в вызове присутствует размер n массива возвращаемых значения. Тогда используется то обстоятельства, что вектор fq=fft(p,n), как отмечено в [1], стр 123, представляет собой массив значений полинома с коэффициентами, составленными из координат вектора p, полученных при подстановке в него чисел . Так что достаточно дополнить векторы a и b нулями до нужного числа. Тогда частотный отклик h получится как h = fft(b,n)./fft(a,n)

Отсюда желательность, чтобы n было степенью 2. Когда же в вызове присутствует частотный вектор w, или n меньше чем max(length(b),length(a)), freqz вычисляет значения полинома, используя схему Горнера, обращаясь к функции polyval.

 

Оглавление