ЦОС. Пытаюсь понять принцип цифровой фильтрации сигнала. Подскажите пожалуйста правильна ли в принципе моя логика.
Код на MATLAB

1. Создадим сигнал
for q=1:1:1024;
mass(q)=sin(2*pi/512*q); % 512 Гц
mass(q)=mass(q)+sin(2*pi/256*q);% 256 Гц
mass(q)=mass(q)+sin(2*pi/128*q);% 128 Гц
mass(q)=mass(q)+sin(2*pi/64*q); % 64 Гц
mass(q)=mass(q)+sin(2*pi/32*q); % 32 Гц
mass(q)=mass(q)+4;
end;
plot (mass);

2. Проанализируем сигнал с помощью ДПФ
dp=fft(mass); % ДПФ
Am=sqrt(real(dp).^2 + imag(dp).^2)/512; % Амплитуда
%(Первая из них видимо – постоянная, а дальше по частотам ? )
plot(Am(1:256));
%title('Frequency content of y'); xlabel('frequency (Hz)');
Как корректно отобразить амплитудно-частотный график ?

3. Создадим произвольный фильтр с помощью fdatool (лишь бы наглядно давил какой-то диапазон)
Через экспортирование (File->Export…) получаем коэффициенты фильтра:
FiL = [
1.0000 -0.0000 -1.0000 1.0000 -0.7044 0.7096 1.0000 1.6212 1.0000 1.0000 -0.0894 0.7528
1.0000 -1.9288 1.0000 1.0000 -1.2281 0.8118 1.0000 1.1771 1.0000 1.0000 0.3301 0.8700
1.0000 -1.8281 1.0000 1.0000 -1.5012 0.9159 1.0000 0.9892 1.0000 1.0000 0.5205 0.9619
1.0000 -1.7785 1.0000 1.0000 -1.6154 0.9770
];
for q=1:1:42; FiL2(q)=FiL(q); end;
for q=43:1:1024; FiL2(q)=0; end;
Здесь я расширяю размерность фильтра до размерности сигнала – это нужно ?

4. Отфильтруем сигнал - сворачиваем сигнал с фильтром
w = conv(mass,FiL2);
dp2=fft(w); % ДПФ
Am2=sqrt(real(dp2).^2 + imag(dp2).^2)/512;
plot(Am2(2:50));
clear d Am dp mass w FiL dp2 Am2;

I. Смысл моего примера (вопроса) в том, что я сгенерировал, как мне кажется, чистейший сигнал: частота дискретизации кратная частотам заложенным в сигнале. Как я думал на амплитудно-частотном графике должен был бы получиться "аккуратный частокол" с единичными амплитудами (а у меня первый коэффициент куда-то ввысь улетает постоянно). Да и формулу нахождения амплитуд я «подгонял под ответ». Затем я собирался создать фильтр (какие коэффициенты выдает тулза Матлаба fdatool я не знаю (кстати в чем (какая прога или тулза) лучше всего фильтры проектировать)), который наглядно «сбивает» какие-то частоты. Коэффициенты фильтра я получил и свернул их с сигналом. Результат стал анализировать с помощью ДПФ – который вывел мне какие-то еще левые частоты, а «сбиваемые» частоты в большинстве случаев «не пострадали». Вот я и спрашиваю: что я не так делаю и может где уже статья есть по данному вопросу (Читать книги типа Айфичера у меня уже нервов не хватает).

II. А первый коэффициент после выполнения ДПФ (для сигнала) является уже частотой или это какая-то постоянная величина не имеющая отношения к спектру, которая как-то отображается на графике (каково ее место на амплитудно-частотном график ) ?

III. Амплитуда высчитывается по такой формуле как у меня (в разных книгах - разные формулы)
Am=sqrt(real(dp).^2 + imag(dp).^2)/512;

Амплитуда = квадратный корень (сумма квадратов действительной и мнимой части ДПФ) / половину частоты дескритизации ?

IV. А почему фильтрацию сигнала не делают по такой схеме: в результате ДПФ получаем значения частот и амплитуд, значит "ненужные" (фильтруемые) частоты можно просто выкинуть, а затем через функцию синуса создать уже отфильтрованный сигнал.

ЗЫ: Я и с MATLAB и с ЦОС работаю впервые
ЗЗЫ: Если есть какая-либо статься написанная по такой схеме (как создать, проанализировать и увидеть ) подскажите пожалуйста.

Уверены, что частота 512 и 256 Гц?

1. Про герцы речь даже не идет. Есть только частота дискретизации в герцах, остальное пляшет относительно нее. Т.е. Ваши 1024 отсчетов могли бы быть получены и с fd=1000 Гц, и с fd=10000000 Гц. Общий вид спектра (при такой генерации сигнала) не изменился бы. Еще Альберт в своей теории доказал,что все относительно.
2. Почему у Вас за период наблюдения (1024 отсчета) влазит 2 периода синуса с частотой 512 Гц, и целых 4 периода синуса с частотой 256 Гц? Поняли в чем смысл? Другими словами, все наоборот: то, что вы считатете 512 Гц - на самом деле fd/512 и т.д.

