從零開始在FPGA上實(shí)現(xiàn)IIR濾波器

描述

濾波可能是任何嵌入式工程師必須設(shè)計(jì)的最常見的 DSP 算法，無(wú)論他們是為 STM32、TI DSP 還是為 FPGA 開發(fā)的算法。濾波非常重要，因?yàn)榇蠖鄶?shù)應(yīng)用程序都連接到現(xiàn)實(shí)世界，因此它們需要捕獲具有所有真實(shí)信號(hào)都具有的不良元素的真實(shí)信號(hào)，即噪聲、偏移……此外，我們需要獲取的信號(hào)，并非總是如此將是我們獲取的主要信號(hào)，例如，對(duì)于生物信號(hào)，我們將獲取與電網(wǎng)頻率相對(duì)應(yīng)的 50 或 60 Hz 的高電平。在這種情況下，我們將需要一個(gè)過(guò)濾器。

作為本博客的讀者，您肯定知道，在數(shù)字系統(tǒng)上，我們有 2 種濾波器，F(xiàn)IR（有限脈沖響應(yīng)）濾波器，它只使用現(xiàn)在輸入的值和過(guò)去的值，以及 IIR（無(wú)限脈沖Response），它使用現(xiàn)在和過(guò)去的輸入值，以及過(guò)去輸出的值。FIR 濾波器很容易為許多應(yīng)用設(shè)計(jì)，但它們的響應(yīng)對(duì)于低階濾波器是有限的。另一方面，使用 IIR 濾波器，我們可以實(shí)現(xiàn)更積極的響應(yīng)，但它們的缺點(diǎn)是濾波器在某些情況下可能不穩(wěn)定。在本文中，我們將了解如何逐步實(shí)現(xiàn)二階濾波器，從濾波器設(shè)計(jì)到 Verilog 實(shí)現(xiàn)。

設(shè)計(jì) IIR 濾波器的過(guò)程始終相同。首先我們要使用它的連續(xù)傳遞函數(shù)來(lái)設(shè)計(jì)濾波器，然后，一旦選擇了濾波器的固有頻率、階數(shù)和品質(zhì)因數(shù)，我們將得到相應(yīng)的 s

功能。為了將該濾波器數(shù)字化，下一步是應(yīng)用一些變換來(lái)用 z 替換變量s ，這個(gè)過(guò)程稱為離散化。為了離散化連續(xù)傳遞函數(shù)，我們可以使用多種方法，其中一種易于使用的方法是Tustin或雙線性變換。

雙線性變換基于將連續(xù)傳遞函數(shù)中的s變量替換為z的函數(shù)。在下一個(gè)示例中，我們可以看到應(yīng)用雙線性變換對(duì)單極點(diǎn)低通濾波器進(jìn)行離散化。接下來(lái)是我們必須在連續(xù)傳遞函數(shù)中執(zhí)行的替換。

?

對(duì)于單極點(diǎn)濾波器，這些將是執(zhí)行雙線性變換的步驟。低通單極濾波器的方程如下：

?

應(yīng)用轉(zhuǎn)換，我們有：

?

與分母運(yùn)算以獲得共同的分母，

?

然后，我們需要將元素與其他元素中的z分開，最后我們將得到一個(gè)z函數(shù)，該函數(shù)取決于濾波器的采樣頻率和固有頻率。

?

在 MATLAB 中，我們可以使用以下命令進(jìn)行此轉(zhuǎn)換hz = c2d(hs,ts,'Tustin')

接下來(lái)是 100 Hz 低通濾波器的整個(gè) MATLAB 代碼。

%% single pole filter discretization

% continuos filter declaration
s = tf('s');

% define characteristics of the filter
fc = 100; 

% compute angular frequency
wc = 2*pi*fc;

% Write the transfer function of a low pass single pole filter
hs = 1/(1+s/wc);

% show its bode diagram
bode(hs)

% discretization using tustin method

% Sampling frequency
fs = 100e3;

ts = 1/fs;

hz = c2d(hs,ts,'Tustin')

赫茲的結(jié)果是

