Страница 1 из 1

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

Добавлено: 06 июл 2017, 17:47
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

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

Добавлено: 13 июл 2017, 02:55
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)

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

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

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

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

Добавлено: 09 сен 2017, 18:07
skameykin22
Полезная информация.