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 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。

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 and 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：对 scipy.signal

中的任何卷积函数使用 mode="full"

我依靠 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 实现

Python 的 MATLAB fftfilt 等价物

MATLAB fftfilt equivalent for Python

python

matlab

filter

fft

fftpack