hz =
 
  0.003132 z + 0.003132
  ---------------------
       z - 0.9937
 
Sample time: 1e-05 seconds
Discrete-time transfer function.

為了驗(yàn)證我們手動(dòng)獲得的方程是否正確，我們必須使用我們之前獲得的值創(chuàng)建一個(gè)傳遞函數(shù)。

hz_man = tf([ts*wc, ts*wc],[ts*wc+2 ts*wc-2], ts)

這種情況下的結(jié)果是

hz_man =
 
  0.006283 z + 0.006283
  ---------------------
     2.006 z - 1.994
 
Sample time: 1e-05 seconds
Discrete-time transfer function.

我們可以看到，將分子和分母除以 2.006，方程與使用c2c命令得到的方程相同。

雙線性變換與其他任何離散化方法一樣，都是一種近似，因此我們將引入一些誤差。該誤差的大小與采樣頻率有關(guān)，因?yàn)樵谀承┣闆r下，非常（非常）高的采樣頻率可能看起來(lái)是連續(xù)的。在實(shí)際系統(tǒng)中，采樣頻率取決于 ADC 和許多實(shí)際組件，我們可以使用高采樣頻率，但在其他情況下，我們會(huì)引入顯著誤差。

在雙線性變換的情況下，我們將引入的誤差之一是頻率扭曲。這種效應(yīng)會(huì)產(chǎn)生濾波器固有頻率的偏移。這種影響在低通濾波器的情況下可以忽略不計(jì)，但在窄帶濾波器的情況下，該誤差可能使濾波器不起作用。

為了找到由雙線性變換引起的誤差，我們需要找到模擬頻率與離散頻率之間的關(guān)系。首先，我們將從雙線性變換方程開始。

?

將z替換為其極坐標(biāo)形式e^(jωd) ，我們有

?

現(xiàn)在，我們將使用半角分割指數(shù)。

?

將分子和分母替換為上述等式，我們有：

?

現(xiàn)在，我們將使用半角分割指數(shù)。

?

將分子和分母替換為上述等式，我們有：

?

然后，我們必須找到與類比的關(guān)系，所以我們要使兩個(gè)方程相等

?

我們可以驗(yàn)證實(shí)部為 0，而虛部wwa直接為

?

請(qǐng)注意，當(dāng)ts趨于零時(shí)，極限是ωa = ωd。

在 MATLAB 中，我們可以很容易地檢查這個(gè)偏差。在這種情況下，我們將聲明一個(gè)二階連續(xù)帶阻濾波器。然后使用 Tustin 方法對(duì)濾波器進(jìn)行離散化。

close all
clear all
clc

% define filter
fc = 3000;
wn = 2*pi*fc;
q = 5;%1/sqrt(2); % butterworth

s = tf('s');

h = (s^2+wn^2)/(s^2+wn/q*s+wn^2)

%bode(h)

% discretize filter
fs = 100e6/2048 % prescale FPGA frequency by factor of 2048
ts = 1/fs;

hz = c2d(h,ts,'tustin');

?

我們可以看到數(shù)字濾波器的頻率是如何移動(dòng)的。然后我們將對(duì)自然頻率應(yīng)用預(yù)變形，并使用這個(gè)新頻率重新離散化濾波器。

% calculate prewarping
wp = 2/ts*tan(wc*ts/2);
hw = (s^2+wp^2)/(s^2+wp/q*s+wp^2)
hwz = c2d(hw,ts,'tustin');

?

我們可以看到在這種情況下，兩個(gè)頻率是相同的，所以誤差是固定的。

為了在現(xiàn)實(shí)世界中驗(yàn)證這種行為，我們將在 Verilog 中實(shí)現(xiàn)這個(gè)二階帶阻濾波器，然后我們將檢查扭曲錯(cuò)誤，然后如何使用預(yù)扭曲來(lái)修復(fù)它。

