В этом разделе представлен пример того, как можно создать фильтр с целочисленными коэффициентами. В этом примере создается фильтр приподнятого косинуса с коэффициентами с плавающей точкой, и коэффициенты фильтра затем преобразуются в целые числа.
Чтобы проиллюстрировать концепции использования целых чисел с фильтрами с фиксированной точкой, в этом примере будет использоваться фильтр приподнятого косинуса:
b = rcosdesign(.25, 12.5, 8, 'sqrt');
b нормированы так, что коэффициент усиления полосы пропускания равен 1 и все меньше 1. В порядок, чтобы сделать их целыми числами, они должны быть масштабированы. Если вы хотели масштабировать их, чтобы использовать 18 биты для каждого коэффициента, область значений возможных значений для коэффициентов становится:
Потому что самый большой коэффициент b положительно, его нужно будет масштабировать как можно ближе к 131071 (без переполнения), порядок минимизировать ошибку квантования. Можно определить экспоненту фактора шкалы путем выполнения:
B = 18; % Number of bits L = floor(log2((2^(B-1)-1)/max(b))); % Round towards zero to avoid overflow bsc = b*2^L;
Кроме того, можно использовать инструмент автомасштабирования чисел с фиксированной точкой следующим образом:
bq = fi(b, true, B); % signed = true, B = 18 bits L = bq.FractionLength;
Это совпадение, что B и L в этом случае их обоих 18, из-за значения наибольшего коэффициента b. Если, например, максимальное значение b было 0,124, L будет 20, в то время как B (количество бит) останется 18.
Сначала создайте фильтр с помощью структуры линии задержки с прямой формой:
h = dfilt.dffir(bsc);
В порядок установки необходимых параметров арифметика должна быть установлена на фиксированную точку:
h.Arithmetic = 'fixed'; h.CoeffWordLength = 18;
Можно проверить, что коэффициенты h все целые числа:
all(h.Numerator == round(h.Numerator)) ans = 1
Теперь можно изучить величину реакцию фильтра, используя fvtool:
fvtool(h, 'Color', 'white')
Это показывает большой коэффициент усиления 117 дБ в полосе пропускания, что связано с большими значениями коэффициентов - это приведет к тому, что выход фильтра будет намного больше, чем вход. Метод решения этой проблемы будет рассмотрен в следующих разделах.
Вам нужно будет задать входные параметры вашего фильтра в соответствующие значения для работы с целыми числами. Для примера, если вход в фильтр от A/D конвертера с 12-битным разрешением, вы должны установить вход следующим образом:
h.InputWordLength = 12; h.InputFracLength = 0;
The info метод возвращает сводные данные настроек фильтра.
info(h)
Discrete-Time FIR Filter (real) ------------------------------- Filter Structure : Direct-Form FIR Filter Length : 101 Stable : Yes Linear Phase : Yes (Type 1) Arithmetic : fixed Numerator : s18,0 -> [-131072 131072) Input : s12,0 -> [-2048 2048) Filter Internals : Full Precision Output : s31,0 -> [-1073741824 1073741824) (auto determined) Product : s29,0 -> [-268435456 268435456) (auto determined) Accumulator : s31,0 -> [-1073741824 1073741824) (auto determined) Round Mode : No rounding Overflow Mode : No overflow
В этом случае все дробные длины теперь равны нулю, что означает, что фильтр h настроен на обработку целых чисел.
Можно сгенерировать входной сигнал для фильтра путем квантования до 12 бит с помощью функции автоматического сканирования, или можно следовать той же процедуре, которая использовалась для коэффициентов, обсуждавшихся ранее. В этом примере создайте сигнал с двумя синусоидами:
n = 0:999; f1 = 0.1*pi; % Normalized frequency of first sinusoid f2 = 0.8*pi; % Normalized frequency of second sinusoid x = 0.9*sin(0.1*pi*n) + 0.9*sin(0.8*pi*n); xq = fi(x, true, 12); % signed = true, B = 12 xsc = fi(xq.int, true, 12, 0);
Для фильтрации входного сигнала, сгенерированного выше, введите следующее:
ysc = filter(h, xsc);
Вот ysc является полным выходом точности, что означает, что в расчете не было сброшено ни одних бит. Это делает ysc наилучший выход, которого вы можете достичь, учитывая 12-битный вход и 18-битные коэффициенты. Это может быть проверено путем фильтрации с помощью двойной точности с плавающей точностью и сравнения результатов двух операций фильтрации:
hd = double(h);
xd = double(xsc);
yd = filter(hd, xd);
norm(yd-double(ysc))
ans =
0
Теперь можно изучить выход по сравнению с входом. Этот пример строит график только последних нескольких выборок, чтобы минимизировать эффект переходных процессов:
idx = 800:950;
xscext = double(xsc(idx)');
gd = grpdelay(h, [f1 f2]);
yidx = idx + gd(1);
yscext = double(ysc(yidx)');
stem(n(idx)', [xscext, yscext]);
axis([800 950 -2.5e8 2.5e8]);
legend('input', 'output');
set(gcf, 'color', 'white');
Трудно сравнить два сигнала на этом рисунке из-за большого различия в шкалах. Это связано с большим усилением фильтра, поэтому вам нужно будет компенсировать усиление фильтра:
stem(n(idx)', [2^18*xscext, yscext]);
axis([800 950 -5e8 5e8]);
legend('scaled input', 'output');
Вы можете увидеть, как сигналы сравниваются намного легче, когда масштабирование было сделано, как видно на вышеприведенном рисунке.
Если вы исследуете выход словосочетание,
ysc.WordLength
ans =
31
вы заметите, что количество бит в выходе значительно больше, чем во входе. Поскольку такой рост количества бит, представляющих данные, может быть нежелательным, вам может потребоваться обрезать словесную длину выхода. Лучший способ сделать это - отбросить наименее значимые биты, порядок минимизировать ошибку. Однако, если вы знаете, что есть неиспользованные биты высокого порядка, вы должны также отбросить эти биты.
Чтобы определить, существуют ли неиспользованные наиболее значимые биты (MSBs), можно посмотреть, где рост WordLength возникает в расчетах. В этом случае рост разрядности происходит, чтобы учесть результаты сложения продуктов входных данных (12 биты) и коэффициентов (18 биты). Каждый из этих продуктов имеет длину 29 бит (можно проверить это с помощью info(h)). Рост разрядности из-за накопления продукта зависит от длины фильтра и значений коэффициента - однако это определение в худшем случае в том смысле, что никакое предположение по входному сигналу не сделано, кроме того, и в результате могут быть неиспользованные MSB. Вы должны будете быть осторожны, так как MSB, которые считаются неиспользованными неправильно, вызовут переполнение.
Предположим, вы хотите сохранить 16 биты для выхода. В этом случае нет роста разрядности из-за сложений, поэтому настройка выхода бита будет 16 для wordlength и -14 для длины дроби.
Поскольку фильтрация уже выполнена, можно отбросить некоторые биты от ysc:
yout = fi(ysc, true, 16, -14);
Также можно задать выходные битовые длины фильтра непосредственно (это полезно, если вы планируете фильтровать много сигналов):
specifyall(h); h.OutputWordLength = 16; h.OutputFracLength = -14; yout2 = filter(h, xsc);
Можно проверить, что результаты совпадают в любом случае:
norm(double(yout) - double(yout2))
ans =
0
Однако, если сравнить это с полным выходом точности, вы заметите, что существует ошибка округления из-за отброшенных бит:
norm(double(yout)-double(ysc))
ans =
1.446323386867543e+005
В этом случае различия трудно заметить при построении графиков данных, как показано ниже:
stem(n(yidx), [double(yout(yidx)'), double(ysc(yidx)')]);
axis([850 950 -2.5e8 2.5e8]);
legend('Scaled Input', 'Output');
set(gcf, 'color', 'white');
Поскольку фильтр в этом примере имеет такой большой коэффициент усиления, выход находится в иной шкале, чем вход. Однако это масштабирование является чисто теоретическим, и вы можете масштабировать данные, как вам нравится. В этом случае у вас есть 16 биты для выхода, но вы можете прикрепить любое масштабирование, которое вы выберете. Было бы естественно переосмыслить выход, чтобы иметь вес 2 ^ 0 (или L = 0) для LSB. Это эквивалентно масштабированию выходного сигнала вниз на множитель 2 ^ (-14). Однако ошибка расчета или округления отсутствует. Это можно сделать, выполнив следующее:
yri = fi(yout.int, true, 16, 0);
stem(n(idx)', [xscext, double(yri(yidx)')]);
axis([800 950 -1.5e4 1.5e4]);
legend('input', 'rescaled output');
Этот график показывает, что выход все еще больше, чем вход. Если бы вы сделали фильтрацию в двойной точности с плавающей точностью, это было бы не так - потому что здесь для выхода используется больше бит, чем для входа, поэтому MSB взвешиваются по-другому. Другой способ можно увидеть, рассмотрев величину ответ масштабированного фильтра:
[H,w] = freqz(h); plot(w/pi, 20*log10(2^(-14)*abs(H)));
Этот график показывает, что коэффициент усиления полосы пропускания все еще выше 0 дБ.
Чтобы поместить вход и выход в одну шкалу, MSB должны быть взвешены одинаково. MSB входа имеет вес 2 ^ 11, в то время как MSB с масштабированным выходом имеет вес 2 ^ (29-14) = 2 ^ 15. Вы должны дать выходу MSB вес 2 ^ 11 следующим образом:
yf = fi(zeros(size(yri)), true, 16, 4);
yf.bin = yri.bin;
stem(n(idx)', [xscext, double(yf(yidx)')]);
legend('input', 'rescaled output');
Эта операция эквивалентна масштабированию коэффициента усиления фильтра на 2 ^ (-18).
[H,w] = freqz(h); plot(w/pi, 20*log10(2^(-18)*abs(H)));
Вышеописанный график показывает усиление 0 дБ в полосе пропускания, при желании.
С помощью этой окончательной версии выхода, yf больше не является целым числом. Однако это связано только с интерпретацией - целыми числами, представленными битами в yf идентичны тем, которые представлены битами в yri. Можно проверить это, сравнив их:
max(abs(yf.int - yri.int))
ans =
0
Установите параметры фильтра, чтобы работать с целыми числами
The set2int способ обеспечивает удобный способ установки параметров фильтра для работы с целыми числами. Метод работает, масштабируя коэффициенты до целочисленных чисел и устанавливая коэффициенты и входную длину дроби в нуль. Это позволяет вам использовать коэффициенты с плавающей точкой непосредственно.
h = dfilt.dffir(b); h.Arithmetic = 'fixed';
Коэффициенты представлены 18 битами, а входной сигнал представлен 12 битами:
g = set2int(h, 18, 12);
g_dB = 20*log10(g)
g_dB =
1.083707984390332e+002
The set2int метод возвращает коэффициент усиления фильтра путем масштабирования коэффициентов до целых чисел, поэтому коэффициент усиления всегда является степенью 2. Вы можете проверить, что усиление, которое мы получаем здесь, согласуется с усилением фильтра ранее. Теперь можно также проверить, что фильтр h правильно настроен для работы с целыми числами:
info(h) Discrete-Time FIR Filter (real) ------------------------------- Filter Structure : Direct-Form FIR Filter Length : 101 Stable : Yes Linear Phase : Yes (Type 1) Arithmetic : fixed Numerator : s18,0 -> [-131072 131072) Input : s12,0 -> [-2048 2048) Filter Internals : Full Precision Output : s31,0 -> [-1073741824 1073741824) (auto determined) Product : s29,0 -> [-268435456 268435456) (auto determined) Accumulator: s31,0 -> [-1073741824 1073741824) (auto determined) Round Mode : No rounding Overflow Mode : No overflow
Здесь можно увидеть, что все дробные длины теперь равны нулю, поэтому этот фильтр настроен правильно для работы с целыми числами.
Повторная интерпретация выхода
Можно сравнить выход с выходным сигналом ссылки с двойной точностью с плавающей точностью и проверить, что расчеты выполняются фильтром h делается в полной точности.
yint = filter(h, xsc);
norm(yd - double(yint))
ans =
0
Затем можно обрезать выход только до 16 бит:
yout = fi(yint, true, 16);
stem(n(yidx), [xscext, double(yout(yidx)')]);
axis([850 950 -2.5e8 2.5e8]);
legend('input', 'output');
Еще раз, график показывает, что вход и выход находятся в разных шкалах. В порядок масштабирования выхода так, чтобы сигналы можно было легче сравнивать на графике, вам нужно будет взвесить MSB соответствующим образом. Можно вычислить новую длину дроби с помощью усиления фильтра, когда коэффициенты были целыми числами:
WL = yout.WordLength;
FL = yout.FractionLength + log2(g);
yf2 = fi(zeros(size(yout)), true, WL, FL);
yf2.bin = yout.bin;
stem(n(idx)', [xscext, double(yf2(yidx)')]);
axis([800 950 -2e3 2e3]);
legend('input', 'rescaled output');
Этот окончательный график показывает фильтрованные данные, повторно масштабированные, чтобы соответствовать входу шкалы.