Создайте КИХ-фильтр Используя целочисленные коэффициенты
Этот раздел обеспечивает пример того, как можно создать фильтр с целочисленными коэффициентами. В этом примере создается фильтр приподнятого косинуса с коэффициентами с плавающей точкой, и коэффициенты фильтра затем преобразованы в целые числа.
Задайте коэффициенты фильтра
Чтобы проиллюстрировать концепции использования целых чисел с фильтрами фиксированной точки, этот пример будет использовать фильтр приподнятого косинуса:
b = rcosdesign(.25, 12.5, 8, 'sqrt');
b нормированы так, чтобы усиление полосы пропускания было равным 1 и было все меньшим, чем 1. Для того, чтобы сделать их целыми числами, они должны будут масштабироваться. Если вы хотели масштабировать их, чтобы использовать 18 битов для каждого коэффициента, область значений возможных значений для коэффициентов становится:
Поскольку самый большой коэффициент b положительно, это должно будет масштабироваться максимально близко к 131 071 (не переполняясь) для того, чтобы минимизировать ошибку квантования. Можно определить экспоненту масштабного коэффициента путем выполнения:
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;
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');
Вы видите, как сигналы выдерживают сравнение намного более легко, если масштабирование было сделано, как замечено на вышеупомянутом рисунке.
Обрежьте Выходной WordLength
Если вы исследуете выход wordlength,
ysc.WordLength
ans =
31
вы заметите, что количество битов в выходе значительно больше во входе. Поскольку такой рост в количестве битов, представляющих данные, не может быть желательным, вы, возможно, должны обрезать wordlength выхода. Лучший способ сделать это должно отбросить младшие значащие биты, для того, чтобы минимизировать ошибку. Однако, если вы знаете, что существуют неиспользованные высокого уровня биты, необходимо отбросить те биты также.
Чтобы определить, существуют ли неиспользованные старшие значащие биты (MSBs), можно посмотреть на то, где рост в WordLength возникает в расчете. В этом случае рост разрядности происходит, чтобы вместить результаты добавляющих продуктов входа (12 битов) и коэффициентов (18 битов). Каждый из этих продуктов 29 битов длиной (можно проверить это использование info(h)). Рост разрядности из-за накопления продукта зависит от длины фильтра и содействующих значений - однако, это - определение худшего случая в том смысле, что никакое предположение на входном сигнале не сделано кроме того, и в результате может быть неиспользованный MSBs. Необходимо будет быть осторожными, хотя, когда MSBs, которые считают неиспользованными неправильно, вызовут переполнение.
Предположим, что вы хотите сохранить 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');
Этот график показывает, что выход еще больше, чем вход. Если бы вы сделали просачивание с двойной точностью с плавающей точкой, это не имело бы место — потому что здесь больше битов используется для выхода, чем для входа, таким образом, MSBs взвешиваются по-другому. Вы видите это иначе путем рассмотрения ответа величины масштабированного фильтра:
[H,w] = freqz(h); plot(w/pi, 20*log10(2^(-14)*abs(H)));
Этот график показывает, что усиление полосы пропускания все еще выше 0 дБ.
Чтобы поместить ввод и вывод на ту же шкалу, MSBs должен быть взвешен одинаково. Вход 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
Сконфигурируйте Параметры Фильтра, чтобы работать с Целыми числами Используя set2int Метод
Усадите параметры фильтра за работу с Целыми числами
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
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');
Еще раз график показывает, что ввод и вывод в различных шкалах. Для того, чтобы масштабировать выход так, чтобы сигналы могли быть сравнены более легко в графике, необходимо будет взвесить MSBs соответственно. Можно вычислить новую дробную длину с помощью усиления фильтра, когда коэффициенты были целыми числами:
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');
Этот итоговый график показывает отфильтрованные данные, перемасштабированные, чтобы совпадать с входной шкалой.