Фильтры в матлабе

[barskuk]

Здаравствуйте.
У меня такой вопрос по матлаб.
У меня есть сигнал, который представляет собой посаженную на пьедистал синусойду с частотой 4,286 МГц, промодулированную пилой.
Этот сигнал необходимо пропустить через некое устройство:

   +--|   1   |-----|   2  |--|     |
__|                           |  4  |------выход
  |                           |     |   
   +--|   3   |---------------|     |

где
1 - полосовой фильтр
2 - ограничитель аплитуды
3 - режекторный фильтр
4 - сумматор

фильтры должны пропускать и задерживать соответственно частоту 4,28 МГц
Помогите мне! Как это реализовать в матлаб!!
Спасибо

[barskuk]

Начнем, а Вы поправьте меня если что не так (я с матлабом до этого не общался, а нужно кое-что для диплома сделать)

close all;
clear all;

     %Подгружаю исходник сигнала который нужно профильтровать
     fid = fopen('c:\MATLAB701\work\signal\0.bin', 'r');

     F = fread(fid);

Добавлено спустя 35 секунд:

теперь расчитываем фильтр. начну с режекторного фильтра. мне необходимо чтобы он вырезал частоту 4,286

function Hd = filter

Fs = 13500;

Fpass1 = 4280;      
Fstop1 = 4286; 
Fstop2 = 4286; 
Fpass2 = 4290; 
Dpass1 = 0.028774368332; 
Dstop  = 0.001;  
Dpass2 = 0.057501127785; 
dens   = 20;  
[N, Fo, Ao, W] = firpmord([Fpass1 Fstop1 Fstop2 Fpass2]/(Fs/2), [1 0 ...
                          1], [Dpass1 Dstop Dpass2]);

b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

[ashik]

Приветствую всех!

У меня возникла следующая проблема:
Мне нужно расчитать коэффициенты IIR фильтров.
При использовании fdatool я получаю не одни и те же результаты, что при получении коэффициентов с коммандной строки.


Например, мне нужны коэффициенты полосового фильтра Чебышева первого типа, 165-235 Hz, 6-го порядка, с неравномерностью в полосе пропускания 3 дБ на частоте 22050Hz:

В окне fdatool я выбираю:

Filter Type - bandpass
Design Method - IIR, Chebywev Type I
Filter order - 6
Frequency Specifications:
Unit - Hz; Fs -22050; Fpass1 - 165; FPass2 - 235;
Magnitude Specifications:
Units - dB; Apass - 3;

коєффициенты преобразовываю в Single Section.
экспортирую, и получаю следующее:

real64_T B[7] = { 2.471229843725e-007, 0,-7.413689531176e-007, 0, 7.413689531176e-007, 0,-2.471229843725e-007};
real64_T A[7] = { 1, -5.97836575023, 14.90172810275, -19.82311671515, 14.84267011882, -5.931073541702, 0.988157816516 };

Это правильные коэффициенты, во всяком случае, у меня свёртка с ними работает правильно.

Когда я делаю, казалось бы, то же самое, но с коммандной строки,

[b,a]=cheby1(3,3,[165 235]/11025)

получаю похожие коэффициенты, но всё же немного отличающиеся, и они не работают:

b =
1.0e-006 *
0.24712242570581 0 -0.74136727711742 0 0.74136727711742 0 -0.24712242570581

a =
1.00000000000000 -5.97836575023038 14.90172810275087 -19.82311671515133 14.84267011881655 -5.93107354170150 0.98815781651604


как видно, А-коеффициенты совпадают, B - нет, хотя очень близки.

Мне нужно вычислять именно с командной строки, т.к. таких фильтров нажно получить около сотни.

Помогите, пожалуйста, разобраться: это проблема точности вычислений, или я просто что-то неверное делаю?
Как мне с коммандной строки получить то же, что и с помощью fdatool?

[Yuri Gendelman]

ashik писал(а):

real64_T B[7] = { 2.471229843725e-007, 0,-7.413689531176e-007, 0, 7.413689531176e-007, 0,-2.471229843725e-007};
real64_T A[7] = { 1, -5.97836575023, 14.90172810275, -19.82311671515, 14.84267011882, -5.931073541702, 0.988157816516 };

b =
1.0e-006 *
0.24712242570581 0 -0.74136727711742 0 0.74136727711742 0 -0.24712242570581

a =
1.00000000000000 -5.97836575023038 14.90172810275087 -19.82311671515133 14.84267011881655 -5.93107354170150 0.98815781651604

У Вас различия в коэффициентах - в 6-м знаке. Мне трудно представить, что такая мизерная разница дает ошибки в расчетах.
При проверке коэффициентов [b, a] Вы не забыли домножить на $10^{-6}$ ВСЕ элементы массива b, а не только первый?

