Least Sqquares FIR

least_squares_fir.py

import numpy as np
from scipy.linalg import solve_toeplitz
######
##
## Algorithm that assumes 
## a fixed length impulse response (length imp_length)
## and calculates the least-squares solution of the 
## convolution equation y = convolve(x, w). With optional l2-regularization.
##
## Arguments:
## reg - l2 regulariztion strength (optional).
## x - input signal length N
## y - output signal length N
## imp_length - desired impulse response length.
##
## Returns:
## w - Learned impulse response (length imp_length)
##
## Remember to zero-pad x/y if normal convolution (not circular) is desired.
## This is typically the case.  
    
def least_squares_fir(x,y,imp_length,reg=0.0):

    n=len(x)
    X = np.fft.rfft(x)
    Y = np.fft.rfft(y)
    
    b = np.fft.irfft(np.conj(X)*Y,n)[0:imp_length] # Normal equation vec b
    
    tmp_a = np.fft.irfft(np.conj(X)*X,n) # Normal Equation matrix A is first (imp_length x imp_length) elements of the circulant matrix made from this vector.
    
    c = tmp_a[0:imp_length]
    r = tmp_a[::-1][0:imp_length-1]
    r = np.concatenate((tmp_a[[0]],r))

    r[0]+=reg 
    c[0]+=reg
    
    return solve_toeplitz( (c,r), b)