Специальные методы
Этот раздел расскажет вам больше о работе с данными, представляющими временные ряды, уделяя особое внимание следующим алгоритмам.
-
Гармонический анализ
-
Фильтрация
-
Амплитудно-частотная характеристика
Основным инструментом современного гармонического анализа является преобразование Фурье. В MATLAB'е этому соответствуют функции fft, ifft, fft2, ifft2, fftn и ifftn. Напомним, что по определению, дискретное преобразование Фурье F(x)={X(k),0<=k<=N-1} последовательности
x={x(k),0<=k<=N-1} вычисляется по формуле:
Вот перевод справочного экрана, получаемого при вызове 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. В некоторых руководствах принято деление на длину входной последовательности относить на прямое преобразование и даже распределять поровну - каждому по 1/(корень N).
Замечание 4. Чем длиннее, тем более частой будет сетка для оси частот в приближении непрерывного преобразования Фурье с помощью дискретных данных. Однако, следует иметь в виду эффект симметрии ДПФ в случае реального входного вектора. Он состоит в том, что, если все x(i) - вещественны, то X(N-i)=conj(x={x(k),0<=k<=N-1}X(i)), где conj(a+i*b)=a-i*b. Поэтому имеет смысл рассматривать только первые N/2 точек ДПФ.
Выполнение операций с данными "на-лету", т.е. по мере поступления в порт ЭВМ последовательных блоков данных требует строгой оптимизации вычислительных алгоритмов по скорости выполнения и памяти. Можно использовать преобразования Фурье, требующее в своем быстром варианте 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.
Функции MATLAB'а butter, 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.
Чтобы наглядно представлять воздействие на сигнал линейной системы, реализованной с помощью функции 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, полученных при подстановке в него чисел exp(-2Пkn/N*j),0<=k<=N-1. Так что достаточно дополнить векторы a и b нулями до нужного числа. Тогда частотный отклик h получится как h = fft(b,n)./fft(a,n)
Отсюда желательность, чтобы n было степенью 2. Когда же в вызове присутствует частотный вектор w, или n меньше чем max(length(b),length(a)), freqz вычисляет значения полинома, используя схему Горнера, обращаясь к функции polyval.