Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Comparing to Matlab's implementation #8

Closed
mischki opened this issue Oct 3, 2019 · 4 comments
Closed

Comparing to Matlab's implementation #8

mischki opened this issue Oct 3, 2019 · 4 comments
Assignees
Labels

Comments

@mischki
Copy link

mischki commented Oct 3, 2019

Hi
I have been trying to implement the padasip library to switch to python from Matlab. In order to do this I took the data generated and predicted by the Matlab example and ran them in python with padasip. As you can see from the attached plot the results are pretty different. (Matlab's output, y_mat (red trace) and padasip's output in purple.). Matlab seems to converge better using the same input parameters.

Could you please take a quick look and tell me if the code is correct? If so, what could be the difference in the algorithm?
I'm also attaching matlab's x, d and y arrays.

Thank you very much.

####### PYTHON CODE
%matplotlib notebook

import numpy as np
import matplotlib.pylab as plt
import padasip as pa 

d = np.load('d.npy')
y_mat = np.load('y_mat.npy')
x = np.load('x.npy')

mu = 0.8
m = 13

xx = pa.input_from_history(x, m) # reshaped in into a stack of filter-sized pieces 

# identification
f = pa.filters.FilterLMS(n=m, mu=mu)
y, e, w = f.run(d[m-1:], xx)

plt.figure(figsize=(18, 6))
plt.plot(x, label='x')
plt.plot(d, label = 'd', lw = 2)
plt.plot(y_mat, label= 'y_mat', lw=2)
plt.plot(y, label = 'y')
plt.tight_layout()
plt.legend()
%%%%% MATLAB CODE from here
filt = dsp.FIRFilter;
filt.Numerator = fircband(12,[0 0.4 0.5 1], [1 1 0 0], [1 0.2],... 
{'w' 'c'});
 
x = 0.1*randn(250,1);
n = 0.01*randn(250,1);
d = filt(x) + n;
 
mu = 0.8;
lms = dsp.LMSFilter(13,'StepSize',mu)
 
 
[y,e,w] = lms(x,d);
plot(1:250, [d,y,e])
title('System Identification of an FIR filter')
legend('Desired','Output','Error')
xlabel('Time index')
ylabel('Signal value')

matlab's d
matlab's x
matlab's y

image

@matousc89
Copy link
Owner

I have no idea.

I would suggest to plot history of the weights (or weights updates). And check if they start from the sam origin (zeros probably).

Also you can try to recalculate increments for few steps manually - print first few errors and weight increments to find if it is consistent with LMS rule.

Please keep me updated if you find anything.

@mischki
Copy link
Author

mischki commented Oct 3, 2019

OK, thanks, I'll try to debug step by step. I just wanted to confirm that the way I set it up was correct, particularly the reshaping of x using pa.input_from_history and f.run(d[m-1:], xx).

@matousc89
Copy link
Owner

Fair point. I have not check the validity of that line (I thought you are sure with this one). And probably this could be really the source of a problem.

You can easily check how is it done in Python, if you print d[m-1:][k] and xx[k] for range of numbers [1, 2, 3, 4, 5, 6] it should look like:

d[k] = 4, xx[k] = [1, 2, 3]
d[k+1] = 5, xx[k+1] = [2, 3, 4]

And it should work the same with your data:

xx[k] = [d[k-n], ... ,d[k-1]]

However, I am not sure how to check this in Matlab.

Let me know what you find, I am really curious.

@matousc89
Copy link
Owner

It seems that the Matlab fir filter you use is completely different (and more sophisticated), than the Padasip LMS implementation. Try comparison against the https://www.mathworks.com/help/dsp/ref/dsp.lmsfilter-system-object.html.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants