matousc89 / padasip Goto Github PK

Hi,

Can these filters be used for the change point detection?
Any examples?

Hi,
just a minor comment on the documentation. It states that there is an helper function to optimize the learning rate, but it is not implemented.

Hi,

I was running your RLS filter on large input and desired output signals (index length c.a. 700k) and noticed that the filter was pretty slow compared to Matlab (c.a. 1-2 hours vs. 5-10 seconds).

After looking through https://github.com/matousc89/padasip/blob/master/padasip/filters/rls.py, I noticed that there is a bottleneck in line 123

self.w_history[k:] = self.w

This line inserts w into the kth row and all the rows after the kth one for each k (time complexity O(N^2) where N is the signal length)

This bottleneck can be solved by adding a comma to line 123:

self.w_history[k, :] = self.w

that is, insert w only into the kth line in w_history (worst case time complexity O(M) where M is the filter length?)

After changing this line, the filter performed much faster (5-10 seconds like Matlab's RLS filter).

Hi,

I am working witch some ecg data and want to try and remove some motion artifacts from the signal using the accelerometers data. I am following this paper, which is written in matlab code. What I tried now is this

ecg = [ 5.12  5.28  5.41  5.56  5.71  5.84  6.01  6.14  6.28  6.43  6.6   6.72
  6.85  7.    7.14  7.3   7.43  7.55  7.7   7.82  7.96  8.42  8.74  8.31
  9.03 18.04 21.15 15.04 15.18 21.36 16.89  9.74  8.64 12.59 10.53 12.15
 26.53 40.02 41.19 28.55 20.38 24.8  18.96 10.35 11.43 18.26 13.43  8.16
  8.32 14.92 11.35  8.61  8.77 17.4  17.86 14.32 18.07 25.76 22.75 16.62
 16.55 23.32 18.36 10.54 10.07 16.75 12.41  9.06  9.19 15.95 13.62  9.39
 10.38 18.6  16.83 12.64 15.53 23.85 22.86 18.33 21.35 29.14 28.06 23.31
 24.62 29.94 23.98 16.58 19.56 27.18 25.26 20.2  20.86 27.13 24.88 19.51
 21.61 28.44 26.99 22.29]
a = [61.4        61.1429408  61.1560402  61.355666   61.6        61.7526485
 61.74060948 61.58011501 61.36       61.19293072 61.15967963 61.27348498
 61.48       61.68981649 61.82470523 61.8533976  61.8        61.72347136
 61.68157705 61.69961973 61.76       61.81628039 61.82203364 61.75755968
 61.64       61.51250152 61.41988386 61.3851854  61.4        61.43297464
 61.45036915 61.43613901 61.4        61.369451   61.37121731 61.41355133
 61.48       61.53830016 61.55913501 61.53373826 61.48       61.43333539
 61.42728887 61.47470248 61.56       61.64694723 61.69766599 61.69250612
 61.64       61.57165315 61.52462153 61.52179827 61.56       61.61231643
 61.64294654 61.62627121 61.56       61.46581015 61.37779666 61.32541007
 61.32       61.35153414 61.39621507 61.42994451 61.44       61.42904883
 61.41041332 61.39842472 61.4        61.41216749 61.42625596 61.4353692
 61.44       61.44810346 61.46960493 61.50887311 61.56       61.60805922
 61.6359193  61.6328438  61.6        61.54975208 61.49917066 61.46144961
 61.44       61.42831401 61.41529711 61.39266702 61.36       61.32459801
 61.29654143 61.28215882 61.28       61.28177902 61.27767054 61.26266298
 61.24       61.21947606 61.21145686 61.22012118]
Hd = pa.filters.FilterLMS(n=32, mu=0.0008)
y, e, w = Hd.run(a, ecg)

Where a is a list of the summation of the accelerometer data and ecg is a list of ecg data (added some of it if someone want to test).
I keep getting a typeError when running Hd.run(a,ecg) (TypeError: object of type 'numpy.int64' has no len())

I am a bit unsure if it is the paper I have misunderstood or the use of the adaptive filter here, but if someone has any inputs that would be nice.

Comparing the code for the R update with the mathematical description, the following lines look fishy.

FilterRLS.adapt at rls.py:L220:

R1 = np.dot(np.dot(np.dot(self.R,x),x.T),self.R)

and FilterRLS.run at rls.py:L268:

R1 = np.dot(np.dot(np.dot(self.R,x[k]),x[k].T),self.R)

These lines are supposed to calculate the conjugation of the outer product between a columns vector x and it's transpose (a row vector), multiplied from both sides by the matrix R. However, for a 1D array x in the first case or x[k] in the second case, x.T and x[k].T are still the exact same 1D array. (Contrary to MATLAB, numpy doesn't silently make all 1D arrays into a column or row vector). Thus, the R-update is likely wrong.

A correct way to write that code could be

R1 = self.R @ (x[:,None]*x[None,:]) @ self.R

and

R1 = self.R @ (x[k,:,None]*x[k,None,:]) @ self.R

respectively.

If that misconception about 1D arrays in numpy made into other algorithms, it could likely be the cause of #8.

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

Suppose I have a pre-recorded sample audio file x and have a live input audio signal y. Can I use padasip in order to cancel x from y in real time?

Hello, I'm really new to signal processing and was referencing this guide: https://github.com/matousc89/Python-Adaptive-Signal-Processing-Handbook/blob/master/notebooks/ALE/ALE.ipynb. However whenever I run my program which I'll post below, I get the error: IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed. So I was wondering how to fix it or maybe how to go about reading wav files differently to work with Padasip.

