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

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.