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;

所以，我在这里尝试将 MATLAB 中定义的这个函数实现到 Python 中

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

您链接到的函数是一个 Python 等价于 Matlab 函数。它只是碰巧坏了。

无论如何，MNE 也有 fftfilt 函数使用的 overlap 和 add 方法的实现。它是库的私有函数，我不确定您是否可以将其称为与 Matlab 函数完全等效的函数，但也许它很有用。你可以在这里找到源代码：https://github.com/mne-tools/mne-python/blob/master/mne/filter.py#L41。

Answer 2

最后，我的代码得到了改进。我将fftfilt（应用两次）替换为scipy.signal.filtfilt（基本相同）。所以我翻译成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) 是 1/sqrt(60)。我的问题是，如何选择过滤器的适当顺序？我尝试从 N=1 到 N=21，大于 21 结果data_phase_unwrapped 都是 NAN。我也试过了，在filtfilt 中给padlen 赋值，但我不太理解。

Answer 3

这有点晚了，但我在翻译自己的一些 matlab 代码时找到了答案。

TLDR：将mode="full" 用于 scipy.signal 中的任何卷积函数

我依靠scipy's recipes 来指导我完成这件事。我的其余答案实际上是该页面的摘要。 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')

过滤愉快！

Answer 4

我在转换 MATLAB 代码时也遇到了问题。我来自这个 MATLAB 代码：

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

到这个python代码：

from scipy.signal import convolve

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

如果您想要比卷积更快的东西，请参阅 this 重叠并添加 fft 实现