MATLAB fftfilt equivalent for Python

Question

I am trying to traslate the following function created in MATLAB into Python,

function output_phase = fix_phasedata180(phase_data_degrees, averaging_length)

x = exp(sqrt(-1)*phase_data_degrees*2/180*pi);
N = averaging_length;
b = 1/sqrt(N)*ones(1,N);
y = fftfilt(b,x);y = fftfilt(b,y(end:-1:1));y = y(end:-1:1); # This is a quick implementation of filtfilt using fftfilt instead of filter
output_phase = (phase_data_degrees-(round(mod(phase_data_degrees/180*pi-unwrap(angle(y))/2,2*pi)*180/pi/180)*180));
temp = mod(output_phase(1),90);
output_phase = output_phase-output_phase(1)+temp;
output_phase = mod(output_phase,360);
s = find(output_phase>= 180);
output_phase(s) = output_phase(s)-360;

So, I am trying to implement this function defined in MATLAB into Python here

def fix_phasedata180(data_phase, averaging_length):
    x = np.exp(1j*data_phase*2./180.*np.pi)
    N = averaging_length
    b = 1./np.sqrt(N)*np.ones(N)
    y = fftfilt(b,x)          
    y = fftfilt(b,y[::-1])
    y = y[::-1]
    output_phase = data_phase - np.array(map(round,((data_phase/180.*np.pi-np.unwrap(np.angle(y))/2.)%(2.*np.pi))*180./np.pi/180.))*180
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

I was thinking that the function fftfilt was a clone of fftfilt in MATLAB, when I run I have the following error

ValueError                                Traceback (most recent call last)
<ipython-input-40-eb6944fd1053> in <module>()
      4 N = averaging_length
      5 b = 1./np.sqrt(N)*np.ones(N)
----> 6 y = fftfilt(b,x)

D:/folder/fftfilt.pyc in fftfilt(b, x, *n)
     66         k = min([i+N_fft,N_x])
     67         yt = ifft(fft(x[i:il],N_fft)*H,N_fft) # Overlap..
---> 68         y[i:k] = y[i:k] + yt[:k-i]            # and add
     69         i += L
     70     return y

ValueError: could not broadcast input array from shape (0,0) into shape (0)

So, my question is: are there any equivalent to MATLAB fftfilt in Python? The aim of my function output_phase is to correct the fast variations in a phase signal and then correct n*90 degrees shifts, showed bellow

Answer 1

I also had issues when converting a MATLAB code. I went from this MATLAB code:

signal_weighted = fftfilt( weight, signal.^2 ) / Ntau;

to this python code:

from scipy.signal import convolve

signal_weighted = convolve(signal**2 ,weightData, 'full', 'direct') / Ntau
signal_weighted = signal_weighted[:len(signal)]

If you want something faster than convolve, see this overlap and add fft implementation

Answer 2

This is a bit late but I found the answer to this while translating some matlab code of my own.

TLDR: Use mode="full" for any of the convolve functions in scipy.signal

I leaned on scipy's recipes to guide me through this. The rest of my answer is effectively a summary of that page. Matlabs fftfilt function can be replaced with any of the convolve functions mentioned in the cookbook (np.convolve, scipy.signal.convolve, .oaconvolve, .fttconvolve), if you pass mode='full'.

import numpy as np
from numpy import convolve as np_convolve
from scipy.signal import fftconvolve, lfilter, firwin
from scipy.signal import convolve as sig_convolve

# Create the m by n data to be filtered.
m = 1
n = 2 ** 18
x = np.random.random(size=(m, n))

ntaps_list = 2 ** np.arange(2, 14)
for ntaps in ntaps_list:
    # Create a FIR filter.
    b = firwin(ntaps, [0.05, 0.95], width=0.05, pass_zero=False)
    conv_result = sig_convolve(x, b[np.newaxis, :], mode='full')

Happy filtering!

Answer 3

Finally, I got an improvement in my code. I replace the fftfilt (twice applied) by the scipy.signal.filtfilt (that is basically the same). So my code traslated into python will be:

import numpy as np
import scipy.signal as sg

AveragingLengthAmp = 10
AveragingLengthPhase = 10
PhaseFixLength = 60
averaging_length = channel_sampling_freq1*PhaseFixLength

def fix_phasedata180(data_phase, averaging_length):
    data_phase = np.reshape(data_phase,len(data_phase))
    x = np.exp(1j*data_phase*2./180.*np.pi)
    N = float(averaging_length)
    b, a = sg.butter(10, 1./np.sqrt(N))
    y = sg.filtfilt(b, a, x)
    output_phase = data_phase - np.array(map(round,((data_phase/180*np.pi-np.unwrap(np.angle(y))/2)%(2*np.pi))*180/np.pi/180))*180
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

out1 = fix_phasedata180(data_phase, averaging_length)

def fix_phasedata90(data_phase, averaging_length):
    data_phase = np.reshape(data_phase,len(data_phase))
    x = np.exp(1j*data_phase*4./180.*np.pi)
    N = float(averaging_length)
    b, a = sg.butter(10, 1./np.sqrt(N))
    y = sg.filtfilt(b, a, x)
    output_phase = data_phase - np.array(map(round,((data_phase/180*np.pi-np.unwrap(np.angle(y))/4)%(2*np.pi))*180/np.pi/90))*90
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    output_phase = output_phase%360
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

offset = 0
data_phase_unwrapped = np.zeros(len(out2))
data_phase_unwrapped[0] = out2[0]
for jj in range(1,len(out2)):
    if out2[jj]-out2[jj-1] > 180:
        offset = offset + 360
    elif out2[jj]-out2[jj-1] < -180:
        offset = offset - 360
    data_phase_unwrapped[jj] = out2[jj] - offset

Here fix_phasedata180 fix the 180-degrees shifts, similarly for fix_phasedata90. The channel_sampling_freq1 is 1/sec.

The result is:

that is mostly right. Only I have some question understanding the scipy.signal.butter and scipy.signal.filtfilt. As you see, I choose:

b, a = sg.butter(10, 1./np.sqrt(N))

Here the order of the filter (N) is 10 and the critical frequency (Wn) is 1/sqrt(60). My question is, How can I choose the appropiated order of the filter? I tried since N=1 until N=21, larger than 21 the result data_phase_unwrapped are all NAN. I tried too, giving values for padlen in filtfilt but I didnt understand it well.

Answer 4

The function you linked to is a Python equivalent to the Matlab function. It just happens to be broken.

Anyway, MNE also has an implementation of the overlap and add method used by the fftfilt function. It's a private function of the library, and I'm not sure if you can call it exactly equivalent to the Matlab function, but maybe it's useful. You can find the source code here: https://github.com/mne-tools/mne-python/blob/master/mne/filter.py#L41.