0

我正在尝试将在 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 度偏移，如下所示 在此处输入图像描述

4

4 回答 4

0

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

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

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

快乐过滤！

于 2020-05-29T14:53:50.667 回答
0

您链接到的函数是与 Matlab 函数等效的 Python。它只是碰巧坏了。

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

于 2016-01-11T15:41:42.343 回答
0

我在转换 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)]

如果您想要比卷积更快的东西，请查看此重叠并添加 fft 实现

于 2021-03-08T17:24:06.393 回答
0

最后，我的代码得到了改进。我用（基本相同）替换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_freq11/秒。

结果是： 在此处输入图像描述

这基本上是对的。只有我对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。我也尝试过，为padlenin赋值，filtfilt但我不太理解。

于 2016-01-18T18:44:42.110 回答