在開始編碼之前，我們必須考慮我們想要?jiǎng)?chuàng)建的模塊是什么。在我的情況下，能夠參數(shù)化模塊是非常重要的，這樣，如果我改變輸入的寬度，或者輸出的寬度，那只代表一個(gè)參數(shù)的改變，而不需要修改模塊的代碼因?yàn)檫@意味著模塊的新測(cè)試。此外，我們必須考慮我們將使用什么數(shù)字格式。對(duì)于這種過(guò)濾器，很明顯我們需要有符號(hào)格式，以及管理十進(jìn)制數(shù)的能力. 為了達(dá)到這個(gè)要求，我們必須在定點(diǎn)和浮點(diǎn)之間進(jìn)行選擇，對(duì)我來(lái)說(shuō)決定很明確，我想要一個(gè)可以單獨(dú)實(shí)例化的代碼，沒(méi)有外部浮點(diǎn)單元，所以格式將是定點(diǎn)的。如果我們加上參數(shù)的要求，和定點(diǎn)格式，我們需要參數(shù)來(lái)定義信號(hào)的寬度，還有小數(shù)部分的寬度，這與我們需要的精度有關(guān)。此外，將選擇精度以獲得與設(shè)計(jì)盡可能相似的實(shí)施濾波器的響應(yīng)。現(xiàn)在，我們需要參數(shù)化多少寬度？我們可以定義一個(gè)寬度，用于輸入和輸出以及系數(shù)，但是這樣，系數(shù)的寬度將根據(jù)數(shù)據(jù)生成器的寬度來(lái)選擇，對(duì)于某些過(guò)濾器，系數(shù)在穩(wěn)定極限上，這可能代表一個(gè)問(wèn)題。因此，至少我們將定義 2 種不同的寬度，一種用于輸入和輸出，根據(jù)數(shù)據(jù)消耗和數(shù)據(jù)源，另一種根據(jù)我們需要的分辨率。另一個(gè)我們要考慮寬度的事情，是內(nèi)部運(yùn)算的分辨率，因?yàn)樵谀承┣闆r下，穩(wěn)定性問(wèn)題與系數(shù)本身的寬度無(wú)關(guān)，而是運(yùn)算分辨率，所以在這一點(diǎn)上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數(shù)的寬度，以及內(nèi)部濾波器操作的寬度。最后，過(guò)濾器參數(shù)將如下所示。一個(gè)用于輸入和輸出，根據(jù)數(shù)據(jù)消耗和數(shù)據(jù)源，另一個(gè)根據(jù)我們需要的分辨率。另一個(gè)我們要考慮寬度的事情，是內(nèi)部運(yùn)算的分辨率，因?yàn)樵谀承┣闆r下，穩(wěn)定性問(wèn)題與系數(shù)本身的寬度無(wú)關(guān)，而是運(yùn)算分辨率，所以在這一點(diǎn)上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數(shù)的寬度，以及內(nèi)部濾波器操作的寬度。最后，過(guò)濾器參數(shù)將如下所示。一個(gè)用于輸入和輸出，根據(jù)數(shù)據(jù)消耗和數(shù)據(jù)源，另一個(gè)根據(jù)我們需要的分辨率。另一個(gè)我們要考慮寬度的事情，是內(nèi)部運(yùn)算的分辨率，因?yàn)樵谀承┣闆r下，穩(wěn)定性問(wèn)題與系數(shù)本身的寬度無(wú)關(guān)，而是運(yùn)算分辨率，所以在這一點(diǎn)上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數(shù)的寬度，以及內(nèi)部濾波器操作的寬度。最后，過(guò)濾器參數(shù)將如下所示。穩(wěn)定性問(wèn)題與系數(shù)本身的寬度無(wú)關(guān)，而與運(yùn)算分辨率有關(guān)，因此在這一點(diǎn)上，將在 MATLAB 或 Python 上生成的系數(shù)的寬度與內(nèi)部濾波器操作的寬度解耦會(huì)很有趣. 最后，過(guò)濾器參數(shù)將如下所示。穩(wěn)定性問(wèn)題與系數(shù)本身的寬度無(wú)關(guān)，而與運(yùn)算分辨率有關(guān)，因此在這一點(diǎn)上，將在 MATLAB 或 Python 上生成的系數(shù)的寬度與內(nèi)部濾波器操作的寬度解耦會(huì)很有趣. 最后，過(guò)濾器參數(shù)將如下所示。

