I'm trying to implement a simulation of an ANC system with python, using this model 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') output_wf.setframerate(RATE) output_wf.setsampwidth(WIDTH) output_wf.setnchannels(CHANNELS) ## 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 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 filt.adapt(y_main_path,input_value) output = np.clip(output, -MAXVALUE, MAXVALUE) output = output.astype(int) binary_date = struct.pack('h', output) stream.write(binary_data) output_wf.writeframes(binary_data) binary_data = wf.readframes(1) stream.stop_stream() stream.close() p.terminate() # Close wavefiles wf.close() output_wf.close()