CIC фильтр по статьям

Все что касается фильтрации
kostar
Сообщения: 5
Зарегистрирован: 28 июн 2017, 19:47

CIC фильтр по статьям

Сообщение kostar »

Начитался прекрасных статей и решил сделать интерполяцию цифрового сигнала по
http://www.dsplib.ru/content/cic/cic.html
http://www.dsplib.ru/content/cicid/cicid.html

Делаю в лоб, по написанному.

Но после ФНЧ (интергаторов) получается полная абракадабра. Поскажите пожалуйста, ЧТЯДНТ?

Код: Выделить всё

close all; clc; clear;

Fs = 1000;
T = 1/Fs;
L = 1000;
t = (0:L-1)*T;
x = 5*sin(2*pi*50*t);

% Коэффициент интерполяции
R = 3;

% Величина задержки
D = 2*R;

% Порядок фильтра
N = 4;

b = 1;
a = 1;
for n = 1 : N
    b = conv(b, [1 0 -1]);
    a = conv(a, [1 -1]);
end
lb = length(b);
la = length(a);

% Усиление
K = 20*N*log10(D);

xR = zeros(1, floor(length(x)*R));
for n = 1 : length(x)
    xR(n*R) = x(max(1,n-lb+1):n) * b(end:-1:max(1,lb-n+1))';
end

xI = zeros(size(xR));
for k = 2 : length(xI)
    xI(k) = xR(k) - xI(max(1,k-la+1):k-1) * a(end-1:-1:max(1,la-k+1))';
end

kostar
Сообщения: 5
Зарегистрирован: 28 июн 2017, 19:47

Re: CIC фильтр по статьям

Сообщение kostar »

Проблема в коде была в том, что свертка начиналась с конца вектора коэффициентов, а не сначала. Поэтому и результат был неправильный.
Вместо

Код: Выделить всё

b(end:-1:max(1,lb-n+1))
a(end-1:-1:max(1,la-k+1))
в коде должно быть

Код: Выделить всё

b(min(lb, k):-1:1)
a(min(la, k):-1:2)
Тогда свертка с импульсной переходной функцией считается правильно.
Можно использовать альтернативное решение от Сергея Бахурина с добавлением нолей слева к исходному сигналу.

Теперь код программы выглядит так:

Код: Выделить всё

% Matlab 8.5.0.197613 (R2015a)
close all; clc; clear;

% Формирование сигнала.
Fs = 1000;
T = 1/Fs;
L = 50;
t = (0:L-1)*T;
x = 5*sin(2*pi*50*t);

% Коэффициент интерполяции
R = 3;

% Величина задержки
D = 4;

% Порядок фильтра
N = 4;

% Усиление
Klg = 20*N*log10(D);
K = D^N;

% Расчет коэффициентов передаточной функции
% H(w) = (1 - z^(-D))^N / (1 - z^-1)^N
b = 1;
a = 1;
for n = 1 : N
    b = linear_conv(b, [1 zeros(1,D-1) -1]);
    a = linear_conv(a, [1 -1]);
end
lb = length(b);
la = length(a);

% АЧХ и ФЧХ полученного фильтра стандартными средствами матлаб.
figure;
freqz(b,a);

% Гребенчатый КИХ фильтр
xC = zeros(size(x));
for n = 1 : length(x)
    xC(n) = x(max(1,n-lb+1):n) * b(min(lb, n):-1:1)';
end

% Добавляем нули для интерполяции
xR = zeros(1,length(xC)*R);
xR(1:R:end) = xC;

% Интергатор - БИХ фильтр
xI = zeros(size(xR));
for k = 1 : length(xI)
    xI(k) = xR(k) - xI(max(1,k-la+1):k-1) * a(min(la, k):-1:2)';
end

% Нормируем на коэффициент усиления
xI = xI / K;

% Применяем фильтр к сигналу стандартными средствами матлаб
y = filter(b,a,x);

% И смотрим результат
figure;
plot(t,x, t,y, (0:R*L-1)*T/R,xI)

xlabel('Время, с')
ylabel('Амлитуда сигнала, у.е.')
legend({'сигнал', 'filter()', 'интерполяция'})
И выдаёт следующий результат:
Изображение

Интерполированный сигнал задержан как относительно исходного сигнала, так и относительно выхода матлабовской фунции filter(). Почему это происходит, пока не понятно. На мой взгляд они должны быть синхронными. Плюс нормирование на коэффициент усиления из из статьи (формула 5)

даёт иную амплитуду сигнала.

Я считаю, что матлаб выдаёт правильный результат и стараюсь понять, где у меня пробелы в понимании происходящего.

Ещё интересные штуки с фильтром происходят при больших значениях и/или .

Аватара пользователя
skameykin22
Сообщения: 1
Зарегистрирован: 09 сен 2017, 08:18
Контактная информация:

Re: CIC фильтр по статьям

Сообщение skameykin22 »

Полезная информация.

Ответить