描述

濾波可能是任何嵌入式工程師必須設計的最常見的 DSP 算法，無論他們是為 STM32、TI DSP 還是為 FPGA 開發的算法。濾波非常重要，因為大多數應用程序都連接到現實世界，因此它們需要捕獲具有所有真實信號都具有的不良元素的真實信號，即噪聲、偏移……此外，我們需要獲取的信號，并非總是如此將是我們獲取的主要信號，例如，對于生物信號，我們將獲取與電網頻率相對應的 50 或 60 Hz 的高電平。在這種情況下，我們將需要一個過濾器。

作為本博客的讀者，您肯定知道，在數字系統上，我們有 2 種濾波器，FIR（有限脈沖響應）濾波器，它只使用現在輸入的值和過去的值，以及 IIR（無限脈沖Response），它使用現在和過去的輸入值，以及過去輸出的值。FIR 濾波器很容易為許多應用設計，但它們的響應對于低階濾波器是有限的。另一方面，使用 IIR 濾波器，我們可以實現更積極的響應，但它們的缺點是濾波器在某些情況下可能不穩定。在本文中，我們將了解如何逐步實現二階濾波器，從濾波器設計到 Verilog 實現。

設計 IIR 濾波器的過程始終相同。首先我們要使用它的連續傳遞函數來設計濾波器，然后，一旦選擇了濾波器的固有頻率、階數和品質因數，我們將得到相應的 s

功能。為了將該濾波器數字化，下一步是應用一些變換來用 z 替換變量s ，這個過程稱為離散化。為了離散化連續傳遞函數，我們可以使用多種方法，其中一種易于使用的方法是Tustin或雙線性變換。

雙線性變換基于將連續傳遞函數中的s變量替換為z的函數。在下一個示例中，我們可以看到應用雙線性變換對單極點低通濾波器進行離散化。接下來是我們必須在連續傳遞函數中執行的替換。

對于單極點濾波器，這些將是執行雙線性變換的步驟。低通單極濾波器的方程如下：

應用轉換，我們有：

與分母運算以獲得共同的分母，

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

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

接下來是 100 Hz 低通濾波器的整個 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')

赫茲的結果是

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

為了驗證我們手動獲得的方程是否正確，我們必須使用我們之前獲得的值創建一個傳遞函數。

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

這種情況下的結果是

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

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

雙線性變換與其他任何離散化方法一樣，都是一種近似，因此我們將引入一些誤差。該誤差的大小與采樣頻率有關，因為在某些情況下，非常（非常）高的采樣頻率可能看起來是連續的。在實際系統中，采樣頻率取決于 ADC 和許多實際組件，我們可以使用高采樣頻率，但在其他情況下，我們會引入顯著誤差。

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

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

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

現在，我們將使用半角分割指數。

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

現在，我們將使用半角分割指數。

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

然后，我們必須找到與類比的關系，所以我們要使兩個方程相等

我們可以驗證實部為 0，而虛部wwa直接為

請注意，當ts趨于零時，極限是ωa = ωd。

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

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');

我們可以看到數字濾波器的頻率是如何移動的。然后我們將對自然頻率應用預變形，并使用這個新頻率重新離散化濾波器。

% 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');

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

為了在現實世界中驗證這種行為，我們將在 Verilog 中實現這個二階帶阻濾波器，然后我們將檢查扭曲錯誤，然后如何使用預扭曲來修復它。

在開始編碼之前，我們必須考慮我們想要創建的模塊是什么。在我的情況下，能夠參數化模塊是非常重要的，這樣，如果我改變輸入的寬度，或者輸出的寬度，那只代表一個參數的改變，而不需要修改模塊的代碼因為這意味著模塊的新測試。此外，我們必須考慮我們將使用什么數字格式。對于這種過濾器，很明顯我們需要有符號格式，以及管理十進制數的能力. 為了達到這個要求，我們必須在定點和浮點之間進行選擇，對我來說決定很明確，我想要一個可以單獨實例化的代碼，沒有外部浮點單元，所以格式將是定點的。如果我們加上參數的要求，和定點格式，我們需要參數來定義信號的寬度，還有小數部分的寬度，這與我們需要的精度有關。此外，將選擇精度以獲得與設計盡可能相似的實施濾波器的響應。現在，我們需要參數化多少寬度？我們可以定義一個寬度，用于輸入和輸出以及系數，但是這樣，系數的寬度將根據數據生成器的寬度來選擇，對于某些過濾器，系數在穩定極限上，這可能代表一個問題。因此，至少我們將定義 2 種不同的寬度，一種用于輸入和輸出，根據數據消耗和數據源，另一種根據我們需要的分辨率。另一個我們要考慮寬度的事情，是內部運算的分辨率，因為在某些情況下，穩定性問題與系數本身的寬度無關，而是運算分辨率，所以在這一點上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數的寬度，以及內部濾波器操作的寬度。最后，過濾器參數將如下所示。一個用于輸入和輸出，根據數據消耗和數據源，另一個根據我們需要的分辨率。另一個我們要考慮寬度的事情，是內部運算的分辨率，因為在某些情況下，穩定性問題與系數本身的寬度無關，而是運算分辨率，所以在這一點上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數的寬度，以及內部濾波器操作的寬度。最后，過濾器參數將如下所示。一個用于輸入和輸出，根據數據消耗和數據源，另一個根據我們需要的分辨率。另一個我們要考慮寬度的事情，是內部運算的分辨率，因為在某些情況下，穩定性問題與系數本身的寬度無關，而是運算分辨率，所以在這一點上，將是有趣的解耦將在 MATLAB 或 Python 上生成的系數的寬度，以及內部濾波器操作的寬度。最后，過濾器參數將如下所示。穩定性問題與系數本身的寬度無關，而與運算分辨率有關，因此在這一點上，將在 MATLAB 或 Python 上生成的系數的寬度與內部濾波器操作的寬度解耦會很有趣. 最后，過濾器參數將如下所示。穩定性問題與系數本身的寬度無關，而與運算分辨率有關，因此在這一點上，將在 MATLAB 或 Python 上生成的系數的寬度與內部濾波器操作的寬度解耦會很有趣. 最后，過濾器參數將如下所示。

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