[ashik]

Эта разница в 6-м знаке имеет значение. Фильтр попросту не работает (у меня). коеффициенты фильтров с разными центральными частотами начинают отличаться только в 6-м знаке.
165-235Hz: {0.24712242570581, 0 -0.74136727711742, 0 0.74136727711742, 0 -0.24712242570581},
265-335Hz: {0.24712275214629, 0 -0.74136825643888, 0 0.74136825643888, 0 -0.24712275214629},
365-435Hz: {0.24712313254425, 0 -0.74136939763274, 0 0.74136939763274, 0 -0.24712313254425}.
и т. д.

Спасибо за ответ,
Антон.

Добавлено спустя 1 минуту 52 секунды:

П.с. домножить не забыл

[linder]

Есть слабый квазипериодический звуковой сигнал ~900Гц , который записывается с приемника, на него действуют различные помехи:Например низкочастотный шум, возможно электро-наводки и т д
В результате получается график: Большая Низкочастотная синусоида, на которой "сидит" наш сигнал с другим высокочастотным шумом.
я предполагаю, что мне нужно:
1) отфильтровать низкие частоты и сгладить эту синусоиду до прямой
2)избавиться от высокочастотного шума
3) построить полученный чистый сигнал

подскажите что мне стоит использовать в в MATLAB в этом случае обычный Filter Design или wavelet'ы?
какие фильтры мне стоит? есть ли пример программ?
Ищу в интернете давно! Литературы полно, а вот примеров мало...
Сейчас нашел книгу Сергиенко "Цифровая обработка сигналов", но пока еще не нашел, того что мне нужно

[Alik]

Для начала я бы попробовал обычный ФНЧ первого порядка, тот самый который состоит из резистора и конденсатора. В цифровом виде это фильтр с конечной импульсной характеристикой (FIR). Тривиальная вещь, но вдруг поможет, ведь сигнал он по сути звук и если выводить его на динамик, это уже фильтр такого рода. Можно конечно сразу перейти к характеристикам помех, задействовать технику "вероятностных" фильтров типа Кальмана, но зачем? Вы ведь еще не знаете как подавляются помехи в простейшем случае, я имею ввиду не только коэффициент гармоник после ФНЧ, но и субъективную оценку качества сигнала на слух.

[denbur]

Имеется Фильтр низких частот Parks-McClellana.
с параметрами:
P1 = 34; % Период отсечки (в барах)
D1 = 23; % Период отсечки переходного процесса (в барах)
A1 = 40; % Затухание в полосе задержки, -Дб
rp = 0.08; % Биения в полосе пропускания, Дб

часть кода была взята отсюда:
http://www.mathworks.com/access/helpdes ... pmord.html


Помогите переделать его в полосовой фильтр Parks-McClellana с параметрами:
P1 = 10;
D1 = 8;
P2 = 40
D2 = 44
A1 = 40;
A2 = 40
rp = 0.08;


заранее спасибо


-------------------------------------------------------------------------------------------------
clear all;

x=[1.9985,1.9885,1.9871,1.9921,1.9929,1.9889,1.9933,1.9778,1.9812,1.979,1.9858,1.9771,1.975,1.9746,1.9705,1.9747,1.986,1.9846,1.9837,1.9827,1.9812,1.9758,1.9803,.9823];

fs = 1; % Частота дискретизации 1 бар в ед.времени
P1 = 34; % Период отсечки (в барах)
D1 = 23; % Период отсечки переходного процесса (в барах)
A1 = 40; % Затухание в полосе задержки, -Дб
rp = 0.08; % Биения в полосе пропускания, Дб
mydelay = 0; % Задержка (в барах)

f = [1/P1 1/D1]; % частоты, соответствующие периодам отсечки
a = [1 0]; % желаемые амплитуды в полосах пропускания (1) и задерживания (0)

dev = [(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-A1/20)];
[n,fo,ao,w] = firpmord(f,a,dev,fs); % вычисление минимального порядка фильтра
% порядок фильтра равен n+1
b = firpm(n,fo,ao) ; % вычисление симметричной импульсной характеристики

h = b((floor((n+1)/2)+1-mydelay):end); % вычисление половины импульсной характеристики (с учетом задержки)

h = h/sum(h); % нормализация урезанной импульсной характеристики (чтобы сумма её отсчетов была равна единице)


y = filter(h,1,x); % вычисление выхода фильтра
t = length(h):length(x);
plot(t,x(length(h):end),'b',t,y(length(h):end),'r'); % рисунок ряда цен и выхода фильтра
pause
freqz(h,1,1000,fs) % частотная характеристика фильтра
pause
impz(h,1); % импульсная характеристика фильтра


----------------------------------------------------------------------------