from scipy.io import wavfile
import padasip as pa

N = 100
MU = 0.01
INPUT_FILE = "voices.wav"
NOISE_FILE = "noise.wav"

def read(fn):
fs, data = wavfile.read(fn)
x = data[:,0].copy()
x = x.astype("float64")
return x

def cancel():
x = read(INPUT_FILE)
v = read(NOISE_FILE)
d = x + v
f = pa.filters.FilterLMS(n=N, mu=MU)
y, e, w = f.run(d, x)

def main():
cancel()

if name == "main":
main()

Could you please help to know if input vector should be resized into fixed batch size during training?

-----------Error observed in the following setup-----------
Setup:
1.Input vector of shape:(140928,)
2.Length of filter: 16
3. Desired vector of shape: (140928,)

In step 3: After resizing either row or column changes size.

###This contradicts either of the requirement
###1. Number of rows equals to number of input shape
###3. Number of columns equals to the length of filter.
###4. length of input vector and desired vector

4. Resizing the input vector: 1. Transpose 2. Resize into (140928, 140928/16)

Observation:
In run() method following

So basically, I presume that before feeding input vector directly, it is necessary to resize into batch size meeting the requirement of step 3.
OR,
What is a correct approach to carry out the simulation?

I believe there is an error in the adapt() and run() methods of FilterRLS (lines 220 and 268 of rls.py). Lines 220-222 are currently written as:

R1 = np.dot(np.dot(np.dot(self.R,x),x.T),self.R)
R2 = self.mu + np.dot(np.dot(x,self.R),x.T)
self.R = 1/self.mu * (self.R - R1/R2)

The matrix R1 should be R*x*x^T*R. In the current implementation, the second product (the product of the vectors R*x and x^T) is a dot product, when it should be a matrix product.

This means the product of the first three matrices (R*x*x^T) is being incorrectly calculated as a scalar, not a matrix, and the matrix R1 will always be a scalar times R, so R1 (and therefore R) will always be diagonal matrices.

Tested with python version 2.7.10 and numpy version 1.13.0rc2.

When using the adapt function of the AdaptiveFilter class

a TypeError is thrown:

---

TypeError Traceback (most recent call last)
/tmp/ipykernel_10111/2229501335.py in
26 print(f'measurement: {x}')
27
---> 28 f.adapt(x, buffer)
29 buffer = np.append(buffer, x)[1:]
30 measurements.append(x)

~/.local/lib/python3.8/site-packages/padasip/filters/base_filter.py in adapt(self, d, x)
164 y = self.predict(x)
165 e = d - y
--> 166 self.w += self.learning_rule(self, e, x)
167
168 def run(self, d, x):

TypeError: learning_rule() takes 3 positional arguments but 4 were given

The self attribute in calling the learning_rule function shouldn't be needed.

Hi,

Thanks for this excellent adaptive signal processing package. I am new to this and I would like to implement Adaptive Line Enhancement (ALE) using the Padasip LMS filter for denoising the data available in a 1D array.

ALE implemented using dsp.LMSFilter in MATLAB:

It would be helpful if you could provide me some pointers for replicating dsp.LMSFilter available in Matlab using Padasip. thanks in advance.

Great module, thank you for awesome job! I'm interested in recursive regression. I wonder if Multi-layer Perceptron can be applied recursively or will get extended to recursive application at some point? As I see, at the moment there's only autoregressive option padasip.ann.mlp.NetworkMLP.update(d) and nothing like padasip.filters.rls.FilterRLS.adapt(d,x).

Hello @matousc89, I am new to signal processing so I have very little grasp on theoretical concepts of signal, speech and audio processing.

Consider this, I have a audio or maybe a live stream of audio which is spitted into 5 seconds of audio data, in the end I have a audio data which needs to noise free for an application. I want noise to removed adaptively. I am loading file with librosa, how can I use your code to integrate my use case?

Hey,
I am using the RLS adaptive filtering technique for noise cancellation from the noisy desired signal with reference noise input.

f = pa.filters.FilterRLS(mu=0.9, n=n)
y, e, w = f.run(d, x)

Could you please clarify the terms y,e, and w? I understand that y is the filtered output of reference noise, e is the output error between filtered y and desired signal d and w is the weights. So, the error signal e is the noiseless desired signal we required right?

Thanks in advance :)

Hi,

I used filters.FilterRLS. The initial y data was 0 and x was an array of 0s. RLS.w was set to 0 before adapting. However, RLS.R and RLS.w became very big after calling adapt(y,x), at the order of 6.29647980e+35 and 6.43140682e+224, respectively. A few adaptive steps later, w became nan. Setting both of them to 0 by RLS.R *= 0 and RLS.w *= 0 before adapting solved this problem, but I am not sure whether or not this leads to poor fit.

I am dealing with this on a remote server so it is not convenient to provide details. Sorry for that. Thanks!

Hi, first of all congratulation for the implementations that you have done ;)

I was analysing your codes, and I saw something that it might be an issue.
For example, in the NLMS I think that the correct implementation should be:

    for k in range(L, N-L):
        y = np.dot(self.w, np.flipud(x[k-L:k]))
         ...

Where L is the size of the filter.

https://en.wikipedia.org/wiki/Finite_impulse_response

Looking forward for your reply