parameter inout_width = 16,
parameter inout_decimal_width = 15,
parameter coefficient_width = 16,
parameter coefficient_decimal_width = 15,
parameter internal_width = 16,
parameter internal_decimal_width = 15

接下來(lái)我們必須考慮接口。在模塊之間以連續(xù)方式傳輸數(shù)據(jù)的應(yīng)用中，AXI4-Stream 將是最佳選擇。在輸入和輸出字段中，我們將定義主從 AXI4-Stream 接口來(lái)獲取和發(fā)送數(shù)據(jù)。盡管輸入和輸出寬度是參數(shù)化的，但如果我們想將模塊連接到現(xiàn)有的 AXI4-Stream IP，則此寬度受限于總線的寬度。

input aclk,
input resetn,

/* slave axis interface */
input [inout_width-1:0] s_axis_tdata,
input s_axis_tlast,
input s_axis_tvalid,
output s_axis_tready,

/* master axis interface */
output reg [inout_width-1:0] m_axis_tdata,
output reg m_axis_tlast,
output reg m_axis_tvalid,
input m_axis_tready,

關(guān)于系數(shù)，插入它們最好是使用 AXI4-Lite 接口，但我想設(shè)計(jì)這個(gè)模塊而不需要使用 Zynq 或 Microblaze，所以插入系數(shù)的簡(jiǎn)單方法是作為輸入。

/* coefficients */
input signed [pw_coefficient_width-1:0] b0,
input signed [pw_coefficient_width-1:0] b1,
input signed [pw_coefficient_width-1:0] b2,
input signed [pw_coefficient_width-1:0] a1,
input signed [pw_coefficient_width-1:0] a2

現(xiàn)在，在定義了所有輸入和輸出之后，我們必須考慮數(shù)據(jù)流。首先，由于我們定義了不同的格式，為了能夠在系數(shù)和數(shù)據(jù)之間進(jìn)行運(yùn)算，我們必須將所有信號(hào)的格式更改為內(nèi)部格式，內(nèi)部格式由內(nèi)部寬度和小數(shù)寬度定義。使用的整數(shù)寬度被定義為一個(gè)localparam，并且對(duì)應(yīng)于寬度和小數(shù)寬度之間的差異。要更改格式，我們將在信號(hào)的低端填充零，直到完成小數(shù)部分。對(duì)于整數(shù)部分，由于使用的格式是有符號(hào)的，我們必須進(jìn)行符號(hào)擴(kuò)展，即填充 MSb 的值，直到整數(shù)部分完成。請(qǐng)注意，這是有效的，因?yàn)閮?nèi)部寬度大于或等于輸入輸出和系數(shù)寬度。

