ЦОС. Пытаюсь понять принцип цифровой фильтрации сигнала. Опции

[b_of_b]

ЦОС. Пытаюсь понять принцип цифровой фильтрации сигнала. Подскажите пожалуйста правильна ли в принципе моя логика.
Код на 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 и с ЦОС работаю впервые
ЗЗЫ: Если есть какая-либо статься написанная по такой схеме (как создать, проанализировать и увидеть ) подскажите пожалуйста.

[Vasiliy Rufitski...]

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

Формула, описывающая гармоническое незатухающее колебание:
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 Гц, с нужной частотой среза,
которая должна быть в диапазоне фильтруемого сигнала, иначе после фильтрации получится почти то же что и до фильтрации.

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

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

[Grt]

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

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

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