Biquad filter produces garbage (solved)

Support
1

I’ve been having trouble with the stability of my biquad filter implementation for quite some time now. Today I finally created the following minimal test to find out what is going on. For some reason, looping with biquad1 works just fine, while the MyHDL implementation (biquad2) produces garbage. Can someone spot what’s wrong here? I’m using MyHDL 0.9.0.

from myhdl import instance, always_seq, intbv, modbv, delay
from myhdl import ResetSignal, Signal, Simulation, StopSimulation
import matplotlib.pyplot as plt
import numpy as np

b0, b1, b2, a1, a2 = 0x00083371, 0x001066e2, 0x00083371, 0xc172a282, 0x1e9590ee
    
def biquad1(y1, y2):
    """Step response for direct form I biquad filter"""
    return int(intbv(((b0 + b1 + b2) >> 1) -
                     ((intbv(a1)[32:].signed() * intbv(int(y1))[32:].signed()) >> 29) -
                     ((intbv(a2)[32:].signed() * intbv(int(y2))[32:].signed()) >> 29))[32:])

def biquad2(clk, reset, y):
    """Step response for direct form I biquad filter (MyHDL)"""
    y1 = Signal(y.val)
    y2 = Signal(y.val)
    
    @always_seq(clk.posedge, reset=reset)
    def logic():
        y1.next = y
        y2.next = y1
        y.next = int(intbv(((b0 + b1 + b2) >> 1) -
                     ((intbv(a1)[32:].signed() * intbv(int(y1))[32:].signed()) >> 29) -
                     ((intbv(a2)[32:].signed() * intbv(int(y2))[32:].signed()) >> 29))[32:])

    return logic

def testbench():
    clk = Signal(bool(0))
    reset = ResetSignal(0, active=1, async=True)
    out = Signal(modbv(0)[32:])
    filter = biquad2(clk, reset, out)
    
    @instance
    def generate():
        y = np.zeros(250, dtype='uint32')
        for i in range(500):
            yield delay(1)
            clk.next = not clk
            if clk:
                y[i // 2] = out
        plt.plot(y / 2**30)
        raise StopSimulation
            
    return generate, filter

if __name__ == "__main__":
    tb = testbench()
    Simulation(tb).run()

    y = np.zeros(252, dtype='uint32')
    for i in range(2, 252):
        y[i] = biquad1(y[i-1], y[i-2])
    plt.plot(y[2:] / 2**30)

I know the y.next assignment looks funny, but it is the way it is to make the two functions as identical as possible.

2

I think this does what you intended:

def biquad2(clk, reset, y):
    """Step response for direct form I biquad filter (MyHDL)"""
    y1 = Signal(y.val)
    y2 = Signal(y.val)
    
    @always_seq(clk.posedge, reset=reset)
    def logic():
#         y1.next = y
        y2.next = y1
        y1.next = int(intbv(((b0 + b1 + b2) >> 1) - 
                     ((intbv(a1)[32:].signed() * intbv(int(y1))[32:].signed()) >> 29) - 
                     ((intbv(a2)[32:].signed() * intbv(int(y2))[32:].signed()) >> 29))[32:])

    @always_comb
    def assign():
        y.next = y1
        
    return logic, assign
3

Thanks! Indeed it does, although I can’t quite wrap my head around why the two aren’t equivalent. Makes me wonder whether I have similar structures in my other modules.

4

You simulation fills the result vector at both clock edges:
This produces two identical plots:

'''
Created on 29 mrt. 2018

@author: mhavu / josy
'''

from myhdl import instance, always_seq, intbv, modbv, delay
from myhdl import ResetSignal, Signal, Simulation, StopSimulation
import matplotlib.pyplot as plt
import numpy as np
from myhdl._always_comb import always_comb

b0, b1, b2, a1, a2 = 0x00083371, 0x001066e2, 0x00083371, 0xc172a282, 0x1e9590ee

    
def biquad1(y1, y2):
    """Step response for direct form I biquad filter"""
    return int(intbv(((b0 + b1 + b2) >> 1) - 
                     ((intbv(a1)[32:].signed() * intbv(int(y1))[32:].signed()) >> 29) - 
                     ((intbv(a2)[32:].signed() * intbv(int(y2))[32:].signed()) >> 29))[32:])


def biquad2(clk, reset, y):
    """Step response for direct form I biquad filter (MyHDL)"""
    y1 = Signal(y.val)
    y2 = Signal(y.val)
    
    @always_seq(clk.posedge, reset=reset)
    def logic():
#         y1.next = y
        y2.next = y1
        y1.next = int(intbv(((b0 + b1 + b2) >> 1) - 
                     ((intbv(a1)[32:].signed() * intbv(int(y1))[32:].signed()) >> 29) - 
                     ((intbv(a2)[32:].signed() * intbv(int(y2))[32:].signed()) >> 29))[32:])

    @always_comb
    def assign():
        y.next = y1
        
    return logic, assign


def testbench():
    clk = Signal(bool(0))
    reset = ResetSignal(0, active=1, async=True)
    out = Signal(modbv(0)[32:])
    
    bqf = biquad2(clk, reset, out)
    
    @instance
    def genclk():
        while True:
            yield delay(1)
            clk.next = not clk
        
    @instance
    def generate():
        y = np.zeros(250, dtype='uint32')
        for i in range(250):
            yield clk.posedge
            y[i] = out
        plt.plot(y / 2 ** 30)
        plt.show()

        raise StopSimulation
            
    return genclk, generate, bqf


if __name__ == "__main__":
    tb = testbench()
    Simulation(tb).run()

    y = np.zeros(250, dtype='uint32')
    for i in range(2, 250):
        y[i] = biquad1(y[i - 1], y[i - 2])
    plt.plot(y / 2 ** 30)
    plt.show()
5

Yeah, I noticed that when I tried your answer. Hence I edited the plotting code before I responded. I’m still having trouble understanding why the @always_comb block is necessary, and why the two aren’t equivalent.

6

In the following picture, the top diagram is a representation of your code. The bottom diagram is the representation of Josy code.
image

7

Thanks for the wonderful diagrams @DrPi! So there was an extra delay there. I suspected it, but couldn’t see where, no matter how hard I thought about it. Goes to show once again how useful schematics and diagrams are.