/* resize signals to internal width */
assign input_int = { {(internal_integer_width-inout_integer_width){s_axis_tdata[inout_width-1]}},
                          s_axis_tdata,
                          {(internal_decimal_width-inout_decimal_width){1'b0}} };
assign b0_int = { {(internal_integer_width-coefficient_integer_width){b0[coefficient_width-1]}},
                          b0,
                          {(internal_decimal_width-coefficient_decimal_width){1'b0}} };
assign b1_int = { {(internal_integer_width-coefficient_integer_width){b1[coefficient_width-1]}},
                          b1,
                          {(internal_decimal_width-coefficient_decimal_width){1'b0}} };
assign b2_int = { {(internal_integer_width-coefficient_integer_width){b2[coefficient_width-1]}},
                          b2,
                          {(internal_decimal_width-coefficient_decimal_width){1'b0}} };
assign a1_int = { {(internal_integer_width-coefficient_integer_width){a1[coefficient_width-1]}},
                          a1,
                          {(internal_decimal_width-coefficient_decimal_width){1'b0}} };
assign a2_int = { {(internal_integer_width-coefficient_integer_width){a2[coefficient_width-1]}},
                          a2,
                          {(internal_decimal_width-coefficient_decimal_width){1'b0}} };

數(shù)字濾波器，無(wú)論是 FIR 還是 IIR，都需要存儲(chǔ)過(guò)去輸入的值，而 IIR 還需要存儲(chǔ)過(guò)去輸出的值，因此我們需要一個(gè)流水線結(jié)構(gòu)來(lái)存儲(chǔ)過(guò)去的值。

/* pipeline registers */
always @(posedge aclk)
  if (!resetn) begin
    input_pipe1 <= 0;
    input_pipe2 <= 0;
    output_pipe1 <= 0;
    output_pipe2 <= 0;
  end
  else
    if (s_axis_tvalid) begin
      input_pipe1 <= input_int;
      input_pipe2 <= input_pipe1;
      output_pipe1 <= output_int;
      output_pipe2 <= output_pipe1;
    end

現(xiàn)在，接下來(lái)是執(zhí)行過(guò)濾器計(jì)算。一個(gè)二階濾波器需要執(zhí)行 5 次乘法，記住這并不意味著模塊使用 5 個(gè) DSP slice。我開發(fā)的代碼執(zhí)行組合乘法。這允許 0 個(gè)時(shí)鐘周期延遲，但會(huì)限制濾波器的速度。如果您的時(shí)序約束沒(méi)有得到滿足，您可以注冊(cè)乘法器的輸入，然后進(jìn)行重新定時(shí)，讓 Vivado 選擇更有效的放置寄存器的位置。

/* combinational multiplications */
assign input_b0 = input_int * b0_int;
assign input_b1 = input_pipe1 * b1_int;
assign input_b2 = input_pipe2 * b2_int;
assign output_a1 = output_pipe1 * a1_int;
assign output_a2 = output_pipe2 * a2_int;

乘法的結(jié)果將在輸入的情況下相加，在輸出的情況下相減，以獲得濾波器輸出。由于乘積的輸出大小是操作數(shù)的兩倍，因此加法必須是相同的寬度。要在乘法上使用輸出，我們必須執(zhí)行移位以刪除額外的小數(shù)位。關(guān)于額外的整數(shù)位置將在分配時(shí)被截?cái)唷?/font>

assign output_2int = input_b0 + input_b1 + input_b2 - output_a1 - output_a2;
assign output_int = output_2int >>> (internal_decimal_width);

最后，輸出的值將被重新格式化為輸入輸出寬度。

assign m_axis_tdata = output_int >>> (internal_decimal_width-inout_decimal_width);

對(duì)于 AXI4-Stream 管理信號(hào)，由于濾波器作為橋接器，有 1 個(gè)周期延遲，因此管理信號(hào)會(huì)以 1 個(gè)周期延遲通過(guò)濾波器。

至此，我們已經(jīng)實(shí)現(xiàn)了 Verilog 模塊并準(zhǔn)備好用作帶阻濾波器。接下來(lái)是使用該tfdata函數(shù)在 MATLAB 中獲取系數(shù)。關(guān)于過(guò)濾器的寬度，我們將使用 18 位分辨率用于小數(shù)部分，2 用于整數(shù)部分，以使過(guò)濾器達(dá)到 +-2。接下來(lái)是獲取定點(diǎn)系數(shù)的 MATLAB 代碼。

% filter implemetation

[num, den] = tfdata(hwz)

b0 = num{1}(1)
b1 = num{1}(2)
b2 = num{1}(3)
a1 = den{1}(2)
a2 = den{1}(3)

% quantize

fracBits = 18;

b0_qq = floor(b0*2^fracBits)
b1_qq = floor(b1*2^fracBits)
b2_qq = floor(b2*2^fracBits)

a1_qq = floor(a1*2^fracBits)
a2_qq = floor(a2*2^fracBits)

這將返回下一個(gè)結(jié)果。

b0_qq =

      252631


b1_qq =

     -468081


b2_qq =

      252631


a1_qq =

     -468081


a2_qq =

      243119

一旦我們有了五個(gè)系數(shù)，我們必須將它們作為輸入引入 Verilog 模塊中的系數(shù)輸入。

axis_biquad_v1_0 #(
  .inout_width(15),
  .inout_decimal_width(13),
  .coefficient_width(20),
  .coefficient_decimal_width(18),
  .internal_width(20),
  .internal_decimal_width(18)
  ) axis_biquad_inst0 (
  .aclk(clk100mhz),
  .resetn(pll_locked),
  /* slave axis interface */
  .s_axis_tdata({axis_data_signal_ch1[13], axis_data_signal_ch1}),
  .s_axis_tlast(),
  .s_axis_tvalid(&filter_prescaler),
  .s_axis_tready(),
  /* master axis interface */
  .m_axis_tdata(axis_data_signal_filt),
  .m_axis_tlast(),
  .m_axis_tvalid(),
  .m_axis_tready(),
  /* coefficients */
  .b0(20'd252631),
  .b1(-20'd468081),
  .b2(20'd252631),
  .a1(-20'd468081),
  .a2(20'd243119)
  );

為了驗(yàn)證這種效果是否發(fā)生，并使用預(yù)變形修復(fù)，我將使用帶有 ZMOD Scope 和 ZMOD AWG的Digilent 的 Eclypse Z7 。此外，為了生成信號(hào)并驗(yàn)證濾波器的響應(yīng)，我將使用全新的 Analog Discovery PRO 5250 ，它具有信號(hào) 125 Msps 信號(hào)發(fā)生器以及能夠獲取高達(dá) 1 Gsps 的兩通道示波器。

?

為了獲得自然離散頻率 ( ωd )，我將使用 ADP 5250 的一個(gè)非常有用的功能，它執(zhí)行頻率掃描，為每個(gè)頻率獲取濾波器的增益，因此結(jié)果將是濾波器的波特圖. 要配置此模式，我們需要打開工具網(wǎng)絡(luò)。

?

打開此工具后，我們需要配置點(diǎn)數(shù)，以及最小和最大頻率。我們還需要配置發(fā)生器將產(chǎn)生的信號(hào)幅度。在這種情況下，我的濾波器的固有頻率是 3 kHz，所以如果幅度為 1V，頻率范圍為 500 Hz 到 10 kHz，我將使用配置。點(diǎn)數(shù)將配置頻率步長(zhǎng)。就我而言，我已根據(jù)需要更改了點(diǎn)數(shù)。

首先，我將雙二階濾波器配置為帶阻巴特沃斯濾波器，固有頻率為 3 kHz。這次我沒(méi)有糾正頻率扭曲。然后，單擊 Single，ADC 5250 將執(zhí)行一個(gè)完整的交換。結(jié)果如下圖所示。

?

您可以看到濾波器的響應(yīng)如何是預(yù)期的，但自然頻率為 2.9407 khz，比預(yù)期頻率低 59.3 Hz。對(duì)頻率應(yīng)用預(yù)變形，濾波器的響應(yīng)是下一個(gè)（注意圖像的縮放不同）。

?

通過(guò)頻率校正，誤差如果為0.0052 Hz，這是一個(gè)卑鄙的誤差。在接下來(lái)的測(cè)試中，我將濾波器的品質(zhì)因數(shù)增加到 5。這將對(duì)陷波的帶寬產(chǎn)生影響，從而獲得更窄帶的濾波器。

?

雙線性變換因其簡(jiǎn)單性而被廣泛用于設(shè)計(jì) IIR 濾波器。我們只需要用 z 函數(shù)替換連續(xù)傳遞函數(shù)中的 s。其他方法如脈沖不變性需要部分分?jǐn)?shù)展開代數(shù)，所以過(guò)程有點(diǎn)復(fù)雜。此外，雙線性變換映射整個(gè) s 平面，因此沒(méi)有混疊誤差。另一方面，它會(huì)在頻率上產(chǎn)生扭曲，但使用預(yù)扭曲，正如我們?cè)谶@篇文章中看到的，我們可以修復(fù)該頻率偏差。