I'm trying to implement a simulation of an ANC system with python, using this model here. enter image description here

My simulation keeps diverging, and I honestly don't know why. I'm using a source for LMS adaptive filter from Mathworks here. When I comment out the LMS update function, the system is stable, and the output is exactly the input. But when I plug in the adaptive filter update, the system starts to diverge. So I thought the problem is with my LMS update. And I implemented the LMS adaptive filter with padasip package. The system is still diverging. Now I honestly don't know how I copied the model from Matlab wrong. Can someone help? Though It's not a syntax error, I pasted my code below for reference. I used all difference equations for every system.

Some values I used in this, The input I'm using is a white noise file I generated from Matlab, 1dB power and sampling frequency of 16000 Hz.

b for S(z) is [0.5,0.5,-.3,-.3,-.2,-.2,]

b for S_est(z) is [0.466,0.533,-0.257,-0.274,-0.231,-0.175]

b for Main path S(z) is [0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300]

All I copied from the Simulink simulation of the same system

#ANC Simulation Main

import pyaudio, wave, struct, math
import numpy as np
from matplotlib import pyplot as plt
import padasip as pa 

## Variables Setup
mu = 0.1
MAXVALUE = 2**15-1  # Maximum allowed output signal value (because WIDTH = 2)

# Initialization of adaptive weight w
w = np.zeros(13)

# Main Path Filters
order_path = 12
b_path = np.array([0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300])
x = np.zeros(13)
y_main_path =0 

#Adapt filter
order_adapt = 12
y_adapt =np.zeros(13)

filt = pa.filters.FilterLMS(13, mu=mu)

#second path
order_sec_path = 12
b_sec_path = [0.5,0.5,-.3,-.3,-.2,-.2,0,0,0,0,0,0,0]

y_secpath = 0

#est sec path
order_est_sec_path = 12
b_est_sec_path = np.array([0.466,0.533,-0.257,-0.274,-0.231,-0.175,0,0,0,0,0,0,0])

est_sec_out = 0
y_estsecpath = np.zeros(13)

# File names
wavfile = 'matlab_1db_wn.wav'
output_wavfile = 'ANC_Result.wav'

## Read WAV file
wf = wave.open(wavfile,'rb')

CHANNELS        = wf.getnchannels()     # Number of channels
RATE            = wf.getframerate()     # Sampling rate (frames/second)
signal_length   = wf.getnframes()       # Signal length
WIDTH           = wf.getsampwidth()     # Number of bytes per sample

print('The file has %d channel(s).'            % CHANNELS)
print('The frame rate is %d frames/second.'    % RATE)
print('The file has %d frames.'                % signal_length)
print('There are %d bytes per sample.'         % WIDTH)

# Read first BLOCKLEN
binary_data = wf.readframes(1)

## Output WAV file
output_wf = wave.open(output_wavfile, 'w')

## Open audio stream
p = pyaudio.PyAudio()
stream = p.open(
   format      = p.get_format_from_width(WIDTH),
   channels    = CHANNELS,
   rate        = RATE,
   input       = False,
   output      = True )

## Main Loop
while len(binary_data) > 0:

   # convert binary data to numbers
   input_block = struct.unpack('h', binary_data)
   input_value = input_block[0]

   x = np.delete(x,-1)
   x = np.insert(x,0,input_value)
   y_main_path = np.dot(b_path,np.transpose(x))
   y_adapt = filt.predict(input_value)

   #y_adapt = np.delete(y_adapt,-1)
   #y_adapt = np.insert(y_adapt,0,adapt_out)

   est_sec_out = np.dot(b_est_sec_path,np.transpose(x))
   y_estsecpath = np.delete(y_estsecpath,-1)
   y_estsecpath = np.insert(y_estsecpath,0,est_sec_out)

   #sec path
   y_secpath = np.dot(b_sec_path,np.transpose(y_adapt))

   output = y_main_path - input_value

   # LMS update

   output = np.clip(output, -MAXVALUE, MAXVALUE)
   output = output.astype(int)
   binary_date = struct.pack('h', output)
   binary_data = wf.readframes(1)

# Close wavefiles

1 Answer 1


Turns out it's just a scaling issue. Because the struck.unpack() gives out a signed 16-bit value, the step size is too large for it. Two ways to go around this, you scale down your input and scale back up when outputting, or you just adjust the step size mu.

I used a mu value of 10e-8 and seems like this value works differently in different environments, reasons unknown. With jupyter-notebook is around 10e-7, python idle is 10e-8, and Matlab is 10e-9. Just tweak around that range with the code I posted above and should give you a converged result.

Hope this helps whoever is reading.

edit; for MatLab, I scaled the input by 32768 so the mu value is 10e-9, if max amplitude is 1, use 0.1.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.