对应Python的MATLAB fftfilt

Question

我正试图将以下在MATLAB中创建的函数转换为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;

所以，我试图在Python中实现这个在MATLAB中定义的函数，如下所示

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

我在想，这个功能 fftfilt 是MATLAB中fftfilt的克隆，当我运行时，我有以下错误

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)

那么，我的问题是：在Python中有没有与MATLAB fftfilt等价的函数？我的函数的目的是 output_phase 是修正相位信号的快速变化，然后修正n*90度的偏移，如图所示。

Answer 1

你链接到的函数 是 相当于 Matlab 函数的 Python。它只是碰巧坏了。

总之。MNE 也有一个执行。叠加 fftfilt函数使用的方法。这是一个库的私有函数，我不知道你是否可以调用它完全等同于Matlab函数，但也许它是有用的。你可以在这里找到源代码。https:/github.commne-toolsmne-pythonblobmastermnefilter.py#L41。.

Answer 2

最后，我的代码有了改进。我把 fftfilt (两次应用)，由 scipy.signal.filtfilt (基本上是一样的)。所以我的代码traslated into python将是。

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

这里 fix_phasedata180 修正180度的偏移，同样地，对于 fix_phasedata90. 该 channel_sampling_freq1 是1秒。

结果是。

基本上是正确的 只不过我有一些问题需要理解 scipy.signal.butter 和 scipy.signal.filtfilt。. 如你所见，我选择。

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

这里滤波器的阶数(N)是10，临界频率(Wn)是1sqrt(60)。我的问题是，我怎样才能选择合适的滤波器阶数？我从N=1开始尝试，直到N=21，大于21的结果是 data_phase_unwrapped 都是NAN。我也试过，给 padlen 在 filtfilt 但我没有很好地理解它。

Answer 3

这是一个有点晚，但我发现这个答案，而翻译一些我自己的matlab代码。

TLDR：使用 mode="full" 中的任何一个函数的卷积。

我靠在 蝎子的食谱 来指导我完成这个任务。我的答案的其余部分实际上是该页的总结。Matlabs 的 fftfilt 函数可以用烹饪书中提到的任何一个卷积函数来代替 (np.convolve, scipy.signal.convolve, .oaconvolve, .fttconvolve)，如果你通过 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')

愉快的过滤！