Далее по тексту: первый коэффициент это не частота (!), как в прочем и остальные коэффициенты. Коэффициенты fft - амплитуда (в вашем случае) определенной частотной составляющей в спектре сигнала. Т.е. первый коэффициент, говоря по-русски, это амплитуда синуса с частотой 0. На счет графика спектра сигнала - очень много возможностей. В matlab'е очень хороший help. Далее: фильтрацию в частотной области (с помощью fft) делают. Увы, Вы не первый...

Кстати, даже если сигнал одни только синусы, спектр же такого сигнала не всегда одни только красивые "палки". При недосточном частотном разрешении гармоники могут накладываться друг на друга и т.д., и т.п.

Книги? Статьи? Их море. А начинать ЦОС рекомендую с бессмертного творения Рабинера и Голда. Есть на ftp/

В общем случае размерность фильтра и размерность сигнала не совпадают, поскольку для КИХ фильтра свертку не имеет смысла выполнять на большем интервале времени, чем длина его импульсной характеристики.
Фильтрация сигнала при помощи БПФ используеться относительно редко, в основном тогда, когда надо реализовать фильтр с очень узкой полосой расфильтровки и/или со специфичной АЧХ, когда реализация обычным способом приведет к большим требуемым ресурсам, чем реализация алгоритмов прямого и обратного преобразования Фурье
В своей небольшой практике пользовал fdatool, вполне удобный инструмент по синтезу фильтров, однако мне больше понравилась его реализация в седьмой версии матлаба.

Скажите это как-нибудь по-другому. Потомки Вас в такой формулировке не поймут. Убежден, что можно фильтровать "длинный" сигнал "коротким" фильтром.

А я что написал ?

Нужно читать книги, чтобы понимать самому что делаешь. Вам на форуме могут дать только совет, не факт что верный или оптимальный, но забесплатно Вашу работу никто делать не будет...


ЯвляеЦЦо, но он имеет смысл только в своем частотном окне...


Если под амплитудой Вы подразумеваете амплитудный спектр то почти... причем там частота дискретизации? То в знаменателе длинна БПФ. Иногда ее опускают.



делают и так, просто иногда требования к АЧХ/ФЧХ или импульсной/переходной характеристике не реализуются БПФ на данной платформе.

Насколько я представляю, на входе у Вас сумма гармоник.

Формула, описывающая гармоническое незатухающее колебание:
A*sin(2*pi*f*t), в случае сэмплированного сигнала:
A*sin(2*pi*f*i/fs), где i - номер отсчёта сэмплированного сигнала.
Так что

sin(2*pi/512*q); и
sin(2*pi/256*q);

частотно соотносятнся вовсе не как 512Гц к 256 Гц, а наоборот...

Попробуем привести всё к единой частоте дискретизации:

mass(q)=sin(1*2*pi/512*q);
mass(q)=mass(q)+sin(2*2*pi/512*q);
mass(q)=mass(q)+sin(4*2*pi/512*q);
mass(q)=mass(q)+sin(8*2*pi/512*q);
mass(q)=mass(q)+sin(16*2*pi/512*q);

Имеем ряд из 1, 2, 4, 8 и 16 Гц, при частоте дискретизации 512 Гц
Ах да, и ещё 0 Гц:
mass(q)=mass(q)+4;

Что мы и имеем возможность лицезреть на графике результата fft.

Стоит отметить, что модуль fft отображает "зеркальную часть модуля спектра" в правой части.
С этим лучше в хелп матлаба, там это подробно описано .
В общем, чтобы привести результат в привычный нам спектр от 0 до fs/2 нужно:
Взять только левую половину.
Сместить все составляющие результата на 1 отсчёт влево.
Амплитуду 0-ой составляющей уполовинить.
Тогда 1 отсчёт будет соответствовать половие герца.

> Как я думал на амплитудно-частотном графике должен был бы получиться "аккуратный частокол" с единичными амплитудами >(а у меня первый коэффициент куда-то ввысь улетает постоянно).
mass(q)=mass(q)+4;

Всё правильно. График зеркальный, амплитуда нулевой частоты удвоена (особенность fft).
Если её упполовинить, то нулевая частота как раз 4 станет.

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

А вообще есть книжка Сергиенко А.Б. Цифровая обработка сигналов. Там всё доступно и с кусками кода из матлаба )

Можно вроде даже из нета закачать.

Или еще есть старая книжка Каппилине К. "Цифровая обработка сигналов" там все детально описано все аспекты.

Вы, наверное, не так дяденьку обозвали! Вот, наверно, то, что требовалось:

Каппелини В., Константинидис А.Дж., Эмилиани Э. Цифровые фильтры и их применение /Пер. с англ. - М.: Энергоатомиздат, 1983. 360 с.
Хорошая книжка. Также очень доходчивот написана и эта книга:

Отнес Р., Эноксон Л. - Прикладной анализ временных рядов. Основные методы.
(http://lib.mexmat.ru/books/13730)

Успехов!

Она есть в электронном виде? Можете дать ссылку?

Посмотрите на dsp-book.narod.ru. Там раньше точно был Константинидис

Дата Oct 23 2007, 16:50

Крайне своевременно. Это для поддержания разговора?