接下來我們必須考慮接口。在模塊之間以連續方式傳輸數據的應用中，AXI4-Stream 將是最佳選擇。在輸入和輸出字段中，我們將定義主從 AXI4-Stream 接口來獲取和發送數據。盡管輸入和輸出寬度是參數化的，但如果我們想將模塊連接到現有的 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,

關于系數，插入它們最好是使用 AXI4-Lite 接口，但我想設計這個模塊而不需要使用 Zynq 或 Microblaze，所以插入系數的簡單方法是作為輸入。

/* 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

現在，在定義了所有輸入和輸出之后，我們必須考慮數據流。首先，由于我們定義了不同的格式，為了能夠在系數和數據之間進行運算，我們必須將所有信號的格式更改為內部格式，內部格式由內部寬度和小數寬度定義。使用的整數寬度被定義為一個localparam，并且對應于寬度和小數寬度之間的差異。要更改格式，我們將在信號的低端填充零，直到完成小數部分。對于整數部分，由于使用的格式是有符號的，我們必須進行符號擴展，即填充 MSb 的值，直到整數部分完成。請注意，這是有效的，因為內部寬度大于或等于輸入輸出和系數寬度。

/* 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}} };

數字濾波器，無論是 FIR 還是 IIR，都需要存儲過去輸入的值，而 IIR 還需要存儲過去輸出的值，因此我們需要一個流水線結構來存儲過去的值。

/* 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

現在，接下來是執行過濾器計算。一個二階濾波器需要執行 5 次乘法，記住這并不意味著模塊使用 5 個 DSP slice。我開發的代碼執行組合乘法。這允許 0 個時鐘周期延遲，但會限制濾波器的速度。如果您的時序約束沒有得到滿足，您可以注冊乘法器的輸入，然后進行重新定時，讓 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;

乘法的結果將在輸入的情況下相加，在輸出的情況下相減，以獲得濾波器輸出。由于乘積的輸出大小是操作數的兩倍，因此加法必須是相同的寬度。要在乘法上使用輸出，我們必須執行移位以刪除額外的小數位。關于額外的整數位置將在分配時被截斷。

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);

對于 AXI4-Stream 管理信號，由于濾波器作為橋接器，有 1 個周期延遲，因此管理信號會以 1 個周期延遲通過濾波器。

至此，我們已經實現了 Verilog 模塊并準備好用作帶阻濾波器。接下來是使用該tfdata函數在 MATLAB 中獲取系數。關于過濾器的寬度，我們將使用 18 位分辨率用于小數部分，2 用于整數部分，以使過濾器達到 +-2。接下來是獲取定點系數的 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)

這將返回下一個結果。

b0_qq =

      252631


b1_qq =

     -468081


b2_qq =

      252631


a1_qq =

     -468081


a2_qq =

      243119

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

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)
  );

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

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

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

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

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

通過頻率校正，誤差如果為0.0052 Hz，這是一個卑鄙的誤差。在接下來的測試中，我將濾波器的品質因數增加到 5。這將對陷波的帶寬產生影響，從而獲得更窄帶的濾波器。

雙線性變換因其簡單性而被廣泛用于設計 IIR 濾波器。我們只需要用 z 函數替換連續傳遞函數中的 s。其他方法如脈沖不變性需要部分分數展開代數，所以過程有點復雜。此外，雙線性變換映射整個 s 平面，因此沒有混疊誤差。另一方面，它會在頻率上產生扭曲，但使用預扭曲，正如我們在這篇文章中看到的，我們可以修復該頻率偏差。